USING: kernel math math.functions math.derivatives accessors
words generalizations sequences generic fry locals
compiler.units continuations quotations combinators
combinators.smart macros effects ;
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 -- )
[ 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>weave ( n -- )
dup dup dup '[
[ [ epsilon-part>> ] _ napply ]
_ nkeep
[ ordinary-part>> ] _ napply
_ nweave
] ;
MACRO: chain-rule ( word -- x )
[ input-length ]
[ "derivative" word-prop ]
[ input-length 1+ ]
tri
'[ [ _ duals>weave _ _ nspread ] sum-outputs ] ;
PRIVATE>
! dual-op is similar to execute, but promotes the word to
! operate on duals. Uses the "derivative" word-prop, which
! holds a list of quotations giving the partial derivative of
! the word with respect to each of its argumets.
MACRO: dual-op ( word -- )
[ ] [ input-length ] [ ] tri
'[ [ _ ordinary-op ] _ nkeep _ chain-rule ] ;
: 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 atan }
[ define-dual-method ] each ] with-compilation-unit
! Inverse methods { asin, acos, asinh, acosh, atanh } are not
! generic, so there is no way to specialize them for dual
! numbers.
! 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 ;