from math import *

syn_month =  29.530589    # synodic month / d
sid_year  = 365.25636042  # reference sidereal year / d
jul_year  = 365.25        # julian year / d
red_ratio =   3.0
#red_ratio =   6.0
zmin      =  12
zmax      =  90

year      = jul_year
#year      = sid_year

#ratio0    = year / syn_month / (red_ratio - 1)
ratio0 = 6.1843

def clamp(x):
  return min(max(x, zmin), zmax)

solutions = []
aerrmin = inf

for z2 in range(clamp(floor(ratio0 * zmin**2 / zmax)), clamp(ceil(ratio0 * zmax**2 / zmin)) + 1):
  for z1 in range(clamp(floor(z2 * zmin / (ratio0 * zmax))), clamp(ceil(z2 * zmax / (ratio0 * zmin))) + 1):
    for z3 in range(clamp(floor(z2 * zmin / (ratio0 * z1))), min(clamp(ceil(z2 * zmax / (ratio0 * z1))), z1) + 1):
      z4 = clamp(round(ratio0 * z1 * z3 / z2))
      if z4 > z2:
        break
      ratio = z2 * z4 / (z1 * z3)
      error = ratio / ratio0 - 1
      aerr = abs(error)
      if aerr <= aerrmin:
        if aerr < aerrmin:
          solutions = []
          aerrmin = aerr
        solutions.append((z1, z2, z3, z4, ratio, error))

if solutions:
  _, _, _, _, ratio, error = solutions[0]
  print(f'{ratio0 = :10.8f}')
  print(f'{ratio  = :10.8f}   {error = :+.2e}')
  print()
  print(f'  z1   z2   z3   z4')
  print(f'-------------------')
  for z1, z2, z3, z4, ratio, error in solutions:
    print(f'{z1:4d} {z2:4d} {z3:4d} {z4:4d}')
  print(f'-------------------')