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 : my-join ( seq glue -- seq ) [ '[ [ _ % ] [ , ] interleave ] ] keep make ; : [derivative-at-point] ( derivative-list -- quot ) [ ] [ length '[ _ nkeep ] ] bi my-join [ call ] compose ; ! MACRO: derivative-at-point ( xincs vals derivative-list -- yincs ) ! [derivative-at-point] ; :: chain-rule ( derivative-list n -- x ) [ [ epsilon-part>> ] n napply ] n nkeep [ ordinary-part>> ] n napply derivative-list [derivative-at-point] call 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 ;