{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## The SOMOS 4 sequence\n",
    "\n",
    "The SOMOS 4 sequence is the sequence defined by the recurrence\n",
    "$$u_{n+4} = \\frac{u_{n+1}u_{n+3} + u_{n+2}^2}{u_n}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "def somos(u0, u1, u2, u3, n):\n",
    "    a, b, c, d = u0, u1, u2, u3\n",
    "    for _ in range(4, n+1):\n",
    "        a, b, c, d = b, c, d, (b*d + c*c) / a\n",
    "        print d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "...000000000000100\n",
      "...000000000001101\n",
      "...000000000110111\n",
      "...101010111010111\n",
      "...1100111101111\n",
      "...1110000010010\n",
      "...0001000111001\n",
      "...0000011111101\n",
      "...1000000110101\n",
      "...101101010011\n",
      "...110000000000\n",
      "...000101011101\n",
      "...001001101011\n",
      "...000011110011\n",
      "...11\n"
     ]
    }
   ],
   "source": [
    "Z2 = Zp(2,40,print_mode='digits')\n",
    "u0 = u1 = u2 = Z2(1,15); u3 = Z2(3,15)\n",
    "somos(u0,u1,u2,u3,18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "...000000000000100\n",
      "...000000000001101\n",
      "...000000000110111\n",
      "...101010111010111\n",
      "...101100111101111\n",
      "...1111110000010010\n",
      "...100001000111001\n",
      "...100000011111101\n",
      "...001000000110101\n",
      "...010101101010011\n",
      "...001110000000000\n",
      "...111000101011101\n",
      "...111001001101011\n",
      "...111000011110011\n",
      "...100000000000111\n"
     ]
    }
   ],
   "source": [
    "Z2 = ZpLC(2,40,print_mode='digits')\n",
    "u0 = u1 = u2 = Z2(1,15); u3 = Z2(3,15)\n",
    "somos(u0,u1,u2,u3,18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Basic arithmetic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2*3 + O(3^4), 2 + O(3^4))"
      ]
     },
     "execution_count": 4,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z3 = ZpLC(3)\n",
    "x = Z3(4, 6)\n",
    "y = Z3(2, 4)\n",
    "u = x + y\n",
    "v = x - y\n",
    "u, v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2 + 2*3 + O(3^6), 1 + 3 + O(3^4))"
      ]
     },
     "execution_count": 5,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(u+v, u-v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3 + 3^2 + O(3^7)"
      ]
     },
     "execution_count": 6,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x+x+x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ 81 648]\n",
       "[  0 729]"
      ]
     },
     "execution_count": 7,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = Z3.precision()\n",
    "L.precision_lattice([u,v])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 8,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L.number_of_diffused_digits([u,v])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Computing with Matrices\n",
    "#### Products of many matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0 0]\n",
       "[0 0]"
      ]
     },
     "execution_count": 9,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z2 = Zp(2,40,print_mode='digits')\n",
    "MS = MatrixSpace(Z2,2)\n",
    "M = MS(1)\n",
    "for _ in range(60):\n",
    "    M *= random_element(MS, prec=8)\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[...10101110000000000 ...00100110000000000]\n",
       "[...11111000000000000 ...11011000000000000]"
      ]
     },
     "execution_count": 10,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z2 = ZpLC(2,40,print_mode='digits')\n",
    "MS = MatrixSpace(Z2,2)\n",
    "M = MS(1)\n",
    "for _ in range(60):\n",
    "    M *= random_element(MS, prec=8)\n",
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Determinants and characteristic polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...1000000000"
      ]
     },
     "execution_count": 11,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z3 = Zp(3,40,print_mode='digits')\n",
    "D = diagonal_matrix([ Z3(1,5), Z3(9,5), Z3(27,5), Z3(81,5) ])\n",
    "D.determinant()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 12,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MS = D.parent()\n",
    "P = random_element(MS, prec=5)\n",
    "Q = random_element(MS, prec=5)\n",
    "M = P*D*Q\n",
    "M.determinant()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(...0000000000000000000000000000000000000001)*x^4 + (...00021)*x^3 + (...22200)*x^2 + (0)*x + (0)"
      ]
     },
     "execution_count": 13,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.characteristic_polynomial()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...2000000000"
      ]
     },
     "execution_count": 14,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z3 = ZpLC(3,40,print_mode='digits')\n",
    "D = diagonal_matrix([ Z3(1,5), Z3(9,5), Z3(27,5), Z3(81,5) ])\n",
    "MS = D.parent()\n",
    "P = random_element(MS, prec=5)\n",
    "Q = random_element(MS, prec=5)\n",
    "M = P*D*Q\n",
    "M.determinant()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...0000000000000000000000000000000000000001*x^4 + ...11111*x^3 + ...00100*x^2 + ...1100000*x + ...2000000000"
      ]
     },
     "execution_count": 15,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.characteristic_polynomial()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Computing with polynomials\n",
    "#### Euclidean algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "def euclidean(A,B):\n",
    "    while B != 0:\n",
    "        A, B = B, A % B\n",
    "    return A.monic()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(...0000000000000000000000000000000000000001)*x^5 + (...11000010)*x^4 + (...10100011)*x^3 + (...10010111)*x^2 + (...11111011)*x + (...11110100)"
      ]
     },
     "execution_count": 17,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2 = Qp(2,40,print_mode='digits')\n",
    "S.<x> = PolynomialRing(Q2)\n",
    "P = random_element(S, degree=10, prec=5, integral=True)\n",
    "Q = random_element(S, degree=10, prec=5, integral=True)\n",
    "D = x^5 + random_element(S, degree=4, prec=8, integral=True)\n",
    "D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "With high probability, $P$ and $Q$ are coprime, implying that the gcd of $DP$ is $DQ$ is $D$.\n",
    "However, we observe that we do not always get this expected result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(...1)*x^11 + (...10)*x^10 + (...1)*x^9 + (...1)*x^8 + (0)*x^7 + (...1000)*x^6 + (...1)*x^5 + (...1)*x^4 + (...1)*x^3 + (...100)*x^2 + (...1000)*x + (...100000)"
      ]
     },
     "execution_count": 18,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "euclidean(D*P, D*Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Here, we stopped too early since the remainder lost its precision.  If we use floating point arithmetic, the opposite problem occurs, since the remainder is not zero soon enough."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...000000000000001*x^5 + ...0000000001010010*x^4 + ...000000011101011*x^3 + ...0000000010011010*x^2 + ...00000000011001100*x + ...0000000000011110"
      ]
     },
     "execution_count": 19,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2 = QpFP(2,15,print_mode='digits')\n",
    "S.<x> = PolynomialRing(Q2)\n",
    "P = random_element(S, degree=10, prec=5, integral=True)\n",
    "Q = random_element(S, degree=10, prec=5, integral=True)\n",
    "D = x^5 + random_element(S, degree=4, prec=8, integral=True)\n",
    "D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...000000000000001"
      ]
     },
     "execution_count": 20,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "euclidean(D*P, D*Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "With lattice precision, we get the right answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...0000000000000000000000000000000000000001*x^5 + ...00010000*x^4 + ...10101011*x^3 + ...11100111*x^2 + ...10100110*x + ...11100100"
      ]
     },
     "execution_count": 21,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2 = QpLC(2,40,print_mode='digits')\n",
    "S.<x> = PolynomialRing(Q2)\n",
    "P = random_element(S, degree=10, prec=5, integral=True)\n",
    "Q = random_element(S, degree=10, prec=5, integral=True)\n",
    "D = x^5 + random_element(S, degree=4, prec=8, integral=True)\n",
    "D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...0000000000000000000000000000000000000001*x^5 + ...00010000*x^4 + ...10101011*x^3 + ...11100111*x^2 + ...10100110*x + ...11100100"
      ]
     },
     "execution_count": 22,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "euclidean(D*P, D*Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Grobner bases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x^3, x*y + ...1100010*x^2, y^2 + ...11001*x^2, z + ...0000000010*x]"
      ]
     },
     "execution_count": 23,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2 = Qp(2,print_mode='digits')\n",
    "R.<x,y,z> = PolynomialRing(Q2, order = 'invlex')\n",
    "F = [ Q2(2,10)*x + Q2(1,10)*z,\n",
    "      Q2(1,10)*x^2 + Q2(1,10)*y^2 - Q2(2,10)*z^2,\n",
    "      Q2(4,10)*y^2 + Q2(1,10)*y*z + Q2(8,10)*z^2 ]\n",
    "from sage.rings.polynomial.toy_buchberger import buchberger_improved\n",
    "g = buchberger_improved(ideal(F));\n",
    "g.sort()\n",
    "g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x^3, x*y + ...111100010*x^2, y^2 + ...1111111001*x^2, z + ...0000000010*x]"
      ]
     },
     "execution_count": 24,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2 = QpLC(2,print_mode='digits')\n",
    "R.<x,y,z> = PolynomialRing(Q2, order = 'invlex')\n",
    "F = [ Q2(2,10)*x + Q2(1,10)*z,\n",
    "      Q2(1,10)*x^2 + Q2(1,10)*y^2 - Q2(2,10)*z^2,\n",
    "      Q2(4,10)*y^2 + Q2(1,10)*y*z + Q2(8,10)*z^2 ]\n",
    "from sage.rings.polynomial.toy_buchberger import buchberger_improved\n",
    "g = buchberger_improved(ideal(F));\n",
    "g.sort()\n",
    "g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "L = Q2.precision()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "227"
      ]
     },
     "execution_count": 26,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L.number_of_diffused_digits()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "77"
      ]
     },
     "execution_count": 27,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L.number_of_tracked_elements()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath with ZpL",
   "name": "sage_ZpL"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}