crc/retroforth
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

examples: add in a few files I ported from the FSL

Browse source 
FossilOrigin-Name: b4efc49cb970beb1aac0bc90f8550948d67a7367690d9b6b315334de6c0e69f0
This commit is contained in:
crc 2019-03-14 12:03:34 +00:00
parent 81a8db3dc3
commit 18ae06fa35
5 changed files with 446 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

71
example/fsl/README.txt Normal file
View file

@ -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

52
example/fsl/cube-rt.forth Normal file
View file

@ -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
~~~

114
example/fsl/elip.forth Normal file
View file

@ -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
```

146
example/fsl/expint.forth Normal file
View file

@ -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
```

63
example/fsl/logistic.forth Normal file
View file

@ -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
```