! Copyright (C) 2009 Jason Merrill.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel math math.functions math.derivatives accessors words
generalizations sequences generic.parser fry locals compiler.units
continuations quotations ;
IN: math.dual
TUPLE: dual ordinary-part epsilon-part ;
C: dual
M: number ordinary-part>> ;
M: number epsilon-part>> drop 0 ;
: unpack-dual ( dual -- ordinary-part epsilon-part )
[ ordinary-part>> ] [ epsilon-part>> ] bi ;
;
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.
! Examines the length of the derivative list (which is equal to the number of
! inputs the word takes), grabs that many duals off the stack, and splits them
! into their ordinary and epsilon parts to send off to (dual-op), which is a
! straightforward implementation of the definition of how functions extend to
! dual numbers.
: dual-op ( duals word -- dual )
dup "derivative" word-prop dup length
'[ _ narray [ [ epsilon-part>> ] map ] [ [ ordinary-part>> ] map ] bi
] 2dip (dual-op) ; inline
: define-dual ( word -- )
[ \ dual swap create-method-in ] keep [ dual-op ] curry define ;
! Specialize math functions to operate on dual numbers.
<< { sqrt exp log sin cos tan sinh cosh tanh atan } [ define-dual ] each >>
! 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 ;