\ pcylfun Parabolic Cylinder Functions U and V,
\ plus related confluent hypergeometric functions
\ Forth Scientific Library Algorithm #20
\ Evaluates the Parabolic Cylinder functions,
\ Upcf U(nu,x)
\ and
\ Vpcf V(nu,x)
\ In addition the following related functions are provided,
\ U() U(a,b,x) Hypergeometric function U for real args
\ M() M(a,b,x) Hypergeometric function M for real args
\ Wwhit W(k,mu,z) Whittaker function W for real args
\ Mwhit M(k,mu,z) Whittaker function M for real args
\ 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. Uses the word 'GAMMA' to evaluate the gamma function
\ 4. The FCONSTANT PI (3.1415926536...)
\ 5. The compilation of the test code is controlled by TRUE/FALSE
\ and the conditional compilation words in the Programming-Tools wordset.
\ 6. Under 4tH this function is neither fast nor precise! (about 5 digits).
\ There is a bit of stack gymnastics in this code particularly in U() and M()
\ but that seems to be in the nature of the algorithm.
\ Baker, L., 1992; C Mathematical Function Handbook, McGraw-Hill,
\ New York, 757 pages, ISBN 0-07-911158-0
\ (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.
\ PCYLFUN V1.1 3 October 1994 EFC
\ include lib/ansfloat.4th
\ include lib/fsl-util.4th
\ Private:
[UNDEFINED] Upcf [IF]
[UNDEFINED] F** [IF] [NEEDS lib/falog.4th] [THEN]
[UNDEFINED] GAMMA [IF] [NEEDS lib/gamma.4th] [THEN]
FLOAT array SUM
FLOAT array TERM
FLOAT array OLD-TERM
\ scratch space for stack conservation
FLOAT array XX-TMP
FLOAT array A-TMP
FLOAT array B-TMP
FLOAT array C-TMP
FLOAT array D-TMP
FLOAT array U-TMP
FLOAT array AV-TMP
FLOAT array BV-TMP
FLOAT array CV-TMP
FLOAT array XV-TMP
FLOAT array NV-TMP
FLOAT array XU-TMP
FLOAT array NU-TMP
FLOAT array AU-TMP
FLOAT array BU-TMP
: FRAC ( -- , f: x -- fractional_part_of_x )
FDUP F>S S>F F-
;
: FAC ( -- , f: x -- z )
1 S>F F+ GAMMA
;
: ?big-x ( -- t/f , f: nu x -- nu x )
FOVER FOVER
4 S>F F> FABS 4 S>F F* FOVER FSWAP F> AND
;
: asymptotic-u ( -- , f: nu x -- U[nu,x] )
S" 0.5" S>FLOAT FOVER FDUP F* F/
A-TMP F!
FDUP FDUP F* S" -0.25" S>FLOAT F*
FEXP D-TMP F!
FOVER S" 0.5" S>FLOAT F+ F** D-TMP F@ FSWAP F/
D-TMP F!
FDUP FDUP
S" 3.5" S>FLOAT F+ FSWAP S" 2.5" S>FLOAT F+ F* A-TMP F@ F* S" 0.5" S>FLOAT F* 1 S>F F-
A-TMP F@ F*
FOVER S" 0.5" S>FLOAT F+ F*
FSWAP S" 1.5" S>FLOAT F+ F*
1 S>F F+
D-TMP F@ F*
;
: simple-u ( -- , f: b a -- z )
FSWAP FDROP XX-TMP F@ FSWAP F**
SUM F@ FSWAP F/
;
: asymptotic-v ( -- , f: nu x -- V[nu,x] )
S" 0.5" S>FLOAT FOVER FDUP F* F/
A-TMP F!
FDUP FDUP F* S" 0.25" S>FLOAT F*
FEXP D-TMP F!
FOVER S" 0.5" S>FLOAT F- F** D-TMP F@ F*
2 S>F PI F/ FSQRT F*
D-TMP F!
FDUP FDUP
S" 3.5" S>FLOAT F- FSWAP S" 2.5" S>FLOAT F- F* A-TMP F@ F* S" 0.5" S>FLOAT F* 1 S>F F+
A-TMP F@ F*
FOVER S" 0.5" S>FLOAT F- F*
FSWAP S" 1.5" S>FLOAT F- F*
1 S>F F+
D-TMP F@ F*
;
: ?bad-M-parms ( -- t/f , f: a b -- a b )
FDUP F0<
FOVER FOVER FSWAP F< AND
FOVER FRAC F0= AND
FDUP FRAC F0= AND
;
\ Public:
: M() ( -- , f: a b x --- u )
XX-TMP F!
?bad-M-parms abort" M() bad parameters "
XX-TMP F@ F0= IF
FDROP FDROP 1 S>F
ELSE
1 S>F TERM F! 1 S>F SUM F!
XX-TMP F@ 10 S>F F>
IF
S" 1.0e30" S>FLOAT OLD-TERM F!
40 0 DO
FOVER FOVER FOVER F- I S>F F+
FSWAP I 1+ S>F FSWAP F- F*
I 1+ S>F XX-TMP F@ F* F/
TERM F@ F* FDUP TERM F!
FABS OLD-TERM F@ F> IF LEAVE THEN
TERM F@ FDUP FABS OLD-TERM F!
SUM F@ F+ SUM F!
OLD-TERM F@ FDUP S" 1.0e-6" S>FLOAT F<
SUM F@ FABS S" 1.0e-6" S>FLOAT F* F<
OR IF LEAVE THEN
LOOP
FOVER FOVER F- XX-TMP F@ FSWAP F**
XX-TMP F@ FEXP F* SUM F@ F*
FSWAP GAMMA F*
FSWAP GAMMA F/
ELSE
40 0 DO
FOVER I S>F F+ XX-TMP F@ F*
FOVER I S>F F+
I 1+ S>F F* F/
TERM F@ F* FDUP TERM F!
FABS
SUM F@ FABS S" 1.0e-6" S>FLOAT F*
F< IF LEAVE THEN
SUM F@ TERM F@ F+ SUM F!
LOOP
FDROP FDROP SUM F@
THEN
THEN
;
\ Private:
: u-small-x ( -- , f: a b -- u )
\ wont work if b is an integer, so tweak it slightly if it is
FDUP FRAC F0= IF S" 1.0e-6" S>FLOAT F+ THEN
BU-TMP F! AU-TMP F!
XX-TMP F@ 0 S>F F> IF \ 0 > x > 5
PI BU-TMP F@ PI F* FSIN F/
AU-TMP F@ BU-TMP F@ F- 1 S>F F+ GAMMA
BU-TMP F@ GAMMA F* U-TMP F!
AU-TMP F@ BU-TMP F@ XX-TMP F@ M() U-TMP F@ F/
U-TMP F!
BU-TMP F@ FNEGATE FDUP BU-TMP F! \ b is now -b
1 S>F F+ FAC
AU-TMP F@ GAMMA F*
XX-TMP F@ BU-TMP F@ 1 S>F F+ F**
FSWAP F/
AU-TMP F@ BU-TMP F@ F+ 1 S>F F+
2 S>F BU-TMP F@ F+
XX-TMP F@ M()
F* FNEGATE
U-TMP F@ F+ F*
ELSE \ -5 < x < 0
PI BU-TMP F@ PI F* FSIN F/
XX-TMP F@ FEXP F*
BU-TMP F@ AU-TMP F@ F-
BU-TMP F@
XX-TMP F@ FNEGATE M()
U-TMP F!
\ NOTE: side effect recovery !!!!
\ M() stores last parameter in XX-TMP
\ which is now the original XX-TMP
\ with it sign changed, so we have
\ to fix it here,
XX-TMP F@ FNEGATE XX-TMP F!
AU-TMP F@ BU-TMP F@ F- 1 S>F F+ GAMMA
BU-TMP F@ GAMMA F* U-TMP F@ FSWAP F/
U-TMP F!
BU-TMP F@ FNEGATE FDUP BU-TMP F! \ b is now -b
2 S>F F+ GAMMA AU-TMP F@ GAMMA F*
BU-TMP F@ 1 S>F F+ PI F* FCOS FSWAP F/
XX-TMP F@ FNEGATE
BU-TMP F@ 1 S>F F+ F** F*
1 S>F AU-TMP F@ F-
BU-TMP F@ 2 S>F F+
XX-TMP F@ M()
F* FNEGATE
U-TMP F@ F+ F*
THEN
;
\ Public:
: U() ( -- , f: a b x --- u )
FDUP XX-TMP F!
FABS 5 S>F F< IF
u-small-x
ELSE
1 S>F TERM F! 1 S>F SUM F!
FSWAP
40 0 DO
FOVER FOVER FSWAP F- I 1+ S>F F+
FOVER I S>F F+ F*
I 1+ S>F XX-TMP F@ F* F/ FNEGATE
FDUP FABS 1 S>F F> IF
TERM F@ SUM F@ F/ FABS S" 1.0e-3" S>FLOAT F<
IF FDROP simple-u
ELSE
FDROP FSWAP u-small-x
THEN
LEAVE
ELSE
TERM F@ F* FDUP TERM F!
SUM F@ F+ SUM F!
TERM F@ SUM F@ F/ FABS S" 1.0e-6" S>FLOAT F<
IF simple-u
LEAVE
THEN
THEN
LOOP
THEN
;
: Upcf ( -- , f: nu x -- U[nu,x] )
?big-x IF
asymptotic-u
ELSE
FSWAP S" 0.5" S>FLOAT F* FSWAP \ nu is now 0.5*nu
FDUP FDUP F* S" -0.25" S>FLOAT F* FEXP A-TMP F!
FOVER S" 0.25" S>FLOAT F+ 2 S>F FSWAP F**
A-TMP F@ FSWAP F/ A-TMP F!
FDUP 0 S>F F> IF
FDUP F* S" 0.5" S>FLOAT F*
FSWAP S" 0.25" S>FLOAT F+ FSWAP
S" 0.5" S>FLOAT FSWAP
U()
ELSE \ x <= 0
PI FSQRT A-TMP F@ F* A-TMP F!
\ saving params in TMPs to conserve stack space
XU-TMP F! FDUP NU-TMP F!
S" 0.75" S>FLOAT F+ GAMMA B-TMP F!
NU-TMP F@ S" 0.25" S>FLOAT F+ GAMMA C-TMP F!
NU-TMP F@ S" 0.25" S>FLOAT F+
S" 0.5" S>FLOAT
XU-TMP F@ FDUP F* S" 0.5" S>FLOAT F*
M() B-TMP F@ F/ B-TMP F!
NU-TMP F@ S" 0.75" S>FLOAT F+
S" 1.5" S>FLOAT
XU-TMP F@ FDUP F* S" 0.5" S>FLOAT F*
M() C-TMP F@ F/ 2 S>F FSQRT F* XU-TMP F@ F* FNEGATE
B-TMP F@ F+
THEN
A-TMP F@ F*
THEN
;
: Vpcf ( -- , f: nu x -- V[nu,x] )
?big-x IF
asymptotic-v
ELSE
\ saving params in TMPs to conserve stack space
XV-TMP F! FDUP NV-TMP F!
S" 0.5" S>FLOAT F+ GAMMA PI F/ AV-TMP F!
NV-TMP F@ PI F* FSIN BV-TMP F!
NV-TMP F@ XV-TMP F@
Upcf CV-TMP F!
NV-TMP F@ XV-TMP F@ FNEGATE
Upcf D-TMP F!
BV-TMP F@ CV-TMP F@ F* D-TMP F@ F+ AV-TMP F@ F*
THEN
;
: Mwhit ( -- , f: k mu z -- M[k,mu,z] )
FOVER FOVER FSWAP S" 0.5" S>FLOAT F+ FSWAP F**
FOVER S" -0.5" S>FLOAT F* FEXP F*
D-TMP F!
FROT FNEGATE S" 0.5" S>FLOAT F+
FROT FSWAP FOVER F+
FSWAP 2 S>F F* 1 S>F F+
FROT M()
D-TMP F@ F*
;
: Wwhit ( -- , f: k mu z -- W[k,mu,z] )
FOVER FOVER FSWAP S" 0.5" S>FLOAT F+ FSWAP F**
FOVER S" -0.5" S>FLOAT F* FEXP F*
D-TMP F!
FROT FNEGATE S" 0.5" S>FLOAT F+
FROT FSWAP FOVER F+
FSWAP 2 S>F F* 1 S>F F+
FROT U()
D-TMP F@ F*
;
[DEFINED] 4TH# [IF]
hide SUM
hide TERM
hide OLD-TERM
hide XX-TMP
hide A-TMP
hide B-TMP
hide C-TMP
hide D-TMP
hide U-TMP
hide AV-TMP
hide BV-TMP
hide CV-TMP
hide XV-TMP
hide NV-TMP
hide XU-TMP
hide NU-TMP
hide AU-TMP
hide BU-TMP
hide FRAC
hide FAC
hide ?big-x
hide asymptotic-u
hide asymptotic-v
hide simple-u
hide ?bad-M-parms
hide u-small-x
[THEN]
[THEN]
false [IF] \ Test Code =============================================
\ compares against test values given in Baker, 1992
: pcf-test ( -- )
CR
." nu x Upcf-actual Vpcf-actual Upcf Vpcf " CR
." -2.0 0.0 -0.6081402 -0.4574753 "
-2 S>F 0 S>F Upcf F. -2 S>F 0 S>F Vpcf F. CR
." 2.0 0.0 0.8108537 0.3431063 "
2 S>F 0 S>F Upcf F. 2 S>F 0 S>F Vpcf F. CR
." -2.0 1.0 0.5156671 -0.5417693 "
-2 S>F 1 S>F Upcf F. -2 S>F 1 S>F Vpcf F. CR
." 2.0 1.0 0.1832067 1.439015 "
2 S>F 1 S>F Upcf F. 2 S>F 1 S>F Vpcf F. CR
;
fclear 10 set-precision
pcf-test
[THEN]
|
