\ 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]