From fe55e939f95b198602ca12659fc13759c33a6fc5 Mon Sep 17 00:00:00 2001 From: Jason Merrill Date: Thu, 12 Feb 2009 23:13:16 -0500 Subject: [PATCH 1/2] Added math.dual and math.derivatives for computing with dual numbers. Also made a few more methods in math.functions generic in order to specialize them on dual numbers. --- basis/math/functions/functions.factor | 8 +- extra/math/derivatives/authors.txt | 1 + .../math/derivatives/derivatives-docs.factor | 11 +++ .../math/derivatives/derivatives-tests.factor | 4 + extra/math/derivatives/derivatives.factor | 34 +++++++ extra/math/derivatives/syntax/authors.txt | 1 + .../derivatives/syntax/syntax-docs.factor | 18 ++++ extra/math/derivatives/syntax/syntax.factor | 10 +++ extra/math/dual/authors.txt | 1 + extra/math/dual/dual-docs.factor | 90 +++++++++++++++++++ extra/math/dual/dual-tests.factor | 14 +++ extra/math/dual/dual.factor | 80 +++++++++++++++++ 12 files changed, 270 insertions(+), 2 deletions(-) create mode 100644 extra/math/derivatives/authors.txt create mode 100644 extra/math/derivatives/derivatives-docs.factor create mode 100644 extra/math/derivatives/derivatives-tests.factor create mode 100644 extra/math/derivatives/derivatives.factor create mode 100644 extra/math/derivatives/syntax/authors.txt create mode 100644 extra/math/derivatives/syntax/syntax-docs.factor create mode 100644 extra/math/derivatives/syntax/syntax.factor create mode 100644 extra/math/dual/authors.txt create mode 100644 extra/math/dual/dual-docs.factor create mode 100644 extra/math/dual/dual-tests.factor create mode 100644 extra/math/dual/dual.factor diff --git a/basis/math/functions/functions.factor b/basis/math/functions/functions.factor index 85b4d711ac..3a1ce18daa 100644 --- a/basis/math/functions/functions.factor +++ b/basis/math/functions/functions.factor @@ -252,10 +252,14 @@ M: real tanh ftanh ; : -i* ( x -- y ) >rect swap neg rect> ; -: asin ( x -- y ) +GENERIC: asin ( x -- y ) foldable + +M: number asin dup [-1,1]? [ fasin ] [ i* asinh -i* ] if ; inline -: acos ( x -- y ) +GENERIC: acos ( x -- y ) foldable + +M: number acos dup [-1,1]? [ facos ] [ asin pi 2 / swap - ] if ; inline diff --git a/extra/math/derivatives/authors.txt b/extra/math/derivatives/authors.txt new file mode 100644 index 0000000000..b6089d8622 --- /dev/null +++ b/extra/math/derivatives/authors.txt @@ -0,0 +1 @@ +Jason W. Merrill \ No newline at end of file diff --git a/extra/math/derivatives/derivatives-docs.factor b/extra/math/derivatives/derivatives-docs.factor new file mode 100644 index 0000000000..4905f260bc --- /dev/null +++ b/extra/math/derivatives/derivatives-docs.factor @@ -0,0 +1,11 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: help.markup help.syntax ; +IN: math.derivatives + +ARTICLE: "math.derivatives" "Derivatives" +"The " { $vocab-link "math.derivatives" } " vocabulary defines the derivative of many of the words in the " { $vocab-link "math" } " and " { $vocab-link "math.functions" } " vocabularies. The derivative for a word is given by a sequence of quotations stored in its " { $snippet "derivative" } " word property that give the partial derivative of the word with respect to each of its inputs." +{ $see-also "math.derivatives.syntax" } +; + +ABOUT: "math.derivatives" diff --git a/extra/math/derivatives/derivatives-tests.factor b/extra/math/derivatives/derivatives-tests.factor new file mode 100644 index 0000000000..f95eb43849 --- /dev/null +++ b/extra/math/derivatives/derivatives-tests.factor @@ -0,0 +1,4 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: tools.test automatic-differentiation.derivatives ; +IN: automatic-differentiation.derivatives.tests diff --git a/extra/math/derivatives/derivatives.factor b/extra/math/derivatives/derivatives.factor new file mode 100644 index 0000000000..8e69cec129 --- /dev/null +++ b/extra/math/derivatives/derivatives.factor @@ -0,0 +1,34 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: kernel math math.functions math.derivatives.syntax ; +IN: math.derivatives + +DERIVATIVE: + [ 2drop ] [ 2drop ] +DERIVATIVE: - [ 2drop ] [ 2drop neg ] +DERIVATIVE: * [ nip * ] [ drop * ] +DERIVATIVE: / [ nip / ] [ sq / neg * ] +! Conditional checks if the epsilon-part of the exponent is +! 0 to avoid getting float answers for integer powers. +DERIVATIVE: ^ [ [ 1 - ^ ] keep * * ] + [ [ dup zero? ] 2dip [ 3drop 0 ] [ [ ^ ] keep log * * ] if ] + +DERIVATIVE: sqrt [ sqrt 2 * / ] + +DERIVATIVE: exp [ exp * ] +DERIVATIVE: log [ / ] + +DERIVATIVE: sin [ cos * ] +DERIVATIVE: cos [ sin neg * ] +DERIVATIVE: tan [ sec sq * ] + +DERIVATIVE: sinh [ cosh * ] +DERIVATIVE: cosh [ sinh * ] +DERIVATIVE: tanh [ sech sq * ] + +DERIVATIVE: asin [ sq neg 1 + sqrt / ] +DERIVATIVE: acos [ sq neg 1 + sqrt neg / ] +DERIVATIVE: atan [ sq 1 + / ] + +DERIVATIVE: asinh [ sq 1 + sqrt / ] +DERIVATIVE: acosh [ [ 1 + sqrt ] [ 1 - sqrt ] bi * / ] +DERIVATIVE: atanh [ sq neg 1 + / ] \ No newline at end of file diff --git a/extra/math/derivatives/syntax/authors.txt b/extra/math/derivatives/syntax/authors.txt new file mode 100644 index 0000000000..b6089d8622 --- /dev/null +++ b/extra/math/derivatives/syntax/authors.txt @@ -0,0 +1 @@ +Jason W. Merrill \ No newline at end of file diff --git a/extra/math/derivatives/syntax/syntax-docs.factor b/extra/math/derivatives/syntax/syntax-docs.factor new file mode 100644 index 0000000000..2273e7b83c --- /dev/null +++ b/extra/math/derivatives/syntax/syntax-docs.factor @@ -0,0 +1,18 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: help.markup help.syntax ; +IN: math.derivatives.syntax + +HELP: DERIVATIVE: +{ $description "Defines the derivative of a word by setting its " { $snippet "derivative" } " word property. Reads a word followed by " { $snippet "n" } " quotations, giving the " { $snippet "n" } " partial derivatives of the word with respect to each of its arguments successively. Each quotation should take " { $snippet "n + 1" } " inputs, where the first input is an increment and the last " { $snippet "n" } " inputs are the point at which to evaluate the derivative. The derivative should be a linear function of the increment, and should have the same number of outputs as the original word." } +{ $examples + { $unchecked-example "USING: math math.functions math.derivatives.syntax ;" + "DERIVATIVE: sin [ cos * ]" + "DERIVATIVE: * [ nip * ] [ drop * ]" "" } +} ; + +ARTICLE: "math.derivatives.syntax" "Derivative Syntax" +"The " { $vocab-link "math.derivatives.syntax" } " vocabulary provides the " { $link POSTPONE: DERIVATIVE: } " syntax for specifying the derivative of a word." +; + +ABOUT: "math.derivatives.syntax" diff --git a/extra/math/derivatives/syntax/syntax.factor b/extra/math/derivatives/syntax/syntax.factor new file mode 100644 index 0000000000..02b0608ed8 --- /dev/null +++ b/extra/math/derivatives/syntax/syntax.factor @@ -0,0 +1,10 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: kernel parser words effects accessors sequences + math.ranges ; + +IN: math.derivatives.syntax + +: DERIVATIVE: scan-object dup stack-effect in>> length [1,b] + [ drop scan-object ] map + "derivative" set-word-prop ; parsing \ No newline at end of file diff --git a/extra/math/dual/authors.txt b/extra/math/dual/authors.txt new file mode 100644 index 0000000000..b6089d8622 --- /dev/null +++ b/extra/math/dual/authors.txt @@ -0,0 +1 @@ +Jason W. Merrill \ No newline at end of file diff --git a/extra/math/dual/dual-docs.factor b/extra/math/dual/dual-docs.factor new file mode 100644 index 0000000000..de3b0749a5 --- /dev/null +++ b/extra/math/dual/dual-docs.factor @@ -0,0 +1,90 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: help.markup help.syntax kernel words math math.functions math.derivatives.syntax ; +IN: math.dual + +HELP: +{ $values + { "ordinary-part" real } { "epsilon-part" real } + { "dual" dual number } +} +{ $description "Creates a dual number from its ordinary and epsilon parts." } ; + +HELP: d* +{ $values + { "x" dual } { "y" dual } + { "x*y" dual } +} +{ $description "Multiply dual numbers." } ; + +HELP: d+ +{ $values + { "x" dual } { "y" dual } + { "x+y" dual } +} +{ $description "Add dual numbers." } ; + +HELP: d- +{ $values + { "x" dual } { "y" dual } + { "x-y" dual } +} +{ $description "Subtract dual numbers." } ; + +HELP: d/ +{ $values + { "x" dual } { "y" dual } + { "x/y" dual } +} +{ $description "Divide dual numbers." } +{ $errors "Throws an error if the ordinary part of " { $snippet "x" } " is zero." } ; + +HELP: d^ +{ $values + { "x" dual } { "y" dual } + { "x^y" dual } +} +{ $description "Raise a dual number to a (possibly dual) power" } ; + +HELP: define-dual-method +{ $values + { "word" word } +} +{ $description "Defines a method on the dual numbers for generic word." } +{ $notes "Uses the derivative word-prop, which holds a list of quotations giving the partial derivatives of the word with respect to each of its arguments. This can be set using " { $link POSTPONE: DERIVATIVE: } "." } ; + +{ define-dual-method dual-op POSTPONE: DERIVATIVE: } related-words + +HELP: dual +{ $class-description "The class of dual numbers with non-zero epsilon part." } ; + +HELP: dual-op +{ $values + { "word" word } +} +{ $description "Similar to " { $link execute } ", but promotes word to operate on duals." } +{ $notes "Uses the derivative word-prop, which holds a list of quotations giving the partial derivatives of the word with respect to each of its arguments. This can be set using " { $link POSTPONE: DERIVATIVE: } ". Once a derivative has been defined for a word, dual-op makes it easy to extend the definition to dual numbers." } +{ $examples + { $unchecked-example "USING: math math.dual math.derivatives.syntax math.functions ;" + "DERIVATIVE: sin [ cos * ]" + "M: dual sin \\sin dual-op ;" "" } + { $unchecked-example "USING: math math.dual math.derivatives.syntax ;" + "DERIVATIVE: * [ drop ] [ nip ]" + ": d* ( x y -- x*y ) \ * dual-op ;" "" } +} ; + +HELP: unpack-dual +{ $values + { "dual" dual } + { "ordinary-part" number } { "epsilon-part" number } +} +{ $description "Extracts the ordinary and epsilon part of a dual number." } ; + +ARTICLE: "math.dual" "Dual Numbers" +"The " { $vocab-link "math.dual" } " vocabulary implements dual numbers, along with arithmetic methods for working with them. Many of the functions in " { $vocab-link "math.functions" } " are extended to work with dual numbers." +$nl +"Dual numbers are ordered pairs " { $snippet ""} "--an ordinary part and an epsilon part--with component-wise addition and multiplication defined by "{ $snippet "* = " } ". They are analagous to complex numbers with " { $snippet "i^2 = 0" } "instead of " { $snippet "i^2 = -1" } ". For well-behaved functions " { $snippet "f" } ", " { $snippet "f() = f(o1) + e1*f'(o1)" } ", where " { $snippet "f'"} " is the derivative of " { $snippet "f" } "." +; + + +ABOUT: "math.dual" diff --git a/extra/math/dual/dual-tests.factor b/extra/math/dual/dual-tests.factor new file mode 100644 index 0000000000..2fe751dd63 --- /dev/null +++ b/extra/math/dual/dual-tests.factor @@ -0,0 +1,14 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: tools.test math.dual kernel accessors math math.functions + math.constants ; +IN: math.dual.tests + +[ 0.0 1.0 ] [ 0 1 sin unpack-dual ] unit-test +[ 1.0 0.0 ] [ 0 1 cos unpack-dual ] unit-test +[ 3 5 ] [ 1 5 2 d+ unpack-dual ] unit-test +[ 0 -1 ] [ 1 5 1 6 d- unpack-dual ] unit-test +[ 2 1 ] [ 2 3 1 -1 d* unpack-dual ] unit-test +[ 1/2 -1/4 ] [ 2 1 1 swap d/ unpack-dual ] unit-test +[ 2 ] [ 1 1 2 d^ epsilon-part>> ] unit-test +[ 2.0 .25 ] [ 4 1 sqrt unpack-dual ] unit-test \ No newline at end of file diff --git a/extra/math/dual/dual.factor b/extra/math/dual/dual.factor new file mode 100644 index 0000000000..214db9b678 --- /dev/null +++ b/extra/math/dual/dual.factor @@ -0,0 +1,80 @@ +! Copyright (C) 2009 Jason W. Merrill. +! See http://factorcode.org/license.txt for BSD license. +USING: kernel math math.functions math.derivatives accessors + macros words effects sequences generalizations fry + combinators.smart generic compiler.units ; + +IN: math.dual + +TUPLE: dual ordinary-part epsilon-part ; + +C: dual + +! Ordinary numbers implement the dual protocol by returning +! themselves as the ordinary part, and 0 as the epsilon part. +M: number ordinary-part>> ; + +M: number epsilon-part>> drop 0 ; + +: unpack-dual ( dual -- ordinary-part epsilon-part ) + [ ordinary-part>> ] [ epsilon-part>> ] bi ; + +> length ; + +MACRO: ordinary-op ( word -- o ) + [ input-length ] keep + '[ [ ordinary-part>> ] _ napply _ execute ] ; + +! Takes N dual numbers ... and weaves +! their ordinary and epsilon parts to produce +! e1 o1 o2 ... oN e2 o1 o2 ... oN ... eN o1 o2 ... oN +! This allows a set of partial derivatives each to be evaluated +! at the same point. +MACRO: duals>nweave ( n -- ) + dup dup dup + '[ + [ [ epsilon-part>> ] _ napply ] + _ nkeep + [ ordinary-part>> ] _ napply + _ nweave + ] ; + +MACRO: chain-rule ( word -- e ) + [ input-length '[ _ duals>nweave ] ] + [ "derivative" word-prop ] + [ input-length 1+ '[ _ nspread ] ] + tri + '[ [ @ _ @ ] sum-outputs ] ; + +PRIVATE> + +MACRO: dual-op ( word -- ) + [ '[ _ ordinary-op ] ] + [ input-length '[ _ nkeep ] ] + [ '[ _ chain-rule ] ] + tri + '[ _ @ @ ] ; + +: define-dual-method ( word -- ) + [ \ dual swap create-method ] keep '[ _ dual-op ] define ; + +! Specialize math functions to operate on dual numbers. +[ { sqrt exp log sin cos tan sinh cosh tanh acos asin atan } + [ define-dual-method ] each ] with-compilation-unit + +! Inverse methods { asinh, acosh, atanh } are not generic, so +! there is no way to specialize them for dual numbers. However, +! they are defined in terms of functions that can operate on +! dual numbers and arithmetic methods, so if it becomes +! possible to make arithmetic operators work directly on dual +! numbers, we will get these for free. + +! Arithmetic methods are not generic (yet?), so we have to +! define special versions of them to operate on dual numbers. +: d+ ( x y -- x+y ) \ + dual-op ; +: d- ( x y -- x-y ) \ - dual-op ; +: d* ( x y -- x*y ) \ * dual-op ; +: d/ ( x y -- x/y ) \ / dual-op ; +: d^ ( x y -- x^y ) \ ^ dual-op ; \ No newline at end of file From 8968093623a2bc7087527fcce6cc7c324148e3af Mon Sep 17 00:00:00 2001 From: Jason Merrill Date: Wed, 18 Feb 2009 21:28:48 -0500 Subject: [PATCH 2/2] Added dual versions of a few more words to math.dual. --- extra/math/derivatives/derivatives.factor | 23 +++++++++++-- extra/math/dual/dual-docs.factor | 42 +++++++++++++++++++++++ extra/math/dual/dual-tests.factor | 4 ++- extra/math/dual/dual.factor | 26 ++++++++++---- 4 files changed, 85 insertions(+), 10 deletions(-) diff --git a/extra/math/derivatives/derivatives.factor b/extra/math/derivatives/derivatives.factor index 8e69cec129..c6a9d1a357 100644 --- a/extra/math/derivatives/derivatives.factor +++ b/extra/math/derivatives/derivatives.factor @@ -1,8 +1,15 @@ ! Copyright (C) 2009 Jason W. Merrill. ! See http://factorcode.org/license.txt for BSD license. -USING: kernel math math.functions math.derivatives.syntax ; +USING: kernel math math.functions math.derivatives.syntax + math.order math.parser summary accessors make combinators ; IN: math.derivatives +ERROR: undefined-derivative point word ; +M: undefined-derivative summary + [ dup "Derivative of " % word>> name>> % + " is undefined at " % point>> # "." % ] + "" make ; + DERIVATIVE: + [ 2drop ] [ 2drop ] DERIVATIVE: - [ 2drop ] [ 2drop neg ] DERIVATIVE: * [ nip * ] [ drop * ] @@ -12,6 +19,15 @@ DERIVATIVE: / [ nip / ] [ sq / neg * ] DERIVATIVE: ^ [ [ 1 - ^ ] keep * * ] [ [ dup zero? ] 2dip [ 3drop 0 ] [ [ ^ ] keep log * * ] if ] +DERIVATIVE: abs + [ 0 <=> + { + { +lt+ [ neg ] } + { +eq+ [ 0 \ abs undefined-derivative ] } + { +gt+ [ ] } + } case + ] + DERIVATIVE: sqrt [ sqrt 2 * / ] DERIVATIVE: exp [ exp * ] @@ -31,4 +47,7 @@ DERIVATIVE: atan [ sq 1 + / ] DERIVATIVE: asinh [ sq 1 + sqrt / ] DERIVATIVE: acosh [ [ 1 + sqrt ] [ 1 - sqrt ] bi * / ] -DERIVATIVE: atanh [ sq neg 1 + / ] \ No newline at end of file +DERIVATIVE: atanh [ sq neg 1 + / ] + +DERIVATIVE: neg [ drop neg ] +DERIVATIVE: recip [ sq recip neg * ] diff --git a/extra/math/dual/dual-docs.factor b/extra/math/dual/dual-docs.factor index de3b0749a5..6c287a8f1e 100644 --- a/extra/math/dual/dual-docs.factor +++ b/extra/math/dual/dual-docs.factor @@ -46,6 +46,48 @@ HELP: d^ } { $description "Raise a dual number to a (possibly dual) power" } ; +HELP: dabs +{ $values + { "x" dual } + { "|x|" dual } +} +{ $description "Absolute value of a dual number." } ; + +HELP: dacosh +{ $values + { "x" dual } + { "y" dual } +} +{ $description "Inverse hyberbolic cosine of a dual number." } ; + +HELP: dasinh +{ $values + { "x" dual } + { "y" dual } +} +{ $description "Inverse hyberbolic sine of a dual number." } ; + +HELP: datanh +{ $values + { "x" dual } + { "y" dual } +} +{ $description "Inverse hyberbolic tangent of a dual number." } ; + +HELP: dneg +{ $values + { "x" dual } + { "-x" dual } +} +{ $description "Negative of a dual number." } ; + +HELP: drecip +{ $values + { "x" dual } + { "1/x" dual } +} +{ $description "Reciprocal of a dual number." } ; + HELP: define-dual-method { $values { "word" word } diff --git a/extra/math/dual/dual-tests.factor b/extra/math/dual/dual-tests.factor index 2fe751dd63..ea46c46124 100644 --- a/extra/math/dual/dual-tests.factor +++ b/extra/math/dual/dual-tests.factor @@ -11,4 +11,6 @@ IN: math.dual.tests [ 2 1 ] [ 2 3 1 -1 d* unpack-dual ] unit-test [ 1/2 -1/4 ] [ 2 1 1 swap d/ unpack-dual ] unit-test [ 2 ] [ 1 1 2 d^ epsilon-part>> ] unit-test -[ 2.0 .25 ] [ 4 1 sqrt unpack-dual ] unit-test \ No newline at end of file +[ 2.0 .25 ] [ 4 1 sqrt unpack-dual ] unit-test +[ 2 -1 ] [ -2 1 dabs unpack-dual ] unit-test +[ -2 -1 ] [ 2 1 dneg unpack-dual ] unit-test \ No newline at end of file diff --git a/extra/math/dual/dual.factor b/extra/math/dual/dual.factor index 214db9b678..36d684bc6d 100644 --- a/extra/math/dual/dual.factor +++ b/extra/math/dual/dual.factor @@ -51,9 +51,9 @@ MACRO: chain-rule ( word -- e ) PRIVATE> MACRO: dual-op ( word -- ) - [ '[ _ ordinary-op ] ] - [ input-length '[ _ nkeep ] ] - [ '[ _ chain-rule ] ] + [ '[ _ ordinary-op ] ] + [ input-length '[ _ nkeep ] ] + [ '[ _ chain-rule ] ] tri '[ _ @ @ ] ; @@ -64,17 +64,29 @@ MACRO: dual-op ( word -- ) [ { sqrt exp log sin cos tan sinh cosh tanh acos asin atan } [ define-dual-method ] each ] with-compilation-unit -! Inverse methods { asinh, acosh, atanh } are not generic, so +! Inverse methods { asinh, acosh, atanh } are not generic, so ! there is no way to specialize them for dual numbers. However, ! they are defined in terms of functions that can operate on ! dual numbers and arithmetic methods, so if it becomes ! possible to make arithmetic operators work directly on dual ! numbers, we will get these for free. -! Arithmetic methods are not generic (yet?), so we have to +! Arithmetic words are not generic (yet?), so we have to ! define special versions of them to operate on dual numbers. : d+ ( x y -- x+y ) \ + dual-op ; -: d- ( x y -- x-y ) \ - dual-op ; +: d- ( x y -- x-y ) \ - dual-op ; : d* ( x y -- x*y ) \ * dual-op ; : d/ ( x y -- x/y ) \ / dual-op ; -: d^ ( x y -- x^y ) \ ^ dual-op ; \ No newline at end of file +: d^ ( x y -- x^y ) \ ^ dual-op ; + +: dabs ( x -- |x| ) \ abs dual-op ; + +! The following words are also not generic, but are defined in +! terms of words that can operate on dual numbers and +! arithmetic. If it becomes possible to implement arithmetic on +! dual numbers directly, these functions can be deleted. +: dneg ( x -- -x ) \ neg dual-op ; +: drecip ( x -- 1/x ) \ recip dual-op ; +: dasinh ( x -- y ) \ asinh dual-op ; +: dacosh ( x -- y ) \ acosh dual-op ; +: datanh ( x -- y ) \ atanh dual-op ; \ No newline at end of file