Paste: math.dual

Author:	Jason Merrill
Mode:	factor
Date:	Thu, 12 Feb 2009 00:49:26

Plain Text |

Delete Paste

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 ;

Annotation: math.duals

Author:	Jason Merrill
Mode:	factor
Date:	Thu, 12 Feb 2009 03:09:39

Plain Text |

Delete Annotation

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 ;

New Annotation

Summary:
Author:
Mode:	textPostgreSQLPowerCenter Parameter Fileactionscriptadaada95ans-forth94antantlrantlr4apacheconfapdlapplescriptaspaspect-jassembly-agcassembly-agsassembly-arm32assembly-m68kassembly-macro32assembly-mcs51assembly-parrotassembly-r2000assembly-x86avroawkbbatchbbjbcelbeanshellbibtexbinsource-agccc#c++cfscriptchillcilclipsclojurecmakecobolcoffeescriptcoldfusioncplex-lpcsscsvcvs-commitddartdjangodockerfiledotdoxygendrawj2ddsssleiffelembperlerlangfactorfhtmlforthfortranfortran90foxprofreemarkergcbasicgettextgherkingnuplotgogradlegroovygsphaskellhaxehexhl7v2hlslhtaccesshtmlhxmli4glicalendariconidlinforminiinno-setupinterlisiojamonjavajava module-infojavaccjavafxjavascriptjcljedit-actionsjflexjhtmljmkjsonjspkotlinlatexlexlilypondlispliterate-haskelllogslogtalklotosluam4macroschedulermailmakefilemaplemarkdownmavenmavscriptmetapostmlmodula3moinmqscmustache-templatemxmlmyghtymysqln3netrexxnqcnsis2objective-cobjectrexxoccamomnimarkopensipsoutlinepascalpatchperlphppikepl-sqlpl-sql9pl1plaintexplantumlpop11postscriptpovraypowerdynamopowershellprogressprologpropertiesprotobufpspptlpurepvwavepyrexpythonqdocrakefilercprdrebolredcoderelax-ng-compactrenderman-ribrestrfcrhtmlroffrpm-specrtfrubyrustrviews#s+sassbtscalaschemesdl/prsgmlshellscriptshtmlsipsippslateslaxsmalltalksmartysmi-mibsparqlsql-loadersqrsquidconfstatasvn-commitswiftswigtcltemplate-toolkittextexinfotldtransact-sqltsptwikityposcripturluscriptvalavbscriptvelocityverilogvhdlvisualbasicvrml2wellknowntextxmlxqxslyabyamlzpt
Body:

Done