Author: | Jason Merrill |
---|---|

Mode: | factor |

Date: | Thu, 12 Feb 2009 00:49:26 |

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> 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 MACRO: ordinary-op ( word n -- ) swap '[ [ ordinary-part>> ] _ napply _ execute ] ; MACRO:: chain-rule ( derivative-list n -- x ) n n n n derivative-list n 1+ n '[ [ [ epsilon-part>> ] _ napply ] _ nkeep [ ordinary-part>> ] _ napply _ nweave _ _ nspread _ nsum ] ; 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 word n '[ _ _ ordinary-op ] n derivative-list n '[ _ _ chain-rule ] '[ _ _ nkeep @ <dual> ] ; 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 ;

Author: | Jason Merrill |
---|---|

Mode: | factor |

Date: | Thu, 12 Feb 2009 03:09:39 |

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> 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 -- ) [ 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 <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 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 ;

Summary: | |
---|---|

Author: | |

Mode: | |

Body: | |