Paste: dual numbers

Author: Jason Merrill
Mode: factor
Date: Fri, 6 Feb 2009 00:16:35
Plain Text | 
USING: kernel math math.functions math.derivatives accessors words 
    generalizations sequences generic.parser fry locals compiler.units
    continuations quotations combinators macros ;

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    

:: chain-rule ( derivative-list n -- x )
    { [ [ epsilon-part>> ] n napply ] [ [ ordinary-part>> ] n napply ] }
    n ncleave
    derivative-list [ n ncurry ] n nwith map 
    spread 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> ] ;
    
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: fry

Author: slava
Mode: factor
Date: Fri, 6 Feb 2009 00:28:47
Plain Text | 
:: [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
        <dual>
    ] ;

New Annotation

Summary:
Author:
Mode:
Body: