{ "cells": [ { "cell_type": "markdown", "id": "55e3d3d0-a548-4187-a422-4f8bdc1a2464", "metadata": {}, "source": [ "This notebook illustrates the use of Pari/GP in SageMath." ] }, { "cell_type": "markdown", "id": "72baa1c5-d3a7-469f-a21d-aec11d553329", "metadata": {}, "source": [ "# 1. Acceleration of Alternating Series" ] }, { "cell_type": "markdown", "id": "6c74e166-efc1-4377-9946-7b3c52a38d4e", "metadata": {}, "source": [ "To approximate $\\pi$ we could use\n", "\\begin{eqnarray*}\n", " \\frac{\\pi}{4} & = & \\sum_{k=0}^\\infty \\frac{(-1)^k}{2 k + 1} \\\\\n", " & = & 1 - \\frac{1}{3} + \\frac{1}{5} - \\frac{1}{7} + \\cdots\n", "\\end{eqnarray*}\n", "but this series converges very slowly.\n", "\n", "Because of slow convergence, we could conclude that $\\pi$\n", "cannot be computed as above, but `gp` proves us wrong." ] }, { "cell_type": "code", "execution_count": 1, "id": "c63ddaa0-769b-478d-ab08-9979bc9e7551", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1415926535897932384626433832795028842\n" ] } ], "source": [ "%%gp\n", "Pi" ] }, { "cell_type": "markdown", "id": "8440cede-76ba-483a-9c8e-39671ccc3ea1", "metadata": {}, "source": [ "Note that `gp` is case sensitive: `Pi` is an approximation of $\\pi$ relative to the working precision, whereas `pi` is just a name." ] }, { "cell_type": "code", "execution_count": 2, "id": "adcc10ef-fdc6-49fe-96af-2a4a39679042", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " realprecision = 38 significant digits\n" ] } ], "source": [ "%%gp\n", "\\p" ] }, { "cell_type": "markdown", "id": "23321598-8308-4a3d-9ca7-a194a3ca8e65", "metadata": {}, "source": [ "We start verifying that the series converges very slowly. First we use 1000 terms." ] }, { "cell_type": "code", "execution_count": 3, "id": "f415645e-74dd-4be2-b7c4-a2b1d5d1c36b", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1425916543395430509011277372522045659\n" ] } ], "source": [ "%%gp\n", "sum4pi1000 = 4*sum(k=0,1000,(-1.0)^k/(2*k+1))" ] }, { "cell_type": "code", "execution_count": 4, "id": "62cbfc55-5e16-4738-97c9-5637b10d8315", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-0.00099900074974981243848435397270168167593\n" ] } ], "source": [ "%%gp\n", "Pi - sum4pi1000" ] }, { "cell_type": "markdown", "id": "f4211617-37ad-4b96-9a5c-db1c6d411e4c", "metadata": {}, "source": [ "Next we use 100,000 terms in the series." ] }, { "cell_type": "code", "execution_count": 5, "id": "13745f71-2eb2-4ea4-bf50-fcf7c54afb36", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1416026534897939884601433645294403920\n" ] } ], "source": [ "%%gp\n", "sum4pi100000 = 4*sum(k=0,100000,(-1.0)^k/(2*k+1))" ] }, { "cell_type": "code", "execution_count": 6, "id": "fea3d492-cc86-4936-9f1d-fd256de108f3", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-9.9999000007499974999812499375078179010 E-6\n" ] } ], "source": [ "%%gp\n", "Pi - sum4pi100000" ] }, { "cell_type": "markdown", "id": "12f32047-024b-4005-a8e8-fd768dd5643e", "metadata": {}, "source": [ "Now let us compute 10 million terms in the series." ] }, { "cell_type": "code", "execution_count": 7, "id": "550f262c-e831-4afe-836a-e16df42198aa", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1415927535897832384633933832545028686\n" ] } ], "source": [ "%%gp\n", "sum4pi10000000 = 4*sum(k=0,10000000,(-1.0)^k/(2*k+1))" ] }, { "cell_type": "code", "execution_count": 8, "id": "be871233-b117-4827-aec2-c628d7cb8ca3", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-9.9999990000000749999974999984397308196 E-8\n" ] } ], "source": [ "%%gp\n", "Pi - sum4pi10000000" ] }, { "cell_type": "markdown", "id": "67a6d0ff-38c6-47ec-89b3-3437a6db0069", "metadata": {}, "source": [ "For two extra decimal places of accuracy we need 100 times more terms." ] }, { "cell_type": "markdown", "id": "d821199b-ad90-46e3-8f23-c9b542c4145f", "metadata": {}, "source": [ "Let us increase the precision to 100 decimal places." ] }, { "cell_type": "code", "execution_count": 9, "id": "bc8faa53-c2ed-4414-9c10-ce58029e54be", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " realprecision = 115 significant digits (100 digits displayed)\n" ] } ], "source": [ "%%gp\n", "\\p 100" ] }, { "cell_type": "code", "execution_count": 10, "id": "6eb23c82-a35e-4f4b-b9f7-66c94c27e540", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068\n" ] } ], "source": [ "%%gp\n", "pi100 = 4*sumalt(k=0, (-1)^k/(2*k+1) )" ] }, { "cell_type": "code", "execution_count": 11, "id": "aa0fbf7c-a64f-4837-992a-bdcde27e7f0d", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-3.045530204778779068 E-115\n" ] } ], "source": [ "%%gp\n", "Pi - pi100" ] }, { "cell_type": "markdown", "id": "3f84d1d5-7901-4782-8023-0f0ae840a097", "metadata": { "tags": [] }, "source": [ "`sumalt` accelerates the convergence of alternating series with\n", "algorithms of Cohen, Villegas, and Zagier." ] }, { "cell_type": "markdown", "id": "121595a2-c530-4dac-98dc-9d024ef52e01", "metadata": {}, "source": [ "# 2. Rational Approximations" ] }, { "cell_type": "markdown", "id": "d2b07436-b026-4de3-89a8-db430190cfad", "metadata": {}, "source": [ "We set the precision to 10 decimal places." ] }, { "cell_type": "code", "execution_count": 12, "id": "a571d167-f8ed-4315-808c-0a17f4710933", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " realprecision = 19 significant digits (10 digits displayed)\n" ] } ], "source": [ "%%gp\n", "\\p 10" ] }, { "cell_type": "code", "execution_count": 13, "id": "fbf7426b-503e-43e0-b1c2-2f2537cbfb30", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3/4\n" ] } ], "source": [ "%%gp\n", "a = 3/4" ] }, { "cell_type": "code", "execution_count": 14, "id": "9da328e8-9cf9-4744-b0c4-f41f6799dbc4", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\"t_FRAC\"\n" ] } ], "source": [ "%%gp\n", "type(a)" ] }, { "cell_type": "markdown", "id": "120650a3-6409-4fff-9a52-20264b2f4c78", "metadata": {}, "source": [ "The type conversion of `a` into a real number goes via the addition." ] }, { "cell_type": "code", "execution_count": 15, "id": "a6d45bd6-c6f5-4e75-bdf5-78124b5d4cad", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.7500000000\n" ] } ], "source": [ "%%gp\n", "b = a + 0.0" ] }, { "cell_type": "code", "execution_count": 16, "id": "04f9c581-17c8-4e91-9bbe-1e03d6299e2a", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\"t_REAL\"\n" ] } ], "source": [ "%%gp\n", "type(b)" ] }, { "cell_type": "markdown", "id": "84fcae3f-f47b-43b1-aeb3-d55fc21aa9bc", "metadata": {}, "source": [ "We can convert a real number to a fraction, via the `bestappr()` function." ] }, { "cell_type": "code", "execution_count": 17, "id": "d9afb4b0-62ec-4fc9-9549-0185eec8e388", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3/4\n" ] } ], "source": [ "%%gp\n", "bestappr(b)" ] }, { "cell_type": "markdown", "id": "53b7d721-4f09-41dc-a4ca-bfbdcbf4c015", "metadata": {}, "source": [ "The second argument to the `bestappr()` function is the upper bound on the size of denominator in the rational approximation of the first argument of `bestappr()`." ] }, { "cell_type": "code", "execution_count": 18, "id": "268b70cc-7d05-46aa-9087-979b5522ac81", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "22/7\n" ] } ], "source": [ "%%gp\n", "bestappr(Pi, 10)" ] }, { "cell_type": "code", "execution_count": 19, "id": "9a6f2e8a-6a16-4f17-934e-ec50d5afaf7e", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "22/7\n" ] } ], "source": [ "%%gp\n", "bestappr(Pi, 100)" ] }, { "cell_type": "code", "execution_count": 20, "id": "4b11aef9-cf7d-4fcb-81c0-d3a8c8f8495d", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "355/113\n" ] } ], "source": [ "%%gp\n", "bestappr(Pi, 1000)" ] }, { "cell_type": "code", "execution_count": 21, "id": "866194ce-6f48-4689-8f85-f79fd910d556", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "355/113\n" ] } ], "source": [ "%%gp\n", "bestappr(Pi, 1000)" ] }, { "cell_type": "markdown", "id": "0f10f2d4-6ae1-42f4-bbc6-1d2fecb60994", "metadata": {}, "source": [ "We can make a sequence of rational approximations for $\\sqrt{2}$, with larger and larger denominators." ] }, { "cell_type": "code", "execution_count": 22, "id": "1417e3f1-825e-4d9e-9265-4dbf7ebb5e67", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[7/5, 99/70, 1393/985, 8119/5741]\n" ] } ], "source": [ "%%gp\n", "apply(k->bestappr(sqrt(2),10^k),[1,2,3,4])" ] }, { "cell_type": "markdown", "id": "8c49d216-a5e3-4191-afe9-cbb054e776e5", "metadata": {}, "source": [ "The above construction is an example of an anonymous function applied to a list." ] }, { "cell_type": "markdown", "id": "8777769d-499a-41b4-91ba-a9bd308ddc36", "metadata": {}, "source": [ "Functions are defined quickly with `gp`." ] }, { "cell_type": "code", "execution_count": 23, "id": "721b902b-84d3-4d43-a14b-ba1c1556f02c", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(x)->x*(x^2-1)\n" ] } ], "source": [ "%%gp \n", "f(x) = x*(x^2-1)" ] }, { "cell_type": "code", "execution_count": 24, "id": "5d586a0e-a4f5-47f8-bbe7-82988508e61c", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6\n" ] } ], "source": [ "%%gp\n", "f(2)" ] }, { "cell_type": "markdown", "id": "7804dcf7-7b34-43f5-9cd4-3a3e184d241b", "metadata": {}, "source": [ "Evaluation of a function in a symbol returns a formula." ] }, { "cell_type": "code", "execution_count": 25, "id": "e4e06cc8-47a2-4c20-a8b8-41d7cf6e315b", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "z^3 - z\n" ] } ], "source": [ "%%gp\n", "p = f(z)" ] }, { "cell_type": "code", "execution_count": 26, "id": "a4b2ae13-55d4-4a08-b837-cc8b0efeec6d", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\"t_CLOSURE\"\n", "\"t_POL\"\n" ] } ], "source": [ "%%gp\n", "type(f)\n", "type(p)" ] }, { "cell_type": "markdown", "id": "a937d6c6-a4a5-44fa-96c0-32869dd05f44", "metadata": {}, "source": [ "We can take the derivative of a polynomial." ] }, { "cell_type": "code", "execution_count": 27, "id": "9da01642-c879-41ff-bd50-c077a3c37b57", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3*z^2 - 1\n" ] } ], "source": [ "%%gp\n", "dfz = deriv(f(z),z)" ] }, { "cell_type": "markdown", "id": "7fb0f3ec-ae2a-4276-a9e8-d5525f16dee6", "metadata": {}, "source": [ "And then use the expression of the derivative in a function." ] }, { "cell_type": "code", "execution_count": 28, "id": "d193dfff-1be4-4a2e-86ec-59afa94c9763", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(x)->substpol(dfz,z,x)\n" ] } ], "source": [ "%%gp\n", "g(x) = substpol(dfz,z,x)" ] }, { "cell_type": "code", "execution_count": 29, "id": "b144b3e3-a963-44a0-9a84-a5e27c421e19", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "11\n" ] } ], "source": [ "%%gp\n", "g(2)" ] }, { "cell_type": "markdown", "id": "11772fd9-8b1c-4194-b6fb-9e91379f931c", "metadata": {}, "source": [ "The above instructions are useful to compute complex roots using the Cauchy integral formula." ] }, { "cell_type": "markdown", "id": "6cb6a041-519d-4205-9f1d-be4008f9a59a", "metadata": {}, "source": [ "# 3. The Cauchy Integral Formula" ] }, { "cell_type": "markdown", "id": "90064cf8-f1a7-4ec1-8b47-ef23d6e835d6", "metadata": {}, "source": [ "The number of roots of $f(z)$ in a disk $C_{a,r}$ in the complex plane\n", "centered at $z=a$ and with radius $r$ is\n", "\n", "$$\n", " \\frac{1}{2 \\pi i} \\oint_{C_{a,r}} \\frac{f'(z)}{f(z)} dx\n", "$$\n", "\n", "We compute circular integrals with `intcirc`." ] }, { "cell_type": "markdown", "id": "4a5e491a-7649-44c0-bbd6-a2fab94957ce", "metadata": {}, "source": [ "This computation works in low precision. We use the functions `f` and `g` defined in the previous section." ] }, { "cell_type": "code", "execution_count": 30, "id": "24f0eb57-2cd9-4024-b47d-d0740edae221", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "%%gp\n", "\\p 5" ] }, { "cell_type": "code", "execution_count": 31, "id": "6bc25c3d-2577-432b-97f5-0e1939f3b739", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "z^3 - z\n" ] } ], "source": [ "%%gp\n", "f(z)" ] }, { "cell_type": "code", "execution_count": 32, "id": "82ebc8a1-ec29-4faa-bfa4-e52e6b2c9aac", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3*z^2 - 1\n" ] } ], "source": [ "%%gp\n", "g(z)" ] }, { "cell_type": "code", "execution_count": 33, "id": "15f4da5a-ac88-4c01-935f-81d8e09ffdb0", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0000 + 0.E-40*I\n" ] } ], "source": [ "%%gp\n", "intcirc(z=0,0.5,g(z)/f(z))" ] }, { "cell_type": "code", "execution_count": 34, "id": "ee6c2045-55d1-4f71-8203-38223841511a", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.0002 + 0.E-40*I\n" ] } ], "source": [ "%%gp\n", "intcirc(z=0,1.5,g(z)/f(z))" ] }, { "cell_type": "code", "execution_count": 35, "id": "aafebe82-9337-46f7-8656-684c046b7239", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0000 + 0.E-40*I\n" ] } ], "source": [ "%%gp\n", "intcirc(z=0.9,0.3,g(z)/f(z))" ] }, { "cell_type": "markdown", "id": "a69a81c1-0758-4e0a-aa09-06679f8bbda3", "metadata": {}, "source": [ "Changing the centers and radii of the circles, we can approximates the complex roots in a domain of interest." ] }, { "cell_type": "markdown", "id": "af572a81-bb5e-48f6-92ee-7cc0af46f775", "metadata": {}, "source": [ "# 4. Expansions and Continued Fractions" ] }, { "cell_type": "markdown", "id": "487ef3b7-d2ed-4375-8bda-6fe11c38e07c", "metadata": {}, "source": [ "Consider the binary expansion of 142." ] }, { "cell_type": "code", "execution_count": 36, "id": "b3d69f01-e5d8-45fb-adf1-793a8e6cde1a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 0, 0, 0, 1, 1, 1, 0]\n" ] } ], "source": [ "%%gp\n", "b = binary(142)" ] }, { "cell_type": "markdown", "id": "0fa0a082-14d7-4b76-a6aa-0f2e27feba20", "metadata": {}, "source": [ "The above coefficients can be computed as below." ] }, { "cell_type": "code", "execution_count": 37, "id": "68d7d9d1-5bd6-4108-9d90-a420ace6bf7b", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 + 2^2 + 2^3 + 2^7 + O(2^8)\n" ] } ], "source": [ "%%gp\n", "142 + O(2^8)" ] }, { "cell_type": "code", "execution_count": 38, "id": "194c2700-1292-4755-93f8-08c77eefb751", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 + 3*4 + 2*4^3 + O(4^6)\n" ] } ], "source": [ "%%gp\n", "c = 142 + O(4^6)" ] }, { "cell_type": "code", "execution_count": 39, "id": "d4c109a7-4b29-4b81-a08c-789cdd8ad05e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 + 3*4 + 2*4^3 + O(4^8)\n" ] } ], "source": [ "%%gp\n", "c = 142 + O(4^8)" ] }, { "cell_type": "code", "execution_count": 40, "id": "86e575cb-e2ad-40df-a7c3-9491b8533073", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "142\n" ] } ], "source": [ "%%gp\n", "d = truncate(c)" ] }, { "cell_type": "code", "execution_count": 41, "id": "ca5e5408-d5a0-4bef-83da-99dcaa6f3406", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 + 6*7 + 2*7^2 + O(7^6)\n" ] } ], "source": [ "%%gp \n", "e = 142 + O(7^6)" ] }, { "cell_type": "markdown", "id": "56409cf4-6601-4eb2-ad7a-6b47a8740341", "metadata": {}, "source": [ "The same syntax applies for power series:" ] }, { "cell_type": "code", "execution_count": 42, "id": "26c5e2b4-1958-4315-bee0-af9ec4d0fabd", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "z + z^2 + 1/3*z^3 - 1/30*z^5 + O(z^6)\n" ] } ], "source": [ "%%gp\n", "exp(z)*sin(z) + O(z^6)" ] }, { "cell_type": "markdown", "id": "4c80fd0f-a414-4c22-ae22-0fe621f77430", "metadata": {}, "source": [ "We can calculate with the number expansions." ] }, { "cell_type": "code", "execution_count": 43, "id": "b0aadb70-f92a-47dd-aeb4-b747a99271d0", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 + 6*7 + 2*7^2 + O(7^6)\n" ] } ], "source": [ "%%gp\n", "a = 142 + O(7^6)" ] }, { "cell_type": "code", "execution_count": 44, "id": "93c1624f-d773-494e-a999-d72b49bdb751", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6 + 6*7^2 + O(7^6)\n" ] } ], "source": [ "%%gp\n", "b = 300 + O(7^6)" ] }, { "cell_type": "code", "execution_count": 45, "id": "28ce07ad-291e-4ea9-b740-bc43901753d5", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 + 2*7^2 + 7^3 + O(7^6)\n" ] } ], "source": [ "%%gp\n", "c = a + b" ] }, { "cell_type": "code", "execution_count": 46, "id": "cdb37013-9708-43d1-9d0d-4e6f12c98555", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "442\n" ] } ], "source": [ "%%gp\n", "lift(c)" ] }, { "cell_type": "markdown", "id": "be09260e-a7d7-4cdb-95d8-21e982861435", "metadata": {}, "source": [ "The `lift(x)` lifts an element `x` of ${\\mathbb Z}_n$ to $\\mathbb Z$." ] }, { "cell_type": "markdown", "id": "ebff4483-8f96-4fcd-804d-80be92f7ae14", "metadata": {}, "source": [ "The code cells below illustrate the calculations with series." ] }, { "cell_type": "code", "execution_count": 47, "id": "77ad5945-433b-4933-8951-4eb09196ad98", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "z + z^2 + 1/3*z^3 - 1/30*z^5 + O(z^6)\n" ] } ], "source": [ "%%gp\n", "p = exp(z) * sin(z) + O(z^6)" ] }, { "cell_type": "code", "execution_count": 48, "id": "6b67c919-86fc-4d0f-a5d6-98ba808a39f3", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 + z - 1/3*z^3 - 1/6*z^4 - 1/30*z^5 + O(z^6)\n" ] } ], "source": [ "%%gp\n", "q = exp(z) * cos(z) + O(z^6)" ] }, { "cell_type": "code", "execution_count": 49, "id": "b8f52fc7-9ba6-482c-b183-42f03789c1e7", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 + 2*z + z^2 - 1/6*z^4 - 1/15*z^5 + O(z^6)\n" ] } ], "source": [ "%%gp\n", "r = p + q" ] }, { "cell_type": "code", "execution_count": 50, "id": "05ddad74-feaa-47b4-a7d3-785934a19b03", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-1/15*z^5 - 1/6*z^4 + z^2 + 2*z + 1\n" ] } ], "source": [ "%%gp\n", "rp = truncate(r)" ] }, { "cell_type": "markdown", "id": "95186d77-b43b-4b98-b489-9aa890012dfb", "metadata": {}, "source": [ "The `truncate(r)` converts a series `r` to a polynomial by truncation." ] }, { "cell_type": "code", "execution_count": 51, "id": "61a618d6-658f-4323-a8a6-85d5c2aaabec", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\"t_SER\"\n", "\"t_POL\"\n" ] } ], "source": [ "%%gp\n", "type(r)\n", "type(rp)" ] }, { "cell_type": "markdown", "id": "3f96ba93-2047-4f51-a9cf-2b1328578152", "metadata": { "tags": [] }, "source": [ "The continued fraction expansion of a rational number is\n", "computed by repeated greatest common divisor computations\n", "of numerator and denominator.\n", "\n", "$$\n", " \\sqrt{21} = 4 + \\frac{1}{1 + \\frac{1}{1 + \\frac{1}{2 \n", " + \\frac{1}{1 + \\frac{1}{1 + \\frac{1}{8 + \\cdots}}}}}}\n", "$$\n", "\n", "We denote the periodicity in $\\sqrt{21}$ by $[4, \\overline{1,1,2,1,1,8}]$." ] }, { "cell_type": "code", "execution_count": 52, "id": "829c6379-7618-4b23-97ec-1fb69f5b2c85", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[4, 1, 1, 2, 1, 1, 8, 1, 1, 2, 1, 1, 8, 1, 1, 2, 1, 1, 8]\n" ] } ], "source": [ "%%gp\n", "L = contfrac(sqrt(21), 19)" ] }, { "cell_type": "markdown", "id": "18622aa8-6d2e-4aa5-b911-c23fc48069c8", "metadata": {}, "source": [ "The $\\sqrt{2}$ has a very simple continued fraction." ] }, { "cell_type": "code", "execution_count": 53, "id": "f939cc2e-5be5-4ca3-a44c-49ea971a5d4c", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]\n" ] } ], "source": [ "%%gp\n", "contfrac(sqrt(2), 20)" ] }, { "cell_type": "markdown", "id": "0c7135fe-de59-4c02-b5d1-ee3d27d6735f", "metadata": {}, "source": [ "# 5. Integer Relation Detection" ] }, { "cell_type": "markdown", "id": "27f8f5d3-cd18-41cd-a6e0-512949f872f3", "metadata": {}, "source": [ "We have $\\log(15) = \\log(5 \\cdot 3) = \\log(5) + \\log(3)$." ] }, { "cell_type": "code", "execution_count": 54, "id": "4397405c-b392-4c6e-8172-049155b10ff2", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "%%gp\n", "\\p 10" ] }, { "cell_type": "code", "execution_count": 55, "id": "1b1a097e-e108-41fd-9a9d-3b6f5abdac81", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2.708050201, 1.609437912, 1.098612289]\n" ] } ], "source": [ "%%gp\n", "u = [log(15.0), log(5.0), log(3.0)]" ] }, { "cell_type": "code", "execution_count": 56, "id": "65576209-f44f-4175-8fdb-ca7d5e69fdee", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-1, 1, 1]~\n" ] } ], "source": [ "%%gp \n", "lindep(u)" ] }, { "cell_type": "markdown", "id": "3e199609-8d6b-4e23-82f1-bf9b2ea230d9", "metadata": {}, "source": [ "We find $(-1) \\log(15) + (+1) \\log(5) + (+1) \\log(3)$." ] }, { "cell_type": "markdown", "id": "ff882936-c6f2-4211-b58c-46deb1e9df1f", "metadata": { "tags": [] }, "source": [ "Calling `lindep(u,-3)` uses the PSLQ Algorithm." ] }, { "cell_type": "markdown", "id": "eab7f343-0533-4f84-a420-27374c3a1aee", "metadata": {}, "source": [ "Note: the type of is `t_COL` (column vector)." ] }, { "cell_type": "markdown", "id": "0a25624b-bc00-4d50-b4f4-f17a0a4800a1", "metadata": {}, "source": [ "# 6. Lattice Basis Reduction" ] }, { "cell_type": "markdown", "id": "0acc2a4e-8572-4036-9ba6-88a0520ec957", "metadata": {}, "source": [ "A lattice is similar to a vector space. Like a vector space,\n", "a lattice is spanned (or generated) by a finite number of basis\n", "vectors. Unlike a vector space, the basis vectors have only\n", "integer entries and only integer linear combinations of the\n", "basis vectors are allowed.\n", "\n", "Given any basis of a lattice, the LLL-algorithm finds a reduced\n", "basis, consisting of short vectors. A short vector has small length\n", "and its entries are small numbers, like the entries in the vectors of\n", "the subset sum problem are either zero or one.\n", "\n", "Consider the following weights: ${\\bf w} = (366,385,392,401,422,437)$\n", "and sum $s = 1215$. This leads to the lattice basis\n", "$(1,0,0,0,0,0,-366)$, $(0,1,0,0,0,0,-385)$, $(0,0,1,0,0,0,-392)$,\n", "$(0,0,0,1,0,0,-401)$, $(0,0,0,0,1,0,-422)$, $(0,0,0,0,0,1,-437)$,\n", "and $(0,0,0,0,0,0,1215)$. The output of the LLL-algorithm gives\n", "a reduced basis, which may contain a solution to\n", "the subset sum problem $s = {\\bf w} \\cdot {\\bf x}$." ] }, { "cell_type": "markdown", "id": "df7efccd-8395-4dd7-8900-cf7f585344ab", "metadata": {}, "source": [ "An application of Lattice Basis Reduction is a method to solve the knapsack problem." ] }, { "cell_type": "code", "execution_count": 57, "id": "1673899f-2bf9-4ce4-b012-480a24609682", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[366, 385, 392, 401, 422, 437]\n" ] } ], "source": [ "%%gp\n", "w = [366, 385, 392, 401, 422, 437]" ] }, { "cell_type": "code", "execution_count": 58, "id": "7af2a622-4365-4ab7-832c-bb474e48bffe", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\"t_VEC\"\n" ] } ], "source": [ "%%gp\n", "type(w)" ] }, { "cell_type": "code", "execution_count": 59, "id": "d4d264c2-568b-47bf-899a-30cbf10d129c", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "List([366, 385, 392, 401, 422, 437])\n" ] } ], "source": [ "%%gp\n", "L = List(w)" ] }, { "cell_type": "code", "execution_count": 60, "id": "77aee602-b4c0-4c06-b694-d858e8240e4a", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1215\n" ] } ], "source": [ "%%gp\n", "s = 1215" ] }, { "cell_type": "code", "execution_count": 61, "id": "bd49c4a6-fc35-43e5-b2cc-a72a9c4a2c28", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-1215\n" ] } ], "source": [ "%%gp\n", "listinsert(L,-s,1)" ] }, { "cell_type": "code", "execution_count": 62, "id": "d4902c9e-8d88-4c7c-8e75-d6e40956b890", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "List([-1215, 366, 385, 392, 401, 422, 437])\n" ] } ], "source": [ "%%gp\n", "L" ] }, { "cell_type": "code", "execution_count": 63, "id": "4f548e19-d7e3-432e-bbb1-363bbf8e1e17", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-1215, 366, 385, 392, 401, 422, 437]\n" ] } ], "source": [ "%%gp\n", "v = Vec(L)" ] }, { "cell_type": "code", "execution_count": 64, "id": "d0952495-f7af-4c4e-a2b4-2e8d5f599ffc", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 0, 0, 1, 1, 1, 0]~\n" ] } ], "source": [ "%%gp\n", "x = lindep(v)" ] }, { "cell_type": "markdown", "id": "35772e3c-c0fe-48c7-97d8-dbd2effb58b8", "metadata": {}, "source": [ "For a floating point number $x$, we can ask to compute\n", "the polynomial that has $x$ as its root." ] }, { "cell_type": "markdown", "id": "eacbb1dc-630b-47d4-aad8-59ca831d726d", "metadata": {}, "source": [ "`algdep`($x,k$) is a call to `lindep`([$1,x,\\ldots,x^k$])" ] }, { "cell_type": "code", "execution_count": 65, "id": "80cf9b98-3d12-4a7e-b295-bc56c563f5e9", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "%%gp\n", "\\p 10" ] }, { "cell_type": "code", "execution_count": 66, "id": "09f31ea0-2a11-465a-869d-09342808c21e", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.200936955\n" ] } ], "source": [ "%%gp\n", "x = 3.0^(1/6)" ] }, { "cell_type": "code", "execution_count": 67, "id": "95f52df6-9774-4589-8c08-b844405dd143", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x^8 - 3*x^2\n" ] } ], "source": [ "%%gp\n", "a = algdep(x,10)" ] }, { "cell_type": "code", "execution_count": 68, "id": "7dfa4804-a944-4f9b-9283-13d36b358447", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " *** at top-level: a=algdep(x,10,-3)\n", " *** ^---------------\n", " *** algdep: domain error in lindep2: accuracy < 0\n" ] } ], "source": [ "%%gp\n", "a = algdep(x,10,-3) " ] }, { "cell_type": "markdown", "id": "9dd7716e-d0e8-4592-9531-906921858c76", "metadata": {}, "source": [ "Using PSLQ with the flag `-3` does not work well." ] }, { "cell_type": "markdown", "id": "1fba768c-3f91-41ef-88cd-2a8a92da7e00", "metadata": {}, "source": [ "# 7. Linear Algebra" ] }, { "cell_type": "markdown", "id": "e918e725-82dc-4826-93d6-ee2ee373e76e", "metadata": {}, "source": [ "Let us make a linear system." ] }, { "cell_type": "code", "execution_count": 69, "id": "b5111dd7-773d-4101-8c65-1c43bb727dcd", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[ 5 -4]\n", "\n", "[-1 1]\n", "\n" ] } ], "source": [ "%%gp\n", "A = [5, -4; -1, 1]" ] }, { "cell_type": "code", "execution_count": 70, "id": "42bcef85-a31b-4cfb-9d48-bb8dccb26f01", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[ 80]\n", "\n", "[100]\n", "\n" ] } ], "source": [ "%%gp\n", "x = [80; 100]" ] }, { "cell_type": "code", "execution_count": 71, "id": "a6665ef4-4646-493b-944f-d63c953992a3", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[ 0]\n", "\n", "[20]\n", "\n" ] } ], "source": [ "%%gp\n", "b = A*x" ] }, { "cell_type": "markdown", "id": "ded90da4-695c-46f0-8420-7d823903cc03", "metadata": {}, "source": [ "The Hermite normal form is computed by `mathnf`:" ] }, { "cell_type": "code", "execution_count": 72, "id": "5aca190f-1fc9-4e0b-a6e1-2c534db80514", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1, 0; 0, 1], [1, 4; 1, 5]]\n" ] } ], "source": [ "%%gp\n", "L = mathnf(A,1)" ] }, { "cell_type": "code", "execution_count": 73, "id": "10dcf7d5-b563-455c-9915-695c25a7186f", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[1 4]\n", "\n", "[1 5]\n", "\n" ] } ], "source": [ "%%gp\n", "U = L[2]" ] }, { "cell_type": "code", "execution_count": 74, "id": "bd662051-d5c2-47ec-bf35-114bfc54c0cd", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[1 0]\n", "\n", "[0 1]\n", "\n" ] } ], "source": [ "%%gp\n", "U*A" ] }, { "cell_type": "markdown", "id": "4136cbb4-7b4d-4101-9d24-023d6ad8205f", "metadata": {}, "source": [ "To solve the linear system using the Hermite normal form,\n", "the right hand side vector is transformed." ] }, { "cell_type": "code", "execution_count": 75, "id": "be3e93f3-d3ab-4c90-a0bd-420127e7a208", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[ 80]\n", "\n", "[100]\n", "\n" ] } ], "source": [ "%%gp\n", "U*b" ] }, { "cell_type": "markdown", "id": "ffaa5730-9dab-417f-81ce-1c38eae871b2", "metadata": {}, "source": [ "Because `U*A` is the identity matrix, `U*b` equals `x`." ] }, { "cell_type": "code", "execution_count": 76, "id": "deab4b59-fbeb-433f-a6b0-6cb6e39dc3bb", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[ 80]\n", "\n", "[100]\n", "\n" ] } ], "source": [ "%%gp\n", "x" ] }, { "cell_type": "markdown", "id": "0e97da0c-99b3-4af2-b8a5-6e93b3f42859", "metadata": {}, "source": [ "To solve the linear system directly, use `matsolve`." ] }, { "cell_type": "code", "execution_count": 77, "id": "819e11c0-4c6e-4863-baa7-9cee998224fb", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[ 80]\n", "\n", "[100]\n", "\n" ] } ], "source": [ "%%gp\n", "y = matsolve(A,b)" ] }, { "cell_type": "markdown", "id": "ba84eef5-05db-4b73-859f-19dedde7a83f", "metadata": {}, "source": [ "# 8. Modular Arithmetic" ] }, { "cell_type": "markdown", "id": "d39f30f5-a298-4357-b965-a81a4b341871", "metadata": {}, "source": [ "Compute modulo a prime is one way to work with exact arithmetic and to avoid expression swell." ] }, { "cell_type": "code", "execution_count": 78, "id": "6d103c2c-76de-488b-8309-7db095d97603", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mod(4, 5)\n" ] } ], "source": [ "%%gp\n", "a = Mod(19,5)" ] }, { "cell_type": "code", "execution_count": 79, "id": "8aee0d0a-2c6e-4759-a2d6-fbce82eff4fe", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mod(3, 5)\n" ] } ], "source": [ "%%gp\n", "b = Mod(23,5)" ] }, { "cell_type": "code", "execution_count": 80, "id": "0ed9d4a7-1a62-4560-a141-38b1552ff99f", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mod(2, 5)\n" ] } ], "source": [ "%%gp\n", "a + b" ] }, { "cell_type": "code", "execution_count": 81, "id": "c0bf4514-1132-4bca-9f78-affdcb26f439", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2\n" ] } ], "source": [ "%%gp\n", "(19 + 23) % 5" ] }, { "cell_type": "code", "execution_count": 82, "id": "169bc06a-e82f-4248-b37e-f2362df2daba", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7\n" ] } ], "source": [ "%%gp\n", "19 % 5 + 23 % 5" ] }, { "cell_type": "code", "execution_count": 83, "id": "c4b2e219-18f9-4f75-96c8-54b0cee5308d", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\"t_INTMOD\"\n" ] } ], "source": [ "%%gp\n", "type(a)" ] }, { "cell_type": "code", "execution_count": 84, "id": "6d9c1712-504f-4dc3-b1fe-7ac70c018a59", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4\n" ] } ], "source": [ "%%gp\n", "c = lift(a)" ] }, { "cell_type": "code", "execution_count": 85, "id": "9165586d-e4ad-451f-8052-a4fb27a3d4d3", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\"t_INT\"\n" ] } ], "source": [ "%%gp\n", "type(c)" ] }, { "cell_type": "code", "execution_count": 86, "id": "ec566b54-0955-49d9-aa4c-fce117f0bc96", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mod(4, 5)\n" ] } ], "source": [ "%%gp\n", "inva = 1/a" ] }, { "cell_type": "code", "execution_count": 87, "id": "cd1a4c90-a398-499f-ad55-e4e2e20117ae", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mod(1, 5)\n" ] } ], "source": [ "%%gp\n", "a*inva" ] }, { "cell_type": "code", "execution_count": 88, "id": "d75f6a6f-cf05-41dc-8e23-ae61eb858c1b", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mod(2, 5)\n" ] } ], "source": [ "%%gp\n", "cba = a^(1/2)" ] }, { "cell_type": "code", "execution_count": 89, "id": "9f2ec96e-923a-4a78-8a85-43490e59aff5", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mod(4, 5)\n" ] } ], "source": [ "%%gp\n", "cba^2" ] }, { "cell_type": "markdown", "id": "f84a793a-6903-43c2-a0ce-6ec5dcb16d12", "metadata": {}, "source": [ "Consider the input data:\n", "\n", "* $M = [m_{ij}]$ is any integral matrix;\n", "\n", "* $D$ a column vector of moduli; and\n", "\n", "* $B$ an integer column vector.\n", "\n", "We look for $\\bf x$, a small integer solution to\n", "\n", "$$\n", " \\sum_{j} m_{ij} x_j \\equiv b_i \\quad \\mod d_i\n", "$$\n", "for $i$ running over all rows of $M$." ] }, { "cell_type": "code", "execution_count": 90, "id": "30c3b1d6-42d7-4edd-8f18-a5151695d6fa", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[1 1]\n", "\n", "[1 2]\n", "\n" ] } ], "source": [ "%%gp\n", "A = [1,1; 1,2]" ] }, { "cell_type": "code", "execution_count": 91, "id": "7d147c19-9489-431b-99c5-4ebde9a174a0", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[3, 7]~\n" ] } ], "source": [ "%%gp\n", "m = [3,7]~" ] }, { "cell_type": "code", "execution_count": 92, "id": "da42caa5-b955-4ec3-b76a-2c614c046de5", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0, 5]~\n" ] } ], "source": [ "%%gp\n", "y = [0,5]~" ] }, { "cell_type": "code", "execution_count": 93, "id": "81d78fbd-f0b8-4845-85ee-7a02f436f9c0", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[12, 0]~\n" ] } ], "source": [ "%%gp\n", "x = matsolvemod(A,m,y)" ] }, { "cell_type": "code", "execution_count": 94, "id": "72154340-56c2-4779-8ca2-fe82cb2f480e", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[12, 12]~\n" ] } ], "source": [ "%%gp\n", "A*x" ] }, { "cell_type": "markdown", "id": "23f4b9ef-603f-4b13-a076-7b187cc17452", "metadata": {}, "source": [ "The technique to solve linear systems of congruences is called \n", "Chinese remaindering and is attributed to Sun-tzu Suan-ching." ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.5", "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.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }