fruit_math_solver.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elliptic Curve defined by y^2 = x^3 + 109*x^2 + 224*x over Rational Field\n"
     ]
    }
   ],
   "source": [
    "# From:\n",
    "# Bremner, A., & MacLeod, A.J. (2014). An unusual cubic representation problem.\n",
    "# http://publikacio.uni-eszterhazy.hu/2858/1/AMI_43_from29to41.pdf\n",
    "\n",
    "# Interested in solutions to (1.1):\n",
    "#   a       b       c  \n",
    "# ----- + ----- + ----- = n\n",
    "# b + c   a + c   a + b\n",
    "\n",
    "n = 4\n",
    "assert (n % 2 == 0)\n",
    "\n",
    "# Homogeneity gives the following cubic \n",
    "# C_N : N(a + b)(b + c)(c + a) = a(a + b)(c + a) + b(b + c)(a + b) + c(c + a)(b + c)\n",
    "# in 2 dimensional projective space\n",
    "\n",
    "# And we get the following elliptic\n",
    "# E_N : y^2 = x^3 + (4N^2 + 12N − 3)x^2 + 32(N + 3)x\n",
    "\n",
    "E = EllipticCurve([0, ((4 * n ** 2) + (12 * n) - 3), 0, (32 * (n + 3)), 0]); print(E)\n",
    "\n",
    "# The accompanying morphisms\n",
    "# Projective cubic to elliptic \n",
    "def cubic_to_elliptic(a, b, c):\n",
    "    x = (-4 * (a + b + 2*c) * (n + 3))/((2 * a + 2 * b - c) + n * (a + b))\n",
    "    y = (4 * (a - b) * (n + 3) * (2 * n + 5))/((2 * a + 2 * b - c) + n * (a + b))\n",
    "    return x, y\n",
    "\n",
    "# The inverse\n",
    "def elliptic_to_cubic(x, y):\n",
    "    a = (8 * (n + 3) - x + y)/(2 * (4 - x) * (n + 3))\n",
    "    b = (8 * (n + 3) - x - y)/(2 * (4 - x) * (n + 3))\n",
    "    c = (-4 * (n + 3) - x * (n + 2))/((4 - x) * (n + 3))\n",
    "    return a, b, c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-100 : 260 : 1)\n",
      "Torsion Subgroup isomorphic to Z/6 associated to the Elliptic Curve defined by y^2 = x^3 + 109*x^2 + 224*x over Rational Field\n",
      "(0 : 1 : 0)\n",
      "(56 : 728 : 1)\n",
      "(4 : 52 : 1)\n",
      "(0 : 0 : 1)\n",
      "(4 : -52 : 1)\n",
      "(56 : -728 : 1)\n"
     ]
    }
   ],
   "source": [
    "# Get generator point for finding solutions\n",
    "try:\n",
    "    Q = E.gen(0); print(Q)\n",
    "except:\n",
    "    # Just in case we have a bad n \n",
    "    print(f\"Error!\\nNo solution exists for n = {n}\")\n",
    "    \n",
    "# Torsion subgroup always isomorphic to Z/Z6 by Lemma 2.1\n",
    "G = E.torsion_subgroup(); print(G);\n",
    "\n",
    "for t in G:\n",
    "    print(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Condition on the x coordinate of the elliptic for a \n",
    "# solution to be positive in the cubic, by Theorem 4.1\n",
    "\n",
    "# There is a minor error in the figure (4.1), the numerator \n",
    "# of the lower bound is missing a parenthesis on the right\n",
    "\n",
    "lbound1 = (3 - 12 * n - 4 * n ** 2 - ((2 * n + 5) * sqrt(4 * n ** 2 + 4 * n - 15))) / 2\n",
    "ubound1 = -2 * (n + 3) * (n + sqrt(n ** 2 - 4))\n",
    "\n",
    "lbound2 = -2 * (n + 3) * (n - sqrt(n ** 2 - 4))\n",
    "ubound2 = -4 * ((n + 3) / (n + 2))\n",
    "\n",
    "def positive_solution(pt):\n",
    "    x = pt[0]          \n",
    "    if (x > lbound1 and x < ubound1):\n",
    "        return True\n",
    "    if (x > lbound2 and x < ubound2):\n",
    "        return True\n",
    "    return False\n",
    "\n",
    "\n",
    "def find_solution():\n",
    "    temp = Q\n",
    "    m = 0;\n",
    "    while True:\n",
    "        # Take our current point, iQ, try torsion\n",
    "        for t in G:\n",
    "            temp_tor = temp + E(t)\n",
    "\n",
    "            if positive_solution(temp_tor):\n",
    "                return temp_tor, m, t\n",
    "\n",
    "        m += 1\n",
    "        temp += Q\n",
    "\n",
    "ell_solution, m, t = find_solution()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found solution with max length: 81\n",
      "Solution = 8Q + (0 : 1 : 0) \n",
      "a=154476802108746166441951315019919837485664325669565431700026634898253202035277999\n",
      "b=36875131794129999827197811565225474825492979968971970996283137471637224634055579\n",
      "c=4373612677928697257861252602371390152816537558161613618621437993378423467772036\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Get (projective) cubic points\n",
    "solution = elliptic_to_cubic(ell_solution[0], ell_solution[1])\n",
    "\n",
    "# LCM of denominators to convert to integers\n",
    "scale_factor = lcm([solution[0].denom(), solution[1].denom(), solution[2].denom()])\n",
    "\n",
    "# Integer solution\n",
    "a, b, c = solution[0] * scale_factor, solution[1] * scale_factor, solution[2] * scale_factor\n",
    "solution = [a, b, c]\n",
    "\n",
    "# Sanity\n",
    "assert(all(int(x) == x for x in solution))\n",
    "assert(all(x > 0 for x in solution))\n",
    "\n",
    "sol_len = max(len(str(a)), len(str(b)), len(str(c)))\n",
    "\n",
    "print('Found solution with max length:', sol_len)\n",
    "\n",
    "# Compare to Table 2\n",
    "print(f'Solution = {m}Q + {t} ')\n",
    "\n",
    "# Printing the solutions\n",
    "print(f'a={a}\\nb={b}\\nc={c}\\n')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.1",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}