From 8785b24e0493e45893f13b07ececc08b1a2f46fe Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Sun, 10 Aug 2008 16:44:17 -0500 Subject: [PATCH 01/11] Now with arbitrary accuracy --- .../math/derivatives/derivatives-docs.factor | 50 +++++++++++- extra/math/derivatives/derivatives.factor | 80 +++++++++++++++++-- 2 files changed, 119 insertions(+), 11 deletions(-) diff --git a/extra/math/derivatives/derivatives-docs.factor b/extra/math/derivatives/derivatives-docs.factor index 23847e82f7..70389f18ad 100644 --- a/extra/math/derivatives/derivatives-docs.factor +++ b/extra/math/derivatives/derivatives-docs.factor @@ -3,7 +3,51 @@ USING: help.markup help.syntax ; IN: math.derivatives HELP: derivative ( x function -- m ) -{ $values { "x" "the x-position on the function" } { "function" "a differentiable function" } } -{ $description "Finds the slope of the tangent line at the given x-position on the given function." } ; +{ $values { "x" "a position on the function" } { "function" "a differentiable function" } } +{ $description + "Approximates the slope of the tangent line by using Ridders' " + "method of computing derivatives, from the chapter \"Accurate computation " + "of F'(x) and F'(x)F''(x)\", from \"Advances in Engineering Software, Vol. 4, pp. 75-76 ." +} +{ $examples + { $example + "USING: math.derivatives prettyprint ;" + "[ sq ] 4 derivative ." + "8" + } + { $notes + "For applied scientists, you may play with the settings " + "in the source file to achieve arbitrary accuracy." + } +} ; -{ derivative-func } related-words +HELP: fast-derivative ( x function -- m ) +{ $values { "x" "a x-position on the function" } { "function" "a differentiable function" } } +{ $description + "Approximates the slope of the tangent line of the provided function " + "by using a secant line with very near points. This implementation is " + "naive and is only provided because it is used in the much more " + "accurate " { $link derivative } " word. Use this word if accuracy " + "is of no importance." +} ; + +HELP: derivative-func ( function -- der ) +{ $values { "func" "a differentiable function" } { "der" "the derivative" } } +{ $description + "Provides the derivative of the function. The implementation simply " + "attaches the " { $link derivative } " word to the end of the function." +} +{ $examples + { $example + "USING: math.derivatives prettyprint ;" + "[ sq ] derivative-func ." + "[ [ sq ] derivative ]" + } +} ; + +ARTICLE: "derivatives" "The Derivative Toolkit" +"A toolkit for computing the derivative of functions." +{ $subsection derivative } +{ $subsection derivative-func } +{ $subsection fast-derivative } ; +ABOUT: "derivatives" diff --git a/extra/math/derivatives/derivatives.factor b/extra/math/derivatives/derivatives.factor index d92066efaf..f77748d0b5 100644 --- a/extra/math/derivatives/derivatives.factor +++ b/extra/math/derivatives/derivatives.factor @@ -1,10 +1,74 @@ -! Copyright © 2008 Reginald Keith Ford II -! Tool for computing the derivative of a function at a point -USING: kernel math math.points math.function-tools ; +! Copyright (c) 2008 Reginald Ford +! Tools for approximating derivatives + +USING: kernel math math.functions locals generalizations float-arrays sequences +math.constants namespaces math.function-tools math.points math.ranges math.order ; IN: math.derivatives -: small-amount ( -- n ) 1.0e-14 ; -: some-more ( x -- y ) small-amount + ; -: some-less ( x -- y ) small-amount - ; -: derivative ( x function -- m ) [ [ some-more ] dip eval ] [ [ some-less ] dip eval ] 2bi slope ; -: derivative-func ( function -- function ) [ derivative ] curry ; \ No newline at end of file +! Ridders' method of a derivative, from the chapter +! "Accurate computation of F'(x) and F'(x)F''(x)", +! From "Advances in Engineering Software, Vol. 4, pp. 75-76 +! \ fast-derivative has been factored out for use by children + +: largest-float ( -- x ) HEX: 7fefffffffffffff bits>double ; foldable +: ntab 10 ; ! max size of tableau (main accuracy setting) +: con 1.41 ; ! stepsize is decreased by this per-iteration +: con2 1.9881 ; ! con^2 +: initial-h 0.02 ; ! distance of the 2 points of the first secant line +: safe 2.0 ; ! return when current err is SAFE worse than the best + ! \ safe probably should not be changed +SYMBOL: i +SYMBOL: j +SYMBOL: err +SYMBOL: errt +SYMBOL: fac +SYMBOL: h +SYMBOL: ans +SYMBOL: matrix + +: (derivative) ( x function -- m ) + [ [ h get + ] dip eval ] + [ [ h get - ] dip eval ] + 2bi slope ; inline +: fast-derivative ( x function -- m ) + over epsilon sqrt * h set + (derivative) ; inline +: init-matrix ( -- ) + ntab [ ntab ] replicate + matrix set ; +: m-set ( value j i -- ) matrix get nth set-nth ; +: m-get ( j i -- n ) matrix get nth nth ; +:: derivative ( x func -- m ) + init-matrix + initial-h h set + x func (derivative) 0 0 m-set + largest-float err set + ntab 1 - [1,b] [| i | + h [ con / ] change + x func (derivative) 0 i m-set + con2 fac set + i [1,b] [| j | + j 1 - i m-get fac get * + j 1 - i 1 - m-get + - + fac get 1 - + / j i m-set + fac [ con2 * ] change + j i m-get j 1 - i m-get - abs + j i m-get j 1 - i 1 - m-get - abs + max errt set + errt get err get <= + [ + errt get err set + j i m-get ans set + ] [ ] + if + ] each + i i m-get i 1 - dup m-get - abs + err get safe * + < + ] all? drop + ans get ; inline +: derivative-func ( function -- function ) [ derivative ] curry ; inline +: fast-derivative-func ( function -- function ) [ fast-derivative ] curry ; inline + From aee8dbdba45abf06eab2e974b2985bc4e9c41ed5 Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Sun, 10 Aug 2008 16:45:13 -0500 Subject: [PATCH 02/11] peer review by myself --- extra/24-game/24-game-docs.factor | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/extra/24-game/24-game-docs.factor b/extra/24-game/24-game-docs.factor index 12a558b2d2..cd82f335d8 100644 --- a/extra/24-game/24-game-docs.factor +++ b/extra/24-game/24-game-docs.factor @@ -31,12 +31,12 @@ HELP: 24-able ( -- vector ) "just using the provided commands and the 4 numbers. The Following are the " "provided commands: " { $link + } ", " { $link - } ", " { $link * } ", " - { $link / } ", and " { $link swap } "." + { $link / } ", " { $link swap } ", and " { $link rot } "." } { $examples { $example "USE: 24-game" - "24-able vector-24-able?" + "24-able vector-24-able? ." "t" } { $notes { $link 24-able? } " is used in " { $link 24-able } "." } From 9d0acc555d51b1aba9918b8fd73269f6a5bd40cb Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Sun, 10 Aug 2008 16:47:52 -0500 Subject: [PATCH 03/11] peer review by myself --- extra/animations/animations-docs.factor | 63 ++++++++++++++++++------- 1 file changed, 47 insertions(+), 16 deletions(-) diff --git a/extra/animations/animations-docs.factor b/extra/animations/animations-docs.factor index 6a1e89a28e..000c0ce4cc 100644 --- a/extra/animations/animations-docs.factor +++ b/extra/animations/animations-docs.factor @@ -1,34 +1,65 @@ USING: help.markup help.syntax ; -IN: extra.animations +IN: animations HELP: animate ( quot duration -- ) + { $values { "quot" "a quot which uses " { $link progress } } { "duration" "a duration of time" } } -{ $description { $link animate } " calls " { $link reset-progress } " , then continously calls the given quot until the duration of time has elapsed. The quot should use " { $link progress } " at least once." } -{ $example - "USING: extra.animations calendar threads prettyprint ;" - "[ 1 sleep progress unparse write \" ms elapsed\" print ] 1/20 seconds animate ;" - "46 ms elapsed\n17 ms elapsed" +{ $description + { $link animate } " calls " { $link reset-progress } + " , then continously calls the given quot until the" + " duration of time has elapsed. The quot should use " + { $link progress } " at least once." +} +{ $examples + { $unchecked-example + "USING: animations calendar threads prettyprint ;" + "[ 1 sleep progress unparse write \" ms elapsed\" print ] " + "1/20 seconds animate ;" + "46 ms elapsed\n17 ms elapsed" + } + { $notes "The amount of time elapsed between these iterations will very." } } ; HELP: reset-progress ( -- ) -{ $description "Initiates the timer. Call this before using a loop which makes use of " { $link progress } "." } ; +{ $description + "Initiates the timer. Call this before using " + "a loop which makes use of " { $link progress } "." +} ; HELP: progress ( -- time ) { $values { "time" "an integer" } } -{ $description "Gives the time elapsed since the last time this word was called, in milliseconds." } -{ $example - "USING: extra.animations threads prettyprint ;" - "reset-progress 3 [ 1 sleep progress unparse write \"ms elapsed\" print ] times ;" - "31 ms elapsed\n18 ms elapsed\n16 ms elapsed" +{ $description + "Gives the time elapsed since the last time" + " this word was called, in milliseconds." +} +{ $examples + { $unchecked-example + "USING: animations threads prettyprint ;" + "reset-progress 3 " + "[ 1 sleep progress unparse write \"ms elapsed\" print ] " + "times ;" + "31 ms elapsed\n18 ms elapsed\n16 ms elapsed" + } + { $notes "The amount of time elapsed between these iterations will very." } } ; -ARTICLE: "extra.animations" "Animations" -"Provides a lightweight framework for properly simulating continuous functions of real time. This framework helps one create animations that use rates which do not change across platforms. The speed of the computer should correlate with the smoothness of the animation, not the speed of the animation!" +ARTICLE: "animations" "Animations" +"Provides a lightweight framework for properly simulating continuous" +" functions of real time. This framework helps one create animations " +"that use rates which do not change across platforms. The speed of the " +"computer should correlate with the smoothness of the animation, not " +"the speed of the animation!" { $subsection animate } { $subsection reset-progress } { $subsection progress } -{ $link progress } " specifically provides the length of time since " { $link reset-progress } " was called, and also calls " { $link reset-progress } " as its last action. This can be directly used when one's quote runs for a specific number of iterations, instead of a length of time. If the animation is like most, and is expected to run for a specific length of time, " { $link animate } " should be used." ; -ABOUT: "extra.animations" \ No newline at end of file +! A little talk about when to use progress and when to use animate + { $link progress } " specifically provides the length of time since " + { $link reset-progress } " was called, and also calls " + { $link reset-progress } " as its last action. This can be directly " + "used when one's quote runs for a specific number of iterations, instead " + "of a length of time. If the animation is like most, and is expected to " + "run for a specific length of time, " { $link animate } " should be used." ; +ABOUT: "animations" \ No newline at end of file From 6ab0f6b09c633f6543b1feee16babefd8aab287a Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Sun, 10 Aug 2008 16:48:54 -0500 Subject: [PATCH 04/11] No one else used middle name --- extra/animations/authors.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extra/animations/authors.txt b/extra/animations/authors.txt index dac0cb42fe..137b1605da 100644 --- a/extra/animations/authors.txt +++ b/extra/animations/authors.txt @@ -1 +1 @@ -Reginald Keith Ford II \ No newline at end of file +Reginald Ford \ No newline at end of file From a44097af93cebbef61d8bee80da94c93618b3b49 Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Sun, 10 Aug 2008 16:49:40 -0500 Subject: [PATCH 05/11] combinators should inline --- extra/math/function-tools/function-tools.factor | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/extra/math/function-tools/function-tools.factor b/extra/math/function-tools/function-tools.factor index 802bf9e14e..ec93a0891a 100644 --- a/extra/math/function-tools/function-tools.factor +++ b/extra/math/function-tools/function-tools.factor @@ -3,7 +3,7 @@ USING: kernel math arrays sequences sequences.lib ; IN: math.function-tools -: difference-func ( func func -- func ) [ bi - ] 2curry ; -: eval ( x func -- pt ) dupd call 2array ; -: eval-inverse ( y func -- pt ) dupd call swap 2array ; -: eval3d ( x y func -- pt ) [ 2dup ] dip call 3array ; +: difference-func ( func func -- func ) [ bi - ] 2curry ; inline +: eval ( x func -- pt ) dupd call 2array ; inline +: eval-inverse ( y func -- pt ) dupd call swap 2array ; inline +: eval3d ( x y func -- pt ) [ 2dup ] dip call 3array ; inline From 6df077805d5a89972aa1679a78565645c3a525ad Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Sun, 10 Aug 2008 18:20:14 -0500 Subject: [PATCH 06/11] minor fixes --- extra/math/derivatives/derivatives-docs.factor | 15 ++------------- extra/math/derivatives/derivatives.factor | 16 +++------------- 2 files changed, 5 insertions(+), 26 deletions(-) diff --git a/extra/math/derivatives/derivatives-docs.factor b/extra/math/derivatives/derivatives-docs.factor index 70389f18ad..0db52adfa5 100644 --- a/extra/math/derivatives/derivatives-docs.factor +++ b/extra/math/derivatives/derivatives-docs.factor @@ -17,20 +17,10 @@ HELP: derivative ( x function -- m ) } { $notes "For applied scientists, you may play with the settings " - "in the source file to achieve arbitrary accuracy." + "in the source file to achieve arbitrary accuracy. " } } ; -HELP: fast-derivative ( x function -- m ) -{ $values { "x" "a x-position on the function" } { "function" "a differentiable function" } } -{ $description - "Approximates the slope of the tangent line of the provided function " - "by using a secant line with very near points. This implementation is " - "naive and is only provided because it is used in the much more " - "accurate " { $link derivative } " word. Use this word if accuracy " - "is of no importance." -} ; - HELP: derivative-func ( function -- der ) { $values { "func" "a differentiable function" } { "der" "the derivative" } } { $description @@ -48,6 +38,5 @@ HELP: derivative-func ( function -- der ) ARTICLE: "derivatives" "The Derivative Toolkit" "A toolkit for computing the derivative of functions." { $subsection derivative } -{ $subsection derivative-func } -{ $subsection fast-derivative } ; +{ $subsection derivative-func } ; ABOUT: "derivatives" diff --git a/extra/math/derivatives/derivatives.factor b/extra/math/derivatives/derivatives.factor index f77748d0b5..96d0fc3a81 100644 --- a/extra/math/derivatives/derivatives.factor +++ b/extra/math/derivatives/derivatives.factor @@ -1,19 +1,13 @@ -! Copyright (c) 2008 Reginald Ford ! Tools for approximating derivatives USING: kernel math math.functions locals generalizations float-arrays sequences math.constants namespaces math.function-tools math.points math.ranges math.order ; IN: math.derivatives -! Ridders' method of a derivative, from the chapter -! "Accurate computation of F'(x) and F'(x)F''(x)", -! From "Advances in Engineering Software, Vol. 4, pp. 75-76 -! \ fast-derivative has been factored out for use by children - : largest-float ( -- x ) HEX: 7fefffffffffffff bits>double ; foldable -: ntab 10 ; ! max size of tableau (main accuracy setting) -: con 1.41 ; ! stepsize is decreased by this per-iteration -: con2 1.9881 ; ! con^2 +: ntab 10 ; ! max size of tableau (main accuracy setting) +: con 1.41 ; ! stepsize is decreased by this per-iteration +: con2 1.9881 ; ! con^2 : initial-h 0.02 ; ! distance of the 2 points of the first secant line : safe 2.0 ; ! return when current err is SAFE worse than the best ! \ safe probably should not be changed @@ -30,9 +24,6 @@ SYMBOL: matrix [ [ h get + ] dip eval ] [ [ h get - ] dip eval ] 2bi slope ; inline -: fast-derivative ( x function -- m ) - over epsilon sqrt * h set - (derivative) ; inline : init-matrix ( -- ) ntab [ ntab ] replicate matrix set ; @@ -70,5 +61,4 @@ SYMBOL: matrix ] all? drop ans get ; inline : derivative-func ( function -- function ) [ derivative ] curry ; inline -: fast-derivative-func ( function -- function ) [ fast-derivative ] curry ; inline From 6060b12ccb2d8e6a9ab4aafbe7d24e5e77cc75bf Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Sun, 10 Aug 2008 18:22:32 -0500 Subject: [PATCH 07/11] minor additions --- extra/animations/animations.factor | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/extra/animations/animations.factor b/extra/animations/animations.factor index 7efd618bbf..db5b3448c1 100644 --- a/extra/animations/animations.factor +++ b/extra/animations/animations.factor @@ -2,11 +2,14 @@ USING: kernel shuffle system locals prettyprint math io namespaces threads calendar ; -IN: extra.animations +IN: animations SYMBOL: last-loop +SYMBOL: sleep-period + : reset-progress ( -- ) millis last-loop set ; : progress ( -- progress ) millis last-loop get - reset-progress ; : set-end ( duration -- end-time ) dt>milliseconds millis + ; -: loop ( quot end -- ) dup millis > [ [ dup call ] dip loop ] [ 2drop ] if ; -: animate ( quot duration -- ) reset-progress set-end loop ; \ No newline at end of file +: loop ( quot end -- ) dup millis > [ [ dup call ] dip loop ] [ 2drop ] if ; inline +: animate ( quot duration -- ) reset-progress set-end loop ; inline +: sample ( revs quot -- avg ) reset-progress dupd times progress swap / ; inline \ No newline at end of file From bd168d06f258860a2d8c074d1e0335b02d6e1ec5 Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Tue, 12 Aug 2008 00:28:22 -0400 Subject: [PATCH 08/11] now with progress-peek --- extra/animations/animations.factor | 2 ++ 1 file changed, 2 insertions(+) diff --git a/extra/animations/animations.factor b/extra/animations/animations.factor index db5b3448c1..803536a51c 100644 --- a/extra/animations/animations.factor +++ b/extra/animations/animations.factor @@ -8,7 +8,9 @@ SYMBOL: last-loop SYMBOL: sleep-period : reset-progress ( -- ) millis last-loop set ; +! : my-progress ( -- progress ) millis : progress ( -- progress ) millis last-loop get - reset-progress ; +: progress-peek ( -- progress ) millis last-loop get - ; : set-end ( duration -- end-time ) dt>milliseconds millis + ; : loop ( quot end -- ) dup millis > [ [ dup call ] dip loop ] [ 2drop ] if ; inline : animate ( quot duration -- ) reset-progress set-end loop ; inline From 2271aae7f070b5c373547f122d07861d239e74d1 Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Tue, 12 Aug 2008 02:42:23 -0400 Subject: [PATCH 09/11] compatible with demos menu --- extra/24-game/24-game.factor | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/extra/24-game/24-game.factor b/extra/24-game/24-game.factor index 52f0cd6833..126215ab13 100644 --- a/extra/24-game/24-game.factor +++ b/extra/24-game/24-game.factor @@ -59,4 +59,5 @@ DEFER: check-status : 24-able? ( vector -- t/f ) [ makes-24? ] with-datastack first ; : 24-able ( -- vector ) build-quad dup 24-able? [ drop build-quad ] unless ; : set-commands ( -- ) { + - * / rot swap q } commands set ; -: play-game ( -- ) set-commands 24-able repeat ; \ No newline at end of file +: play-game ( -- ) set-commands 24-able repeat ; +MAIN: play-game \ No newline at end of file From 359bff5f154caaaaa9ed5037e78d43033cdf48b3 Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Tue, 12 Aug 2008 11:24:00 -0400 Subject: [PATCH 10/11] Derivatives without dynamics OR locals --- extra/math/derivatives/authors.txt | 3 +- .../math/derivatives/derivatives-docs.factor | 66 ++++++- extra/math/derivatives/derivatives.factor | 175 ++++++++++++------ 3 files changed, 181 insertions(+), 63 deletions(-) diff --git a/extra/math/derivatives/authors.txt b/extra/math/derivatives/authors.txt index 137b1605da..3be8a6d4d3 100644 --- a/extra/math/derivatives/authors.txt +++ b/extra/math/derivatives/authors.txt @@ -1 +1,2 @@ -Reginald Ford \ No newline at end of file +Reginald Ford +Eduardo Cavazos \ No newline at end of file diff --git a/extra/math/derivatives/derivatives-docs.factor b/extra/math/derivatives/derivatives-docs.factor index 0db52adfa5..a78a697b76 100644 --- a/extra/math/derivatives/derivatives-docs.factor +++ b/extra/math/derivatives/derivatives-docs.factor @@ -1,5 +1,4 @@ -USING: help.markup help.syntax ; - +USING: help.markup help.syntax math.functions ; IN: math.derivatives HELP: derivative ( x function -- m ) @@ -21,6 +20,46 @@ HELP: derivative ( x function -- m ) } } ; +HELP: (derivative) ( x function h err -- m ) +{ $values + { "x" "a position on the function" } + { "function" "a differentiable function" } + { + "h" "distance between the points of the first secant line used for " + "approximation of the tangent. This distance will be divided " + "constantly, by " { $link con } ". See " { $link init-hh } + " for the code which enforces this. H should be .001 to .5 -- too " + "small can cause bad convergence. Also, h should be small enough " + "to give the correct sgn(f'(x)). In other words, if you're expecting " + "a positive derivative, make h small enough to give the same " + "when plugged into the academic limit definition of a derivative. " + "See " { $link update-hh } " for the code which performs this task." + } + { + "err" "maximum tolerance of increase in error. For example, if this " + "is set to 2.0, the program will terminate with its nearest answer " + "when the error multiplies by 2. See " { $link check-safe } " for " + "the enforcing code." + } +} +{ $description + "Approximates the slope of the tangent line by using Ridders' " + "method of computing derivatives, from the chapter \"Accurate computation " + "of F'(x) and F'(x)F''(x)\", from \"Advances in Engineering Software, " + "Vol. 4, pp. 75-76 ." +} +{ $examples + { $example + "USING: math.derivatives prettyprint ;" + "[ sq ] 4 derivative ." + "8" + } + { $notes + "For applied scientists, you may play with the settings " + "in the source file to achieve arbitrary accuracy. " + } +} ; + HELP: derivative-func ( function -- der ) { $values { "func" "a differentiable function" } { "der" "the derivative" } } { $description @@ -30,8 +69,27 @@ HELP: derivative-func ( function -- der ) { $examples { $example "USING: math.derivatives prettyprint ;" - "[ sq ] derivative-func ." - "[ [ sq ] derivative ]" + "60 deg>rad [ sin ] derivative-func call ." + "0.5000000000000173" + } + { $notes + "Without a heavy algebraic system, derivatives must be " + "approximated. With the current settings, there is a fair trade of " + "speed and accuracy; the first 12 digits " + "will always be correct with " { $link sin } " and " { $link cos } + ". The following code performs a minumum and maximum error test." + { $code + "USING: kernel math math.functions math.trig sequences sequences.lib ;" + "360" + "[" + " deg>rad" + " [ [ sin ] derivative-func call ]" + " ! Note: the derivative of sin is cos" + " [ cos ]" + " bi - abs" + "] map minmax" + + } } } ; diff --git a/extra/math/derivatives/derivatives.factor b/extra/math/derivatives/derivatives.factor index 96d0fc3a81..ad8d944bfe 100644 --- a/extra/math/derivatives/derivatives.factor +++ b/extra/math/derivatives/derivatives.factor @@ -1,64 +1,123 @@ -! Tools for approximating derivatives -USING: kernel math math.functions locals generalizations float-arrays sequences -math.constants namespaces math.function-tools math.points math.ranges math.order ; +USING: kernel continuations combinators sequences math + math.order math.ranges accessors float-arrays ; + IN: math.derivatives +TUPLE: state x func h err i j errt fac hh ans a done ; + : largest-float ( -- x ) HEX: 7fefffffffffffff bits>double ; foldable -: ntab 10 ; ! max size of tableau (main accuracy setting) -: con 1.41 ; ! stepsize is decreased by this per-iteration -: con2 1.9881 ; ! con^2 -: initial-h 0.02 ; ! distance of the 2 points of the first secant line -: safe 2.0 ; ! return when current err is SAFE worse than the best - ! \ safe probably should not be changed -SYMBOL: i -SYMBOL: j -SYMBOL: err -SYMBOL: errt -SYMBOL: fac -SYMBOL: h -SYMBOL: ans -SYMBOL: matrix +: ntab ( -- val ) 8 ; +: con ( -- val ) 1.6 ; +: con2 ( -- val ) con con * ; +: big ( -- val ) largest-float ; +: safe ( -- val ) 2.0 ; -: (derivative) ( x function -- m ) - [ [ h get + ] dip eval ] - [ [ h get - ] dip eval ] - 2bi slope ; inline -: init-matrix ( -- ) - ntab [ ntab ] replicate - matrix set ; -: m-set ( value j i -- ) matrix get nth set-nth ; -: m-get ( j i -- n ) matrix get nth nth ; -:: derivative ( x func -- m ) - init-matrix - initial-h h set - x func (derivative) 0 0 m-set - largest-float err set - ntab 1 - [1,b] [| i | - h [ con / ] change - x func (derivative) 0 i m-set - con2 fac set - i [1,b] [| j | - j 1 - i m-get fac get * - j 1 - i 1 - m-get - - - fac get 1 - - / j i m-set - fac [ con2 * ] change - j i m-get j 1 - i m-get - abs - j i m-get j 1 - i 1 - m-get - abs - max errt set - errt get err get <= - [ - errt get err set - j i m-get ans set - ] [ ] - if - ] each - i i m-get i 1 - dup m-get - abs - err get safe * - < - ] all? drop - ans get ; inline -: derivative-func ( function -- function ) [ derivative ] curry ; inline +! Yes, this was ported from C code. +: a[i][i] ( state -- elt ) [ i>> ] [ i>> ] [ a>> ] tri nth nth ; +: a[j][i] ( state -- elt ) [ i>> ] [ j>> ] [ a>> ] tri nth nth ; +: a[j-1][i] ( state -- elt ) [ i>> ] [ j>> 1 - ] [ a>> ] tri nth nth ; +: a[j-1][i-1] ( state -- elt ) [ i>> 1 - ] [ j>> 1 - ] [ a>> ] tri nth nth ; +: a[i-1][i-1] ( state -- elt ) [ i>> 1 - ] [ i>> 1 - ] [ a>> ] tri nth nth ; +: check-h ( state -- state ) + dup h>> 0 = [ "h must be nonzero in dfridr" throw ] when ; +: init-a ( state -- state ) ntab [ ntab ] replicate >>a ; +: init-hh ( state -- state ) dup h>> >>hh ; +: init-err ( state -- state ) big >>err ; +: update-hh ( state -- state ) dup hh>> con / >>hh ; +: reset-fac ( state -- state ) con2 >>fac ; +: update-fac ( state -- state ) dup fac>> con2 * >>fac ; + +! If error is decreased, save the improved answer +: error-decreased? ( state -- state ? ) [ ] [ errt>> ] [ err>> ] tri <= ; +: save-improved-answer ( state -- state ) + dup err>> >>errt + dup a[j][i] >>ans ; + +! If higher order is worse by a significant factor SAFE, then quit early. +: check-safe ( state -- state ) + dup + [ [ a[i][i] ] [ a[i-1][i-1] ] bi - abs ] [ err>> safe * ] bi >= + [ t >>done ] + when ; +: x+hh ( state -- val ) [ x>> ] [ hh>> ] bi + ; +: x-hh ( state -- val ) [ x>> ] [ hh>> ] bi - ; +: limit-approx ( state -- val ) + [ + [ [ x+hh ] [ func>> ] bi call ] + [ [ x-hh ] [ func>> ] bi call ] + bi - + ] + [ hh>> 2.0 * ] + bi / ; +: a[0][0]! ( state -- state ) + { [ ] [ limit-approx ] [ drop 0 ] [ drop 0 ] [ a>> ] } cleave nth set-nth ; +: a[0][i]! ( state -- state ) + { [ ] [ limit-approx ] [ i>> ] [ drop 0 ] [ a>> ] } cleave nth set-nth ; +: a[j-1][i]*fac ( state -- val ) [ a[j-1][i] ] [ fac>> ] bi * ; +: new-a[j][i] ( state -- val ) + [ [ a[j-1][i]*fac ] [ a[j-1][i-1] ] bi - ] + [ fac>> 1.0 - ] + bi / ; +: a[j][i]! ( state -- state ) + { [ ] [ new-a[j][i] ] [ i>> ] [ j>> ] [ a>> ] } cleave nth set-nth ; + +: update-errt ( state -- state ) + dup + [ [ a[j][i] ] [ a[j-1][i] ] bi - abs ] + [ [ a[j][i] ] [ a[j-1][i-1] ] bi - abs ] + bi max + >>errt ; + +: not-done? ( state -- state ? ) dup done>> not ; + +: derive ( state -- state ) + init-a + check-h + init-hh + a[0][0]! + init-err + 1 ntab [a,b) + [ + >>i + not-done? + [ + update-hh + a[0][i]! + reset-fac + 1 over i>> [a,b] + [ + >>j + a[j][i]! + update-fac + update-errt + error-decreased? [ save-improved-answer ] when + ] + each + check-safe + ] + when + ] + each ; + +: derivative-state ( x func h err -- state ) + state new + swap >>err + swap >>h + swap >>func + swap >>x ; + +! For scientists: +! h should be .001 to .5 -- too small can cause bad convergence, +! h should be small enough to give the correct sgn(f'(x)) +! err is the max tolerance of gain in error for a single iteration- +: (derivative) ( x func h err -- ans error ) + derivative-state + derive + [ ans>> ] + [ errt>> ] + bi ; + +: derivative ( x func -- m ) 0.01 2.0 (derivative) drop ; +: derivative-func ( func -- der ) [ derivative ] curry ; \ No newline at end of file From 6f767add2c4fada948553e2639888a6539d7ae1c Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Tue, 12 Aug 2008 12:00:54 -0400 Subject: [PATCH 11/11] documentation for scientists --- extra/math/derivatives/derivatives-docs.factor | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/extra/math/derivatives/derivatives-docs.factor b/extra/math/derivatives/derivatives-docs.factor index a78a697b76..15dd954b1c 100644 --- a/extra/math/derivatives/derivatives-docs.factor +++ b/extra/math/derivatives/derivatives-docs.factor @@ -96,5 +96,6 @@ HELP: derivative-func ( function -- der ) ARTICLE: "derivatives" "The Derivative Toolkit" "A toolkit for computing the derivative of functions." { $subsection derivative } -{ $subsection derivative-func } ; +{ $subsection derivative-func } +{ $subsection (derivative) } ; ABOUT: "derivatives"