docs/Newton-Raphson.ipynb

   1 {
   2  "cells": [
   3   {
   4    "cell_type": "markdown",
   5    "metadata": {},
   6    "source": [
   7     "# [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)\n",
   8     "Let's use the Newton-Raphson method for finding the root of an equation to write a function that can compute the square root of a number.\n",
   9     "\n",
  10     "Cf. [\"Why Functional Programming Matters\" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf)"
  11    ]
  12   },
  13   {
  14    "cell_type": "markdown",
  15    "metadata": {},
  16    "source": [
  17     "## A Generator for Approximations\n",
  18     "\n",
  19     "To make a generator that generates successive approximations let’s start by assuming an initial approximation and then derive the function that computes the next approximation:\n",
  20     "\n",
  21     "       a F\n",
  22     "    ---------\n",
  23     "        a'"
  24    ]
  25   },
  26   {
  27    "cell_type": "markdown",
  28    "metadata": {},
  29    "source": [
  30     "### A Function to Compute the Next Approximation\n",
  31     "\n",
  32     "This is the equation for computing the next approximate value of the square root:\n",
  33     "\n",
  34     "$a_{i+1} = \\frac{(a_i+\\frac{n}{a_i})}{2}$"
  35    ]
  36   },
  37   {
  38    "cell_type": "markdown",
  39    "metadata": {},
  40    "source": [
  41     "Starting with $\\frac{(a_i+\\frac{n}{a_i})}{2}$ we can derive the Joy expression to compute it using abstract dummy variables to stand in for actual values.  First undivide by two:\n",
  42     "\n",
  43     "$(a_i+\\frac{n}{a_i})$ `2 /`"
  44    ]
  45   },
  46   {
  47    "cell_type": "markdown",
  48    "metadata": {},
  49    "source": [
  50     "Then unadd terms:\n",
  51     "\n",
  52     "$a_i$ $\\frac{n}{a_i}$ `+ 2 /`"
  53    ]
  54   },
  55   {
  56    "cell_type": "markdown",
  57    "metadata": {},
  58    "source": [
  59     "Undivide again:\n",
  60     "\n",
  61     "$a_i$ $n$ $a_i$ `/ + 2 /`"
  62    ]
  63   },
  64   {
  65    "cell_type": "markdown",
  66    "metadata": {},
  67    "source": [
  68     "Finally deduplicate the $a_i$ term:\n",
  69     "\n",
  70     "$a_i$ $n$ `over / + 2 /`"
  71    ]
  72   },
  73   {
  74    "cell_type": "markdown",
  75    "metadata": {},
  76    "source": [
  77     "Let's try out this function `over / + 2 /` on an example:\n",
  78     "\n",
  79     "    F == over / + 2 /"
  80    ]
  81   },
  82   {
  83    "cell_type": "code",
  84    "execution_count": 1,
  85    "metadata": {},
  86    "outputs": [
  87     {
  88      "name": "stdout",
  89      "output_type": "stream",
  90      "text": []
  91     }
  92    ],
  93    "source": [
  94     "[F over / + 2 /] inscribe"
  95    ]
  96   },
  97   {
  98    "cell_type": "markdown",
  99    "metadata": {},
 100    "source": [
 101     "In order to use this function `F` we have to provide an initial estimate for the value of the square root, and we want to keep the input value `n` handy for iterations (we don't want the user to have to keep reentering it.)"
 102    ]
 103   },
 104   {
 105    "cell_type": "code",
 106    "execution_count": 2,
 107    "metadata": {},
 108    "outputs": [
 109     {
 110      "name": "stdout",
 111      "output_type": "stream",
 112      "text": [
 113       "  5 36 • F\n",
 114       "  5 36 • over / + 2 /\n",
 115       "5 36 5 • / + 2 /\n",
 116       "   5 7 • + 2 /\n",
 117       "    12 • 2 /\n",
 118       "  12 2 • /\n",
 119       "     6 • \n",
 120       "\n",
 121       "6"
 122      ]
 123     }
 124    ],
 125    "source": [
 126     "clear\n",
 127     "\n",
 128     "5 36 [F] trace"
 129    ]
 130   },
 131   {
 132    "cell_type": "markdown",
 133    "metadata": {},
 134    "source": [
 135     "The initial estimate can be 2, and we can `cons` the input value onto a quote with `F`:"
 136    ]
 137   },
 138   {
 139    "cell_type": "code",
 140    "execution_count": 3,
 141    "metadata": {},
 142    "outputs": [
 143     {
 144      "name": "stdout",
 145      "output_type": "stream",
 146      "text": [
 147       "6"
 148      ]
 149     }
 150    ],
 151    "source": [
 152     "[F1 2 swap [F] cons] inscribe"
 153    ]
 154   },
 155   {
 156    "cell_type": "code",
 157    "execution_count": 4,
 158    "metadata": {},
 159    "outputs": [
 160     {
 161      "name": "stdout",
 162      "output_type": "stream",
 163      "text": [
 164       "      36 • F1\n",
 165       "      36 • 2 swap [F] cons\n",
 166       "    36 2 • swap [F] cons\n",
 167       "    2 36 • [F] cons\n",
 168       "2 36 [F] • cons\n",
 169       "2 [36 F] • \n",
 170       "\n",
 171       "2 [36 F]"
 172      ]
 173     }
 174    ],
 175    "source": [
 176     "clear\n",
 177     "\n",
 178     "36 [F1] trace"
 179    ]
 180   },
 181   {
 182    "cell_type": "code",
 183    "execution_count": 5,
 184    "metadata": {},
 185    "outputs": [
 186     {
 187      "name": "stdout",
 188      "output_type": "stream",
 189      "text": [
 190       "     2 • 36 F\n",
 191       "  2 36 • F\n",
 192       "  2 36 • over / + 2 /\n",
 193       "2 36 2 • / + 2 /\n",
 194       "  2 18 • + 2 /\n",
 195       "    20 • 2 /\n",
 196       "  20 2 • /\n",
 197       "    10 • \n",
 198       "\n",
 199       "10"
 200      ]
 201     }
 202    ],
 203    "source": [
 204     "trace"
 205    ]
 206   },
 207   {
 208    "cell_type": "code",
 209    "execution_count": 6,
 210    "metadata": {},
 211    "outputs": [
 212     {
 213      "name": "stdout",
 214      "output_type": "stream",
 215      "text": [
 216       "6"
 217      ]
 218     }
 219    ],
 220    "source": [
 221     "36 F"
 222    ]
 223   },
 224   {
 225    "cell_type": "code",
 226    "execution_count": 7,
 227    "metadata": {},
 228    "outputs": [
 229     {
 230      "name": "stdout",
 231      "output_type": "stream",
 232      "text": []
 233     }
 234    ],
 235    "source": [
 236     "clear"
 237    ]
 238   },
 239   {
 240    "cell_type": "code",
 241    "execution_count": 8,
 242    "metadata": {},
 243    "outputs": [
 244     {
 245      "name": "stdout",
 246      "output_type": "stream",
 247      "text": [
 248       "[2 [36 F] codireco]"
 249      ]
 250     }
 251    ],
 252    "source": [
 253     "36 F1 make_generator"
 254    ]
 255   },
 256   {
 257    "cell_type": "code",
 258    "execution_count": 9,
 259    "metadata": {},
 260    "outputs": [
 261     {
 262      "name": "stdout",
 263      "output_type": "stream",
 264      "text": [
 265       "6"
 266      ]
 267     }
 268    ],
 269    "source": [
 270     "x x x first"
 271    ]
 272   },
 273   {
 274    "cell_type": "code",
 275    "execution_count": 10,
 276    "metadata": {},
 277    "outputs": [
 278     {
 279      "name": "stdout",
 280      "output_type": "stream",
 281      "text": [
 282       "6 12"
 283      ]
 284     }
 285    ],
 286    "source": [
 287     "144 F1 make_generator x x x x first"
 288    ]
 289   },
 290   {
 291    "cell_type": "code",
 292    "execution_count": 11,
 293    "metadata": {},
 294    "outputs": [
 295     {
 296      "name": "stdout",
 297      "output_type": "stream",
 298      "text": [
 299       "2 [36 F]"
 300      ]
 301     }
 302    ],
 303    "source": [
 304     "clear\n",
 305     "\n",
 306     "2 [36 F]"
 307    ]
 308   },
 309   {
 310    "cell_type": "code",
 311    "execution_count": 12,
 312    "metadata": {
 313     "scrolled": true
 314    },
 315    "outputs": [
 316     {
 317      "name": "stdout",
 318      "output_type": "stream",
 319      "text": [
 320       "2 [36 F] false"
 321      ]
 322     }
 323    ],
 324    "source": [
 325     "[first] [pop sqr] fork - abs 3 <"
 326    ]
 327   },
 328   {
 329    "cell_type": "code",
 330    "execution_count": 13,
 331    "metadata": {},
 332    "outputs": [
 333     {
 334      "name": "stdout",
 335      "output_type": "stream",
 336      "text": [
 337       "10"
 338      ]
 339     }
 340    ],
 341    "source": [
 342     "pop i"
 343    ]
 344   },
 345   {
 346    "cell_type": "code",
 347    "execution_count": 14,
 348    "metadata": {
 349     "scrolled": true
 350    },
 351    "outputs": [
 352     {
 353      "name": "stdout",
 354      "output_type": "stream",
 355      "text": [
 356       "10 [36 F] false"
 357      ]
 358     }
 359    ],
 360    "source": [
 361     "[36 F] [first] [pop sqr] fork - abs 3 <"
 362    ]
 363   },
 364   {
 365    "cell_type": "code",
 366    "execution_count": 15,
 367    "metadata": {},
 368    "outputs": [
 369     {
 370      "name": "stdout",
 371      "output_type": "stream",
 372      "text": [
 373       "6"
 374      ]
 375     }
 376    ],
 377    "source": [
 378     "pop i"
 379    ]
 380   },
 381   {
 382    "cell_type": "code",
 383    "execution_count": 16,
 384    "metadata": {
 385     "scrolled": true
 386    },
 387    "outputs": [
 388     {
 389      "name": "stdout",
 390      "output_type": "stream",
 391      "text": [
 392       "6 [36 F] true"
 393      ]
 394     }
 395    ],
 396    "source": [
 397     "[36 F] [first] [pop sqr] fork - abs 3 <"
 398    ]
 399   },
 400   {
 401    "cell_type": "code",
 402    "execution_count": null,
 403    "metadata": {},
 404    "outputs": [],
 405    "source": []
 406   },
 407   {
 408    "cell_type": "code",
 409    "execution_count": 17,
 410    "metadata": {},
 411    "outputs": [
 412     {
 413      "name": "stdout",
 414      "output_type": "stream",
 415      "text": [
 416       "2"
 417      ]
 418     }
 419    ],
 420    "source": [
 421     "clear\n",
 422     "\n",
 423     "2"
 424    ]
 425   },
 426   {
 427    "cell_type": "code",
 428    "execution_count": 18,
 429    "metadata": {},
 430    "outputs": [
 431     {
 432      "name": "stdout",
 433      "output_type": "stream",
 434      "text": [
 435       "6"
 436      ]
 437     }
 438    ],
 439    "source": [
 440     "[] true [i [36 F] [first] [pop sqr] fork - abs 3 >] loop pop"
 441    ]
 442   },
 443   {
 444    "cell_type": "code",
 445    "execution_count": 19,
 446    "metadata": {},
 447    "outputs": [
 448     {
 449      "name": "stdout",
 450      "output_type": "stream",
 451      "text": [
 452       "12"
 453      ]
 454     }
 455    ],
 456    "source": [
 457     "clear\n",
 458     "\n",
 459     "7 [] true [i [144 F] [first] [pop sqr] fork - abs 3 >] loop pop"
 460    ]
 461   },
 462   {
 463    "cell_type": "code",
 464    "execution_count": 20,
 465    "metadata": {},
 466    "outputs": [
 467     {
 468      "name": "stdout",
 469      "output_type": "stream",
 470      "text": [
 471       "120"
 472      ]
 473     }
 474    ],
 475    "source": [
 476     "clear\n",
 477     "\n",
 478     "7 [] true [i [14400 F] [first] [pop sqr] fork - abs 3 >] loop pop"
 479    ]
 480   },
 481   {
 482    "cell_type": "markdown",
 483    "metadata": {},
 484    "source": [
 485     "broken due to no float div\n",
 486     "\n",
 487     "    clear\n",
 488     "\n",
 489     "    7 [] true [i [1000 F] [first] [pop sqr] fork - abs 10 >] loop pop"
 490    ]
 491   },
 492   {
 493    "cell_type": "markdown",
 494    "metadata": {},
 495    "source": [
 496     "### Make it into a Generator\n",
 497     "\n",
 498     "Our generator would be created by:\n",
 499     "\n",
 500     "    a [dup F] make_generator\n",
 501     "\n",
 502     "With n as part of the function F, but n is the input to the sqrt function we’re writing. If we let 1 be the initial approximation:\n",
 503     "\n",
 504     "    1 n 1 / + 2 /\n",
 505     "    1 n/1   + 2 /\n",
 506     "    1 n     + 2 /\n",
 507     "    n+1       2 /\n",
 508     "    (n+1)/2\n",
 509     "\n",
 510     "The generator can be written as:\n",
 511     "\n",
 512     "    23 1 swap  [over / + 2 /] cons [dup] swoncat make_generator\n",
 513     "    1 23       [over / + 2 /] cons [dup] swoncat make_generator\n",
 514     "    1       [23 over / + 2 /]      [dup] swoncat make_generator\n",
 515     "    1   [dup 23 over / + 2 /]                    make_generator"
 516    ]
 517   },
 518   {
 519    "cell_type": "code",
 520    "execution_count": null,
 521    "metadata": {
 522     "scrolled": true
 523    },
 524    "outputs": [],
 525    "source": [
 526     "define('gsra 1 swap [over / + 2 /] cons [dup] swoncat make_generator')"
 527    ]
 528   },
 529   {
 530    "cell_type": "code",
 531    "execution_count": null,
 532    "metadata": {},
 533    "outputs": [],
 534    "source": [
 535     "J('23 gsra')"
 536    ]
 537   },
 538   {
 539    "cell_type": "markdown",
 540    "metadata": {},
 541    "source": [
 542     "Let's drive the generator a few time (with the `x` combinator) and square the approximation to see how well it works..."
 543    ]
 544   },
 545   {
 546    "cell_type": "code",
 547    "execution_count": null,
 548    "metadata": {},
 549    "outputs": [],
 550    "source": [
 551     "J('23 gsra 6 [x popd] times first sqr')"
 552    ]
 553   },
 554   {
 555    "cell_type": "markdown",
 556    "metadata": {},
 557    "source": [
 558     "## Finding Consecutive Approximations within a Tolerance\n",
 559     "\n",
 560     "From [\"Why Functional Programming Matters\" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf):\n",
 561     "\n",
 562     "\n",
 563     "> The remainder of a square root finder is a function _within_, which takes a tolerance and a list of approximations and looks down the list for two successive approximations that differ by no more than the given tolerance.\n",
 564     "\n",
 565     "(And note that by “list” he means a lazily-evaluated list.)\n",
 566     "\n",
 567     "Using the _output_ `[a G]` of the above generator for square root approximations, and further assuming that the first term a has been generated already and epsilon ε is handy on the stack...\n",
 568     "\n",
 569     "       a [b G] ε within\n",
 570     "    ---------------------- a b - abs ε <=\n",
 571     "          b\n",
 572     "\n",
 573     "\n",
 574     "       a [b G] ε within\n",
 575     "    ---------------------- a b - abs ε >\n",
 576     "       b [c G] ε within\n",
 577     "\n"
 578    ]
 579   },
 580   {
 581    "cell_type": "markdown",
 582    "metadata": {},
 583    "source": [
 584     "### Predicate\n",
 585     "\n",
 586     "    a [b G]             ε [first - abs] dip <=\n",
 587     "    a [b G] first - abs ε                   <=\n",
 588     "    a b           - abs ε                   <=\n",
 589     "    a-b             abs ε                   <=\n",
 590     "    abs(a-b)            ε                   <=\n",
 591     "    (abs(a-b)<=ε)"
 592    ]
 593   },
 594   {
 595    "cell_type": "code",
 596    "execution_count": null,
 597    "metadata": {},
 598    "outputs": [],
 599    "source": [
 600     "define('_within_P [first - abs] dip <=')"
 601    ]
 602   },
 603   {
 604    "cell_type": "markdown",
 605    "metadata": {},
 606    "source": [
 607     "### Base-Case\n",
 608     "\n",
 609     "    a [b G] ε roll< popop first\n",
 610     "      [b G] ε a     popop first\n",
 611     "      [b G]               first\n",
 612     "       b"
 613    ]
 614   },
 615   {
 616    "cell_type": "code",
 617    "execution_count": null,
 618    "metadata": {},
 619    "outputs": [],
 620    "source": [
 621     "define('_within_B roll< popop first')"
 622    ]
 623   },
 624   {
 625    "cell_type": "markdown",
 626    "metadata": {},
 627    "source": [
 628     "### Recur\n",
 629     "\n",
 630     "    a [b G] ε R0 [within] R1\n",
 631     "\n",
 632     "1. Discard a.\n",
 633     "2. Use `x` combinator to generate next term from `G`.\n",
 634     "3. Run `within` with `i` (it is a \"tail-recursive\" function.)\n",
 635     "\n",
 636     "Pretty straightforward:\n",
 637     "\n",
 638     "    a [b G]        ε R0           [within] R1\n",
 639     "    a [b G]        ε [popd x] dip [within] i\n",
 640     "    a [b G] popd x ε              [within] i\n",
 641     "      [b G]      x ε              [within] i\n",
 642     "    b [c G]        ε              [within] i\n",
 643     "    b [c G]        ε               within\n",
 644     "\n",
 645     "    b [c G] ε within"
 646    ]
 647   },
 648   {
 649    "cell_type": "code",
 650    "execution_count": null,
 651    "metadata": {},
 652    "outputs": [],
 653    "source": [
 654     "define('_within_R [popd x] dip')"
 655    ]
 656   },
 657   {
 658    "cell_type": "markdown",
 659    "metadata": {},
 660    "source": [
 661     "### Setting up\n",
 662     "\n",
 663     "The recursive function we have defined so far needs a slight preamble: `x` to prime the generator and the epsilon value to use:\n",
 664     "\n",
 665     "    [a G] x ε ...\n",
 666     "    a [b G] ε ..."
 667    ]
 668   },
 669   {
 670    "cell_type": "code",
 671    "execution_count": null,
 672    "metadata": {},
 673    "outputs": [],
 674    "source": [
 675     "define('within x 0.000000001 [_within_P] [_within_B] [_within_R] tailrec')\n",
 676     "define('sqrt gsra within')"
 677    ]
 678   },
 679   {
 680    "cell_type": "markdown",
 681    "metadata": {},
 682    "source": [
 683     "Try it out..."
 684    ]
 685   },
 686   {
 687    "cell_type": "code",
 688    "execution_count": null,
 689    "metadata": {
 690     "scrolled": true
 691    },
 692    "outputs": [],
 693    "source": [
 694     "J('36 sqrt')"
 695    ]
 696   },
 697   {
 698    "cell_type": "code",
 699    "execution_count": null,
 700    "metadata": {
 701     "scrolled": true
 702    },
 703    "outputs": [],
 704    "source": [
 705     "J('23 sqrt')"
 706    ]
 707   },
 708   {
 709    "cell_type": "markdown",
 710    "metadata": {},
 711    "source": [
 712     "Check it."
 713    ]
 714   },
 715   {
 716    "cell_type": "code",
 717    "execution_count": null,
 718    "metadata": {
 719     "scrolled": true
 720    },
 721    "outputs": [],
 722    "source": [
 723     "4.795831523312719**2"
 724    ]
 725   },
 726   {
 727    "cell_type": "code",
 728    "execution_count": null,
 729    "metadata": {},
 730    "outputs": [],
 731    "source": [
 732     "from math import sqrt\n",
 733     "\n",
 734     "sqrt(23)"
 735    ]
 736   }
 737  ],
 738  "metadata": {
 739   "kernelspec": {
 740    "display_name": "Joypy",
 741    "language": "",
 742    "name": "thun"
 743   },
 744   "language_info": {
 745    "file_extension": ".joy",
 746    "mimetype": "text/plain",
 747    "name": "Joy"
 748   }
 749  },
 750  "nbformat": 4,
 751  "nbformat_minor": 2
 752 }