1 ! Copyright (C) 2009 Jason W. Merrill.
2 ! See http://factorcode.org/license.txt for BSD license.
3 USING: kernel math math.functions math.derivatives accessors
4 macros words effects sequences generalizations fry
5 combinators.smart generic compiler.units ;
9 TUPLE: dual ordinary-part epsilon-part ;
13 ! Ordinary numbers implement the dual protocol by returning
14 ! themselves as the ordinary part, and 0 as the epsilon part.
15 M: number ordinary-part>> ;
17 M: number epsilon-part>> drop 0 ;
19 : unpack-dual ( dual -- ordinary-part epsilon-part )
20 [ ordinary-part>> ] [ epsilon-part>> ] bi ;
24 : input-length ( word -- n ) stack-effect in>> length ;
26 MACRO: ordinary-op ( word -- o )
28 '[ [ ordinary-part>> ] _ napply _ execute ] ;
30 ! Takes N dual numbers <o1,e1> <o2,e2> ... <oN,eN> and weaves
31 ! their ordinary and epsilon parts to produce
32 ! e1 o1 o2 ... oN e2 o1 o2 ... oN ... eN o1 o2 ... oN
33 ! This allows a set of partial derivatives each to be evaluated
35 MACRO: duals>nweave ( n -- )
38 [ [ epsilon-part>> ] _ napply ]
40 [ ordinary-part>> ] _ napply
44 MACRO: chain-rule ( word -- e )
45 [ input-length '[ _ duals>nweave ] ]
46 [ "derivative" word-prop ]
47 [ input-length 1+ '[ _ nspread ] ]
49 '[ [ @ _ @ ] sum-outputs ] ;
53 MACRO: dual-op ( word -- )
54 [ '[ _ ordinary-op ] ]
55 [ input-length '[ _ nkeep ] ]
60 : define-dual-method ( word -- )
61 [ \ dual swap create-method ] keep '[ _ dual-op ] define ;
63 ! Specialize math functions to operate on dual numbers.
64 [ { sqrt exp log sin cos tan sinh cosh tanh acos asin atan }
65 [ define-dual-method ] each ] with-compilation-unit
67 ! Inverse methods { asinh, acosh, atanh } are not generic, so
68 ! there is no way to specialize them for dual numbers. However,
69 ! they are defined in terms of functions that can operate on
70 ! dual numbers and arithmetic methods, so if it becomes
71 ! possible to make arithmetic operators work directly on dual
72 ! numbers, we will get these for free.
74 ! Arithmetic words are not generic (yet?), so we have to
75 ! define special versions of them to operate on dual numbers.
76 : d+ ( x y -- x+y ) \ + dual-op ;
77 : d- ( x y -- x-y ) \ - dual-op ;
78 : d* ( x y -- x*y ) \ * dual-op ;
79 : d/ ( x y -- x/y ) \ / dual-op ;
80 : d^ ( x y -- x^y ) \ ^ dual-op ;
82 : dabs ( x -- |x| ) \ abs dual-op ;
84 ! The following words are also not generic, but are defined in
85 ! terms of words that can operate on dual numbers and
86 ! arithmetic. If it becomes possible to implement arithmetic on
87 ! dual numbers directly, these functions can be deleted.
88 : dneg ( x -- -x ) \ neg dual-op ;
89 : drecip ( x -- 1/x ) \ recip dual-op ;
90 : dasinh ( x -- y ) \ asinh dual-op ;
91 : dacosh ( x -- y ) \ acosh dual-op ;
92 : datanh ( x -- y ) \ atanh dual-op ;