diff --git a/example/fsl/README.txt b/example/fsl/README.txt new file mode 100644 index 0000000..06a943a --- /dev/null +++ b/example/fsl/README.txt @@ -0,0 +1,71 @@ +# Forth Scientific Library + +This is a (partial) port of the Forth Scientific Library to RETRO. +It doesn't provide everything at this point, and may not in the +future. I am slowly porting the files as time allows, but this is +not a high priority project for me. + +The original files can be found at: https://www.taygeta.com/fsl/scilib.html + +I claim no copyright on the ported code; see the individual files +for copyright and license conditions. + +## Introduction + +The Forth Scientific Library was created by Skip Carter with the +help of volunteer contributers and reviewers of code. It provides +a collection of routines that may be useful to those doing +scientific or numerically intensive work using Forth. Below is +Skip's original announcement and proposal for the Project. + +Currently the Project's activities are being coordinated by +Charles Montgomery, who would welcome comments, suggestions, code +contributions, or whatever, at cgm@physics.utoledo.edu + + +## Skip's Original Introduction + +The Forth language is at an important cross-road with regard to +its use as a general scientific programming language. The new +FORTRAN-90 is just becoming available and long time FORTRAN +programmers are finding it different enough that many are +wondering if they might as well learn a new language instead of +sticking with FORTRAN. If the Forth community plays its hand +right, that alternate language could be Forth. To do so Forth +needs to overcome the standard complaints of the FORTRAN +community: + +(1) Its not standardized, so how can I port my software ? +(2) I have lots of pre-existing FORTRAN code that works perfectly + well and I am not in any hurry to re-write it. Can Forth + interface with my FORTRAN code ? +(3) There are no scientific libraries in Forth. + +The recently adopted ANS Forth handily addresses #1, and in fact +it is the adoption of the standard that makes issues #2 and #3 +worth addressing. + +With regard to #2, I think the adoption of the standard will help +here, since the interface to other software is the kind of feature +that will distinguish one vendors ANS Forth from anothers. While +the standard does not address such interfaces, I don't think that +there will be too much divergence on how this is done. The Unix +world has no such standard, and I have only encountered 2 different +C-FORTRAN conventions in over 15 years of using Unix. + +So #1 is now solved, and the vendors will (I hope!) address #2. The +third point can be addressed by the Forth community itself. Several +potential scientific users of Forth discussed these issues at the +1994 Rochester Forth Conference. It was decided that we should +undertake the project of writing a scientific library in ANS Forth. + +The plan is to write a set of Forth words to implement such libraries +as the ACM libraries, BLAS, LINPACK, etc. The libraries will be +publicly available in source form (in some sort of "public" release: +public domain, copyleft, copyrighted but freely distributable, etc). +.... + + Everett (Skip) Carter Phone: 831-641-0645 FAX: 831-641-0647 + Taygeta Scientific Inc. INTERNET: skip@taygeta.com + 1340 Munras Ave., Suite 314 UUCP: ...!uunet!taygeta!skip + Monterey, CA. 93940 WWW: http://www.taygeta.com/skip.html diff --git a/example/fsl/cube-rt.forth b/example/fsl/cube-rt.forth new file mode 100644 index 0000000..fc70d61 --- /dev/null +++ b/example/fsl/cube-rt.forth @@ -0,0 +1,52 @@ +\ Cube root of real number by Newton's method +\ ANS compatible version V1.2 10/6/1994 + +\ Forth Scientific Library Algorithm #5 + +\ This code conforms with ANS requiring: +\ The FLOAT and FLOAT EXT word sets +\ Non STANDARD words +\ : FTUCK ( F: x y -- y x y) FSWAP FOVER ; +\ : F2* ( F: x -- 2x ) FDUP F+ ; +\ : F**2 FDUP F* ; + + +\ (c) Copyright 1994 Julian V. Noble. Permission is granted +\ by the author to use this software for any application provided +\ the copyright notice is preserved. + + +3 S>D D>F FCONSTANT F=3 + +: X' ( F: N x -- x') + FTUCK F**2 F/ FSWAP F2* F+ F=3 F/ ; + +\ The magic number 1E-8 needs no change, even when extended (80-bit) precision +\ is needed. +: CONVERGED? ( F: x' x x' --) ( -- f) + F- FOVER F/ FABS 1.0E-8 F< ; + +: FCBRT ( F: N -- N^1/3) + FDUP F0< FABS ( F: -- |N|) ( -- f) + FDUP FSQRT ( F: -- N x0 ) + BEGIN FOVER FOVER X' FTUCK CONVERGED? UNTIL + X' IF FNEGATE THEN ; + +~~~ +:x' (f:nx-X) + f:tuck f:dup f:* f:/ f:swap f:dup f:+ f:+ .3 f:/ ; + +:converged? (f:XxX-) (-f) + f:- f:over f:/ f:abs .1.0e-8 f:lt? ; + +:fcbrt (f:n-[n^1/3]) + f:dup .0 f:lt? f:abs + f:dup f:sqrt + [ f:over f:over x' f:tuck converged? ] until + x' n:-zero? [ f:negate ] if ; +~~~ + +~~~ +.9 fcbrt f:put nl +~~~ + diff --git a/example/fsl/elip.forth b/example/fsl/elip.forth new file mode 100644 index 0000000..c2add24 --- /dev/null +++ b/example/fsl/elip.forth @@ -0,0 +1,114 @@ +\ elip Complete Elliptic Integral ACM Algorithm #149 + +\ Forth Scientific Library Algorithm #2 + +\ Evaluates the Complete Elliptic Integral, +\ Elip[a, b] = int_0^{\pi/2} 1/Sqrt{a^2 cos^2(t) + b^2 sin^2(t)} dt + +\ This function can be used to evaluate the complete elliptic integral +\ of the first kind, by using the relation K[m] = a Elip[a,b], m = 1 - b^2/a^2 + +\ This code conforms with ANS requiring: +\ 1. The Floating-Point word set +\ 2. The immediate word '%' which takes the next token +\ and converts it to a floating-point literal +\ 3. The FCONSTANT PI (3.1415926536...) +\ +\ Both a recursive form and an iterative form are given, but because of the +\ large stack consumption the recursive form is probably not of much +\ practical use. + +\ Caution: this code can potentially go into an endless loop +\ for certain values of the parameters. + +\ Collected Algorithms from ACM, Volume 1 Algorithms 1-220, +\ 1980; Association for Computing Machinery Inc., New York, +\ ISBN 0-89791-017-6 + +\ (c) Copyright 1994 Everett F. Carter. Permission is granted by the +\ author to use this software for any application provided the +\ copyright notice is preserved. + +CR .( ELIP V1.2 21 September 1994 EFC ) + + +: elip1 ( --, f: a b -- elip[a,b] ) \ recursive form + + FOVER FOVER FOVER F- FABS + FSWAP % 1.0e-8 F* + F< IF + FDROP + pi % 2.0 F/ FSWAP F/ + ELSE + FOVER FOVER F+ % 2.0 F/ + FROT FROT F* FSQRT + RECURSE + THEN +; + +: elip2 ( --, f: a b -- elip[a,b] ) \ nonrecursive version + + BEGIN + FOVER FOVER F+ % 2.0 F/ + FROT FROT F* FSQRT + + FOVER FOVER FOVER F- FABS + FSWAP % 1.0e-8 F* + F< UNTIL + + FDROP + + pi % 2.0 F/ FSWAP F/ + +; + + +\ test driver, calculates the complete elliptic integral of the first +\ kind (K(m)) using the relation: K[m] = a Elip[a,b], m = 1 - b^2/a^2 +\ compare with Abramowitz & Stegun, Handbook of Mathematical Functions, +\ Table 17.1 + +: elip_test ( -- ) + CR + ." m K(m) exact a Elip1[a,b] a Elip2[a,b] " CR + + ." 0.0 1.57079633 " + % 1000.0 % 1000.0 elip1 % 1000.0 F* F. ." " + % 1000.0 % 1000.0 elip2 % 1000.0 F* F. CR + + ." 0.44 1.80632756 " + % 400.0 % 299.33259 elip1 % 400.0 F* F. ." " + % 400.0 % 299.33259 elip2 % 400.0 F* F. CR + + ." 0.75 2.15651565 " + % 1000.0 % 500.0 elip1 % 1000.0 F* F. ." " + % 1000.0 % 500.0 elip2 % 1000.0 F* F. CR + + ." 0.96 3.01611249 " + % 500.0 % 100.0 elip1 % 500.0 F* F. ." " + % 500.0 % 100.0 elip2 % 500.0 F* F. CR + + +; + +# The Code + +~~~ +:elip1 (-,f:ab-elip[a,b]) + f:over f:over f:over f:- f:abs + f:swap .1.0e-8 f:* + f:lt? [ f:drop + f:PI .2.0 f:/ f:swap f:/ + ] + [ f:over f:over f:+ .2.0 f:/ + f:rot f:rot f:* f:sqrt + elip1 + ] choose +; +~~~ + +# Test + +``` +.1000 .1000 elip1 .1000 f:* f:put nl +``` diff --git a/example/fsl/expint.forth b/example/fsl/expint.forth new file mode 100644 index 0000000..b31a607 --- /dev/null +++ b/example/fsl/expint.forth @@ -0,0 +1,146 @@ +\ expint Real Exponential Integral ACM Algorithm #20 + +\ Forth Scientific Library Algorithm #1 + +\ Evaluates the Real Exponential Integral, +\ E1(x) = - Ei(-x) = int_x^\infty exp^{-u}/u du for x > 0 +\ using a rational approximation + +\ This code conforms with ANS requiring: +\ 1. The Floating-Point word set +\ 2. The immediate word '%' which takes the next token +\ and converts it to a floating-point literal +\ + +\ Collected Algorithms from ACM, Volume 1 Algorithms 1-220, +\ 1980; Association for Computing Machinery Inc., New York, +\ ISBN 0-89791-017-6 + +\ (c) Copyright 1994 Everett F. Carter. Permission is granted by the +\ author to use this software for any application provided the +\ copyright notice is preserved. + + +CR .( EXPINT V1.1 21 September 1994 EFC ) + + +: expint ( --, f: x -- expint[x] ) + FDUP + % 1.0 F< IF + FDUP % 0.00107857 F* % 0.00976004 F- + FOVER F* + % 0.05519968 F+ + FOVER F* + % 0.24991055 F- + FOVER F* + % 0.99999193 F+ + FOVER F* + % 0.57721566 F- + FSWAP FLN F- + ELSE + FDUP % 8.5733287401 F+ + FOVER F* + % 18.059016973 F+ + FOVER F* + % 8.6347608925 F+ + FOVER F* + % 0.2677737343 F+ + + FOVER + FDUP % 9.5733223454 F+ + FOVER F* + % 25.6329561486 F+ + FOVER F* + % 21.0996530827 F+ + FOVER F* + % 3.9584969228 F+ + + FSWAP FDROP + F/ + FOVER F/ + FSWAP % -1.0 F* FEXP + F* + + THEN +; + + +\ test code generates a small table of selected E1 values. +\ most comparison values are from Abramowitz & Stegun, +\ Handbook of Mathematical Functions, Table 5.1 + +: expint_test ( -- ) + + CR + ." x E1(x) exact ExpInt[x] " CR + + + ." 0.5 0.5597736 " + % 0.5 expint F. CR + + ." 1.0 0.2193839 " + % 1.0 expint F. CR + + ." 2.0 0.0489005 " + % 2.0 expint F. CR + + ." 5.0 0.001148296 " + % 5.0 expint F. CR + + ." 10.0 0.4156969e-5 " + % 10.0 expint F. CR + +; + +# The Code + +~~~ +:expint (-,f:x-expint[x]) + f:dup + .1.0 f:lt? [ + f:dup .0.00107857 f:* .0.00976004 f:- + f:over f:* + .0.05519968 f:+ + f:over f:* + .0.24991055 f:- + f:over f:* + .0.99999193 f:+ + f:over f:* + .0.57721566 f:- + f:swap f:E f:log f:- + ] [ + f:dup .8.5733287401 f:+ + f:over f:* + .18.059016973 f:+ + f:over f:* + .8.6347608925 f:+ + f:over f:* + .0.2677737343 f:+ + + f:over + f:dup .9.5733223454 f:+ + f:over f:* + .25.6329561486 f:+ + f:over f:* + .21.0996530827 f:+ + f:over f:* + .3.9584969228 f:+ + + f:swap f:drop + f:/ + f:over f:/ + f:swap .-1.0 f:* f:E f:swap f:power + f:* + ] choose +; +~~~ + +# Tests + +``` +.0.5 expint f:put nl +.1.0 expint f:put nl +.2.0 expint f:put nl +.5.0 expint f:put nl +.10.0 expint f:put nl +``` diff --git a/example/fsl/logistic.forth b/example/fsl/logistic.forth new file mode 100644 index 0000000..b39ecd0 --- /dev/null +++ b/example/fsl/logistic.forth @@ -0,0 +1,63 @@ +\ logistic The Logistic function and its first derivative +\ logistic = Exp( c + a x ) / (1 + Exp( c + a x ) ) +\ d_logistic = a Exp( c + a x ) / (1 + Exp( c + a x ) )^2 + +\ Forth Scientific Library Algorithm #4 + +\ This code conforms with ANS requiring: +\ 1. The Floating-Point word set +\ + +\ (c) Copyright 1994 Everett F. Carter. Permission is granted by the +\ author to use this software for any application provided this +\ copyright notice is preserved. + + +cr .( Logistic V1.2 17 October 1994 EFC ) + + +: logistic ( --, f: x a c -- z ) + FROT FROT + F* F+ + FEXP + FDUP 1.0e0 F+ + F/ +; + +: d_logistic ( -- , f: x a c -- z ) + FSWAP FROT + FOVER F* FROT F+ + FEXP + + FDUP 1.0e0 F+ FDUP F* + F/ F* +; + +\ Examples % 1.0 % 1.0 % 0.0 logistic f. 0.731059 +\ % 3.2 % 1.5 % 0.2 logistic f. 0.993307 +\ % 3.2 % 1.5 % 0.2 d_logistic f. 0.00997209 + +# The Code + +~~~ +:logistic (-,f:xac-z) + f:rot f:rot + f:* f:+ + f:E f:swap f:power + f:dup .1.0e0 f:+ + f:/ ; + +:d_logistic (-,f:xac-z) + f:swap f:rot + f:over f:* f:rot f:+ + f:E f:swap f:power + f:dup .1.0e0 f:+ f:dup f:* + f:/ f:* ; +~~~ + +# Tests + +``` +.1.0 .1.0 .0.0 logistic f:put nl +.3.2 .1.5 .0.2 logistic f:put nl +```