extra/math/dual/dual.factor

   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 ;
   6
   7 IN: math.dual
   8
   9 TUPLE: dual ordinary-part epsilon-part ;
  10
  11 C: <dual> dual
  12
  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>> ;
  16
  17 M: number epsilon-part>> drop 0 ;
  18
  19 : unpack-dual ( dual -- ordinary-part epsilon-part )
  20     [ ordinary-part>> ] [ epsilon-part>> ] bi ;
  21
  22 <PRIVATE
  23
  24 : input-length ( word -- n ) stack-effect in>> length ;
  25
  26 MACRO: ordinary-op ( word -- o )
  27     [ input-length ] keep
  28     '[ [ ordinary-part>> ] _ napply _ execute ] ;
  29
  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
  34 ! at the same point.
  35 MACRO: duals>nweave ( n -- )
  36    dup dup dup
  37    '[
  38        [ [ epsilon-part>> ] _ napply ]
  39        _ nkeep
  40        [ ordinary-part>> ] _ napply
  41        _ nweave
  42     ] ;
  43
  44 MACRO: chain-rule ( word -- e )
  45     [ input-length '[ _ duals>nweave ] ]
  46     [ "derivative" word-prop ]
  47     [ input-length 1+ '[ _ nspread ] ]
  48     tri
  49     '[ [ @ _ @ ] sum-outputs ] ;
  50
  51 PRIVATE>
  52
  53 MACRO: dual-op ( word -- )
  54     [ '[ _ ordinary-op ] ]
  55     [ input-length '[ _ nkeep ] ]
  56     [ '[ _ chain-rule ] ]
  57     tri
  58     '[ _ @ @ <dual> ] ;
  59
  60 : define-dual-method ( word -- )
  61     [ \ dual swap create-method ] keep '[ _ dual-op ] define ;
  62
  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
  66
  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.
  73
  74 ! Arithmetic methods 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 ;