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.db4
Jason Merrill
parent
489bb32a98
commit
fe55e939f9
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1 @@
|
|||
Jason W. Merrill
|
||||
|
|
@ -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"
|
||||
|
|
@ -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
|
||||
|
|
@ -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 + / ]
|
||||
|
|
@ -0,0 +1 @@
|
|||
Jason W. Merrill
|
||||
|
|
@ -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"
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1 @@
|
|||
Jason W. Merrill
|
||||
|
|
@ -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: <dual>
|
||||
{ $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 "<o,e>"} "--an ordinary part and an epsilon part--with component-wise addition and multiplication defined by "{ $snippet "<o1,e1>*<o2,e2> = <o1*o2,e1*o2 + e2*o1>" } ". 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(<o1,e1>) = f(o1) + e1*f'(o1)" } ", where " { $snippet "f'"} " is the derivative of " { $snippet "f" } "."
|
||||
;
|
||||
|
||||
|
||||
ABOUT: "math.dual"
|
||||
|
|
@ -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 <dual> sin unpack-dual ] unit-test
|
||||
[ 1.0 0.0 ] [ 0 1 <dual> cos unpack-dual ] unit-test
|
||||
[ 3 5 ] [ 1 5 <dual> 2 d+ unpack-dual ] unit-test
|
||||
[ 0 -1 ] [ 1 5 <dual> 1 6 <dual> d- unpack-dual ] unit-test
|
||||
[ 2 1 ] [ 2 3 <dual> 1 -1 <dual> d* unpack-dual ] unit-test
|
||||
[ 1/2 -1/4 ] [ 2 1 <dual> 1 swap d/ unpack-dual ] unit-test
|
||||
[ 2 ] [ 1 1 <dual> 2 d^ epsilon-part>> ] unit-test
|
||||
[ 2.0 .25 ] [ 4 1 <dual> sqrt unpack-dual ] unit-test
|
||||
|
|
@ -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> 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 ;
|
||||
|
||||
<PRIVATE
|
||||
|
||||
: input-length ( word -- n ) stack-effect in>> length ;
|
||||
|
||||
MACRO: ordinary-op ( word -- o )
|
||||
[ input-length ] keep
|
||||
'[ [ ordinary-part>> ] _ napply _ execute ] ;
|
||||
|
||||
! Takes N dual numbers <o1,e1> <o2,e2> ... <oN,eN> 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
|
||||
'[ _ @ @ <dual> ] ;
|
||||
|
||||
: 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 ;
|
||||
Loading…
Reference in New Issue