Now with arbitrary accuracy

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" } }
+{ $values { "x" "a 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." } ;
+{ $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"

extra/math/derivatives/derivatives.factor

@@ -1,10 +1,74 @@
-! Copyright © 2008 Reginald Keith Ford II
+! Copyright (c) 2008 Reginald Ford
-! Tool for computing the derivative of a function at a point 
+! Tools for approximating derivatives
-USING: kernel math math.points math.function-tools ;
+
+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 ;
+! Ridders' method of a derivative, from the chapter
-: some-more ( x -- y ) small-amount + ;
+! "Accurate computation of F'(x) and F'(x)F''(x)",
-: some-less ( x -- y ) small-amount - ;
+! From "Advances in Engineering Software, Vol. 4, pp. 75-76
-: derivative ( x function -- m ) [ [ some-more ] dip eval ] [ [ some-less ] dip eval ] 2bi slope ;
+! \ fast-derivative has been factored out for use by children
-: derivative-func ( function -- function ) [ derivative ] curry ;
+
+: 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 <float-array> ] 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
+