Paste: dual numbers

Author: Jason Merrill
Mode: factor
Date: Thu, 5 Feb 2009 01:01:24
Plain Text | 
! Copyright (C) 2009 Jason Merrill.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel math math.functions math.derivatives accessors words 
    generalizations sequences generic.parser fry locals compiler.units
    continuations quotations ;

IN: math.dual

TUPLE: dual ordinary-part epsilon-part ;

C: <dual> dual

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 ( epsilon-list ordinary-list derivative-list -- d )
    epsilon-list derivative-list
    [ '[ ordinary-list _ prefix _ with-datastack first ] call ] 
    2map sum ; inline
    
:: (dual-op) ( epsilon-list ordinary-list word derivative-list -- dual )
    ordinary-list word 1quotation with-datastack first
    epsilon-list ordinary-list derivative-list chain-rule
    <dual> ;    

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.  
! Examines the length of the derivative list (which is equal to the number of 
! inputs the word takes), grabs that many duals off the stack, and splits them
! into their ordinary and epsilon parts to send off to (dual-op), which is a
! straightforward implementation of the definition of how functions extend to
! dual numbers.
: dual-op ( duals word -- dual )
    dup "derivative" word-prop dup length 
    '[ _ narray [ [ epsilon-part>> ] map ] [ [ ordinary-part>> ] map ] bi 
    ] 2dip (dual-op) ; inline

: define-dual ( word -- )
    [ \ dual swap create-method-in ] keep [ dual-op ] curry define ;

! Specialize math functions to operate on dual numbers.
<< { sqrt exp log sin cos tan sinh cosh tanh atan } [ define-dual ] 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.derivatives

Author: Jason Merrill
Mode: factor
Date: Thu, 5 Feb 2009 01:02:43
Plain Text | 
! Copyright (C) 2009 Jason Merrill.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel math math.functions math.derivatives.parser ;
IN: math.derivatives

DERIVATIVE: + [ 2drop ] [ 2drop ]
DERIVATIVE: - [ 2drop ] [ 2drop neg ]
DERIVATIVE: * [ nip * ] [ drop * ]
DERIVATIVE: / [ nip / ] [ sq / neg * ]
! Conditional checks if the epsilon-part of the exponent is 0 to avoid getting
! float answers for integer powers.
DERIVATIVE: ^ [ [ 1 - ^ ] keep * * ] 
    [ [ dup zero? ] 2dip [ 3drop 0 ] [ [ ^ ] keep log * * ] if ]

DERIVATIVE: sqrt [ 2 * / ]

DERIVATIVE: exp [ exp * ]
DERIVATIVE: log [ / ]

DERIVATIVE: sin [ cos * ]
DERIVATIVE: cos [ sin neg * ]
DERIVATIVE: tan [ sec sq * ]

DERIVATIVE: sinh [ cosh * ]
DERIVATIVE: cosh [ sinh * ]
DERIVATIVE: tanh [ sech sq * ]

DERIVATIVE: asin [ sq neg 1 + sqrt / ]
DERIVATIVE: acos [ sq neg 1 + sqrt neg / ]
DERIVATIVE: atan [ sq 1 + / ]

DERIVATIVE: asinh [ sq 1 + sqrt / ]
DERIVATIVE: acosh [ [ 1 + sqrt ] [ 1 - sqrt ] bi * / ]
DERIVATIVE: atanh [ sq neg 1 + / ]

Annotation: math.derivatives.parser

Author: Jason Merrill
Mode: factor
Date: Thu, 5 Feb 2009 01:03:34
Plain Text | 
! Copyright (C) 2009 Jason Merrill.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel parser words effects accessors sequences math.ranges ;
    
IN: math.derivatives.parser

: DERIVATIVE: scan-object dup stack-effect in>> length [1,b] 
    [ drop scan-object ] map "derivative" set-word-prop ; parsing

New Annotation

Summary:
Author:
Mode:
Body: