! 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 combinators macros make ;
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 ;
> ] n napply ] [ [ ordinary-part>> ] n napply ] }
! n ncleave
! derivative-list [ n ncurry ] n nwith map
! spread n narray sum ; inline
MACRO: weave ( n -- quot )
[ '[ [ ] _ ncurry ] ]
[ [ '[ _ ndip ] ] map ]
bi
'[ @ _ cleave ] ;
MACRO: nspread ( quots n -- quot )
over empty? [ 2drop [ ] ] [
[ [ but-last ] dip ] [ [ peek ] dip ] 2bi swap
'[ [ _ _ nspread ] _ ndip @ ]
] if ;
:: chain-rule ( derivative-list n -- x )
[ [ epsilon-part>> ] n napply ] n nkeep [ ordinary-part>> ] n napply
n weave
derivative-list n 1+ nspread
n narray sum ; inline
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.
! :: [dual-op] ( word -- quot )
! word "derivative" word-prop :> derivative-list
! derivative-list length :> n
! [ [ [ ordinary-part>> ] n napply word execute ] n nkeep
! derivative-list n chain-rule
! ] ;
:: [dual-op] ( word -- quot )
word "derivative" word-prop :> derivative-list
derivative-list length :> n
n word '[ [ ordinary-part>> ] _ napply _ execute ] n
derivative-list n
'[
_ _ nkeep
_ _ chain-rule
] ;
MACRO: dual-op ( word -- ) [dual-op] ;
: define-dual-method ( word -- )
[ \ dual swap create-method-in ] 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 >>
! 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 ;