Re: Differentiable Forth

Hi,
Here is a program that gives examples of using dual numbers to calculate
derivatives
(it is not optimized,..., just a proof of concept)
The code begins here
: -frot frot frot ;
: f>dl   ( f: a -- a 0) 0e ;
: dl>rl  ( f: a b -- a) fdrop ;
: dl>eps ( f: a b -- b) fnip ;
: dl. ( f: a b -- ) fswap f. ." + eps " f. ;
: dl(,). ( f: a b -- ) ." (" fswap f. ." , " f. ." )" ;
: dl+ ( f: a b c d -- a+c b+d)
    frot ( f: a c d b)
    f+   ( f: a c b+d)
    -frot ( f: b+d a c)
    f+   ( f: b+d a+c)
    fswap ( f: a+c b+d) ;
: dl- ( f: a b c d -- a-c b-d)
    frot ( f: a c d b)
    fswap f-   ( f: a c b-d)
    -frot ( f: b-d a c)
    f-   ( f: b-d a-c)
    fswap ( f: a-c b-d) ;
fvariable b*c
: dl* ( f: a b c d -- a*c  a*d+b*c)
    -frot ( f: a d b c)
    ftuck f* ( f: a d c b*c)
    b*c f! ( f: a d c)
    frot ( f: d c a)
    ftuck ( f: d a c a)
    f*    ( f: d a a*c)
    -frot  ( f: a*c d a)
    f* b*c f@ f+ ( f: a*c a*d+b*c)
;
: 1/dl ( f: a b -- 1/a -b*1/a^2)
    fswap 1/f ( f: b 1/a)
    ftuck     ( f: 1/a b 1/a)
    fdup f*   ( f: 1/a b 1/a^2)
    f* fnegate ( f: 1/a -b/a^2)
;
: dl/ ( f: a b c d -- a/c rest)
    1/dl dl*
;
: dl^f ( f: a b c -- a^c b*c)
    ftuck ( f: a c b c)
    f* -frot ( f: b*c a c)
    f** fswap ( f: a^c b*c)
;
\
: dldup ( f: a b -- a b a b) fover fover ;
: dlnip ( f: a b c d -- c d) frot fdrop frot fdrop ;
: dldrop ( f: a b -- ) fdrop fdrop ;
fvariable dlswap_temp
: dlswap ( f: a b c d -- c d a b)
    frot dlswap_temp f! ( f: a c d)
    frot ( f: c d a)
    dlswap_temp f@  ;
fvariable dlover_temp1
fvariable dlover_temp2
: dlover ( f: a b c d -- a b c d a b)
    dlswap ( f: a b c d -- c d a b)
    dlover_temp2 f! dlover_temp1 f! ( f: c d)
    dlover_temp1 f@ dlover_temp2 f@ ( f: c d a b)
    dlswap ( f: a b c d)
    dlover_temp1 f@ dlover_temp2 f@ ( f: a b c d a b)
;
: dltuck dlswap dlover ;
\ dual number funuctions of dula number variables
: dlvariable create 2 floats allot ;
: dl! ( dlvar -- ) ( f: f1 f2 -- ) dup float+ f! f! ;
: dl@ ( dlvar -- ) ( f: -- f1 f2)  dup f@ float+ f@ ;
: dlident ( f: a b -- a b) ;
: dlsin ( f: a b -- c d) fswap fdup fsin fswap fcos frot f* ;
: dlcos ( f: a b -- c d) fswap fdup fcos fswap fsin fnegate frot f* ;
: dlexp ( f: a b -- c d) fswap fdup fexp fswap fexp frot f* ;
: dlln ( f: a b -- c d)  fswap fdup fln  fswap 1/f frot f* ;
\ ..... add others
\ derivatives
variable func
: der 1e ' func ! func @ execute dl>eps ;
\ examples
1 [if]
: dlf() 1e 0e dl+ ; \ f(x) = x + 1
: dlf2() dldup dl* ; \ f2(x) = x^2
: dlf3() dldup dlf2() dlswap dlf() dl+ ; \ f3(x) = x^2 + x + 1
: der_f3() ( f: x -- y) 2e f* 1e f+ ; \ d/dx(f3) = 2*x + 1
cr 1e der dlf3() f. \ 3. ok
cr 1e der_f3() f.   \ 3. ok
cr cr
: dlf4() dlf3() dlexp ; \ f3(x) = exp(x^2+x+1)
: der_f4() ( f: x -- y) \ d/dx(f4) = (2*x+1)*exp(x^2+x+1)
    fdup 2e f* 1e f+ fswap
    fdup fdup f* f+ 1e f+ fexp f*
;
cr 2e der dlf4() f. \ 5483.16579214229  ok, calculated at 2e
cr 2e der_f4() f.   \ 5483.16579214229  ok
cr cr
: dlf5() 1/dl ;
: der_f5() ( f: x) \ d/dx(f5) = -1/x^2
    fdup f* 1/f fnegate
;
cr 2e der dlf5() f. \ -0.25  ok calculated at 2e
cr 2e der_f5() f.   \ -0.25
cr cr
: dlf6() dldup dldup dl* dlswap dlsin dl+ 1/dl ; \ f6(x) =
1/(x^2+sin(x))
\ using the derivative calculated analytically d/dx (f6(x)) =
-(2*x+cos(x))/(x^2+sin(x))^2
: der_f6() ( f: x -- y) fdup fdup fdup f* fswap fsin f+ fdup f* 1/f
fswap fdup 2e f* fswap fcos f+ f* fnegate ;
cr 1e der dlf6() f. \ -0.749127330692909  ok calculated at x = 1
cr 1e der_f6() f.   \ -0.749127330692909  ok
cr cr
: dlf7() dldup dldup dl* dlcos dl* ; \ f7(x) = x*cos(x^2),
: der_f7() ( f: x --y)  \ its deriv: d/dx(f7) = x*cos(x^2)-2*x*sin(x^2)
    fdup fdup f* ( f: x x^2)
    fsincos ( f: x s c )
    frot ( f: s c x)
    ftuck ( f: s x c x)
    f* ( f: s x cx)
    -frot ( f: cx s x)
    f* 2e f* f-
;
cr 1e der_f7() f.   \ -1.14263966374765  ok calculated at 1e
cr 1e der dlf7() f. \ -1.14263966374765  ok
cr cr cr
[then]
The code ends here.
Ahmed