753 lines
15 KiB
Plaintext
753 lines
15 KiB
Plaintext
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"# [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)\n",
|
||
"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",
|
||
"\n",
|
||
"Cf. [\"Why Functional Programming Matters\" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## A Generator for Approximations\n",
|
||
"\n",
|
||
"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",
|
||
"\n",
|
||
" a F\n",
|
||
" ---------\n",
|
||
" a'"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### A Function to Compute the Next Approximation\n",
|
||
"\n",
|
||
"This is the equation for computing the next approximate value of the square root:\n",
|
||
"\n",
|
||
"$a_{i+1} = \\frac{(a_i+\\frac{n}{a_i})}{2}$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"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",
|
||
"\n",
|
||
"$(a_i+\\frac{n}{a_i})$ `2 /`"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Then unadd terms:\n",
|
||
"\n",
|
||
"$a_i$ $\\frac{n}{a_i}$ `+ 2 /`"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Undivide again:\n",
|
||
"\n",
|
||
"$a_i$ $n$ $a_i$ `/ + 2 /`"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Finally deduplicate the $a_i$ term:\n",
|
||
"\n",
|
||
"$a_i$ $n$ `over / + 2 /`"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Let's try out this function `over / + 2 /` on an example:\n",
|
||
"\n",
|
||
" F == over / + 2 /"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 1,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": []
|
||
}
|
||
],
|
||
"source": [
|
||
"[F over / + 2 /] inscribe"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"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.)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 2,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
" 5 36 • F\n",
|
||
" 5 36 • over / + 2 /\n",
|
||
"5 36 5 • / + 2 /\n",
|
||
" 5 7 • + 2 /\n",
|
||
" 12 • 2 /\n",
|
||
" 12 2 • /\n",
|
||
" 6 • \n",
|
||
"\n",
|
||
"6"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"clear\n",
|
||
"\n",
|
||
"5 36 [F] trace"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"The initial estimate can be 2, and we can `cons` the input value onto a quote with `F`:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"6"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"[F1 2 swap [F] cons] inscribe"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 4,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
" 36 • F1\n",
|
||
" 36 • 2 swap [F] cons\n",
|
||
" 36 2 • swap [F] cons\n",
|
||
" 2 36 • [F] cons\n",
|
||
"2 36 [F] • cons\n",
|
||
"2 [36 F] • \n",
|
||
"\n",
|
||
"2 [36 F]"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"clear\n",
|
||
"\n",
|
||
"36 [F1] trace"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
" 2 • 36 F\n",
|
||
" 2 36 • F\n",
|
||
" 2 36 • over / + 2 /\n",
|
||
"2 36 2 • / + 2 /\n",
|
||
" 2 18 • + 2 /\n",
|
||
" 20 • 2 /\n",
|
||
" 20 2 • /\n",
|
||
" 10 • \n",
|
||
"\n",
|
||
"10"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"trace"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 6,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"6"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"36 F"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 7,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": []
|
||
}
|
||
],
|
||
"source": [
|
||
"clear"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 8,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"[2 [36 F] codireco]"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"36 F1 make_generator"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 9,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"6"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"x x x first"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"6 12"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"144 F1 make_generator x x x x first"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 11,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"2 [36 F]"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"clear\n",
|
||
"\n",
|
||
"2 [36 F]"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 12,
|
||
"metadata": {
|
||
"scrolled": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"2 [36 F] false"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"[first] [pop sqr] fork - abs 3 <"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 13,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"10"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"pop i"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 14,
|
||
"metadata": {
|
||
"scrolled": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"10 [36 F] false"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"[36 F] [first] [pop sqr] fork - abs 3 <"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 15,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"6"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"pop i"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 16,
|
||
"metadata": {
|
||
"scrolled": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"6 [36 F] true"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"[36 F] [first] [pop sqr] fork - abs 3 <"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 17,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"2"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"clear\n",
|
||
"\n",
|
||
"2"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 18,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"6"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"[] true [i [36 F] [first] [pop sqr] fork - abs 3 >] loop pop"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 19,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"12"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"clear\n",
|
||
"\n",
|
||
"7 [] true [i [144 F] [first] [pop sqr] fork - abs 3 >] loop pop"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 20,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"120"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"clear\n",
|
||
"\n",
|
||
"7 [] true [i [14400 F] [first] [pop sqr] fork - abs 3 >] loop pop"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"broken due to no float div\n",
|
||
"\n",
|
||
" clear\n",
|
||
"\n",
|
||
" 7 [] true [i [1000 F] [first] [pop sqr] fork - abs 10 >] loop pop"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Make it into a Generator\n",
|
||
"\n",
|
||
"Our generator would be created by:\n",
|
||
"\n",
|
||
" a [dup F] make_generator\n",
|
||
"\n",
|
||
"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",
|
||
"\n",
|
||
" 1 n 1 / + 2 /\n",
|
||
" 1 n/1 + 2 /\n",
|
||
" 1 n + 2 /\n",
|
||
" n+1 2 /\n",
|
||
" (n+1)/2\n",
|
||
"\n",
|
||
"The generator can be written as:\n",
|
||
"\n",
|
||
" 23 1 swap [over / + 2 /] cons [dup] swoncat make_generator\n",
|
||
" 1 23 [over / + 2 /] cons [dup] swoncat make_generator\n",
|
||
" 1 [23 over / + 2 /] [dup] swoncat make_generator\n",
|
||
" 1 [dup 23 over / + 2 /] make_generator"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {
|
||
"scrolled": true
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"define('gsra 1 swap [over / + 2 /] cons [dup] swoncat make_generator')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"J('23 gsra')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Let's drive the generator a few time (with the `x` combinator) and square the approximation to see how well it works..."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"J('23 gsra 6 [x popd] times first sqr')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Finding Consecutive Approximations within a Tolerance\n",
|
||
"\n",
|
||
"From [\"Why Functional Programming Matters\" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf):\n",
|
||
"\n",
|
||
"\n",
|
||
"> 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",
|
||
"\n",
|
||
"(And note that by “list” he means a lazily-evaluated list.)\n",
|
||
"\n",
|
||
"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",
|
||
"\n",
|
||
" a [b G] ε within\n",
|
||
" ---------------------- a b - abs ε <=\n",
|
||
" b\n",
|
||
"\n",
|
||
"\n",
|
||
" a [b G] ε within\n",
|
||
" ---------------------- a b - abs ε >\n",
|
||
" b [c G] ε within\n",
|
||
"\n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Predicate\n",
|
||
"\n",
|
||
" a [b G] ε [first - abs] dip <=\n",
|
||
" a [b G] first - abs ε <=\n",
|
||
" a b - abs ε <=\n",
|
||
" a-b abs ε <=\n",
|
||
" abs(a-b) ε <=\n",
|
||
" (abs(a-b)<=ε)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"define('_within_P [first - abs] dip <=')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Base-Case\n",
|
||
"\n",
|
||
" a [b G] ε roll< popop first\n",
|
||
" [b G] ε a popop first\n",
|
||
" [b G] first\n",
|
||
" b"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"define('_within_B roll< popop first')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Recur\n",
|
||
"\n",
|
||
" a [b G] ε R0 [within] R1\n",
|
||
"\n",
|
||
"1. Discard a.\n",
|
||
"2. Use `x` combinator to generate next term from `G`.\n",
|
||
"3. Run `within` with `i` (it is a \"tail-recursive\" function.)\n",
|
||
"\n",
|
||
"Pretty straightforward:\n",
|
||
"\n",
|
||
" a [b G] ε R0 [within] R1\n",
|
||
" a [b G] ε [popd x] dip [within] i\n",
|
||
" a [b G] popd x ε [within] i\n",
|
||
" [b G] x ε [within] i\n",
|
||
" b [c G] ε [within] i\n",
|
||
" b [c G] ε within\n",
|
||
"\n",
|
||
" b [c G] ε within"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"define('_within_R [popd x] dip')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Setting up\n",
|
||
"\n",
|
||
"The recursive function we have defined so far needs a slight preamble: `x` to prime the generator and the epsilon value to use:\n",
|
||
"\n",
|
||
" [a G] x ε ...\n",
|
||
" a [b G] ε ..."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"define('within x 0.000000001 [_within_P] [_within_B] [_within_R] tailrec')\n",
|
||
"define('sqrt gsra within')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Try it out..."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {
|
||
"scrolled": true
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"J('36 sqrt')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {
|
||
"scrolled": true
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"J('23 sqrt')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Check it."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {
|
||
"scrolled": true
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"4.795831523312719**2"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"from math import sqrt\n",
|
||
"\n",
|
||
"sqrt(23)"
|
||
]
|
||
}
|
||
],
|
||
"metadata": {
|
||
"kernelspec": {
|
||
"display_name": "Joypy",
|
||
"language": "",
|
||
"name": "thun"
|
||
},
|
||
"language_info": {
|
||
"file_extension": ".joy",
|
||
"mimetype": "text/plain",
|
||
"name": "Joy"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 2
|
||
}
|