Fruit Math

fruitmath.py

#!/usr/bin/env python3

'''
"An unusual cubic representation problem" Solver

Fruitmath  Copyright (C) 2019  Ryoichi Tanaka

This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.

'''

'''

It Solves equation:
    a/(b+c) + b(a+c) + c/(a+b) = N

References:

[1]  "An unusual cubic representation problem",  Andrew Bremnera , Allan Macleodb
     http://ami.ektf.hu/uploads/papers/finalpdf/AMI_43_from29to41.pdf

[2]  https://www.quora.com/How-do-you-find-the-positive-integer-solutions-to-frac-x-y+z-+-frac-y-z+x-+-frac-z-x+y-4
     Alon Amit 

'''



import fractions
import numpy as np
import sys

# lazy me, type less 
fr = fractions.Fraction

# stupid brute force solver
def demo1(N, M):
  for a in range(1, M-2):
    for b in range(a + 1, M-1):
      for c in range(b + 1, M):
        if a+b == 0 or b+c == 0 or a+c == 0:
          pass
        else:
          s = fr(a, b+c) + fr(b, a+c) + fr(c, a+b)
          #print (a, b, c, s)
          if N == s:
            print ("bingo!", a, b, c)


class Curve(object):
  def __init__(s, N):
    if N < 2:
      raise ValueError('N must >= 2')
    s.N = N
    s.A = (4*s.N**2 + 12*s.N - 3)
    s.B = 32*(s.N + 3)
    s.x = 0
    s.y = 0
    
    print("y^2 = x^3 + {}x^2 + {}x".format(s.A, s.B))

    # the "egg" part of curve
    # solve x^3 + Ax^2 + Bx > 0, one of  zero cross points is 0, 
    # the other 2 are roots of:
    #  x^3 + Ax^2 + Bx
    s.PosL = 0
    s.PosR = 0
    det = s.A * s.A - 4 * s.B
    if det > 0:
      d, rem = isqrt(det)
      x1 = int((-s.A - d) / 2)
      x2 = int((-s.A + d) / 2)
      print (" positive range ({},{})".format(x1, x2))
      s.PosL = x1 
      if (rem != 0):
        s.PosR = x2 + 2
      else:
        s.PosR = x2 + 1

  def xy2abc(s):
    if 0 == (4 - s.x) * (s.N + 3):
      return -1, -1, -1
    a = fr( 8*(s.N + 3) - s.x + s.y       , 2 * (4 - s.x) * (s.N + 3))
    b = fr( 8*(s.N + 3) - s.x - s.y       , 2 * (4 - s.x) * (s.N + 3))
    c = fr(-4*(s.N + 3) - (s.N + 2) * s.x ,     (4 - s.x) * (s.N + 3))
    return a,b,c

  def setxy(s, x, y):
    if y * y == x * x * x + s.A * x * x + s.B * x:
      s.x = x
      s.y = y
    else:
      print('Point({},{}) not on curve!'.format(x,y))
      raise ValueError('Point not on curve!')
  
  def y2fromx(s, x):
    y2 = x * x * x + s.A * x * x + s.B * x
    return y2

  def pd(s, x1, y1): # point double
    L = fr(3 * x1 * x1 + 2 * s.A * x1 + s.B, 2 * y1)
    xr = L * L - s.A - 2 * x1
    yr = L * (x1 - xr) - y1
    # print ("double ", x1, y1, L, xr, yr)
    return xr, yr
  
  def pa(s, x1, y1, x2, y2): # point add
    if (x1 == x2): # and (y1 == y2):
      return s.pd(x1, y1)
    L = fr(y2 - y1, x2 - x1)
    x = L * L - s.A - (x1 + x2)
    y = L * (x1 - x) - y1
    return x,y

# Newton's method works perfectly well on integers:
def isqrt(n):
  x = n
  y = (x + 1) // 2
  while y < x:
    x = y
    y = (x + n // x) // 2
  return x, (n - x * x)

def findfrompoint(cc, x, y, NC):
  # print("G", x, y) 
  try:
    cc.setxy(x,y)
  except:
    return
  i = 1
  while i < NC:
    a, b, c = cc.xy2abc()
    # print("round [{}] a = {}, b = {}, c = {} ".format(i, a, b, c))
    if (a > 0) and (b > 0) and (c > 0):
      lm = np.lcm.reduce([d.denominator for d in [a, b, c]])
      r = [(n * lm).numerator for n in [a , b, c]]
      print("bingo after {} rounds, x={}, y={}\na:{}\nb:{}\nc:{}".format(i, x, y, r[0], r[1], r[2]))
      break
    i += 1
    xr, yr = cc.pa(cc.x, cc.y,  x, y)
    cc.setxy(xr,yr)

def main():
  N = 4
  
  if len(sys.argv) > 1:
    N = int(sys.argv[1])
  
  if len(sys.argv) > 2:
    NC = int(sys.argv[2])
  cc = Curve(N)

  print("Finding integer points for x between {}, {}".format(max(cc.PosL, -10000000), cc.PosR))
  for x in range (max(cc.PosL, -10000000), cc.PosR):
    y2 = cc.y2fromx(x)
    if y2 > 0:
      y, rem = isqrt(y2)
      if 0 == rem and y > 0:
        print("x y2 y rem", x, y2, y, rem) 
        findfrompoint(cc, x, y, NC)

if __name__ == "__main__":
  main()