Plan 9 from Bell Labs’s /usr/web/sources/contrib/fgb/root/sys/src/cmd/4th/lib/pcylfun.4th

Copyright © 2021 Plan 9 Foundation.
Distributed under the MIT License.
Download the Plan 9 distribution.


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



Bell Labs OSI certified Powered by Plan 9

(Return to Plan 9 Home Page)

Copyright © 2021 Plan 9 Foundation. All Rights Reserved.
Comments to [email protected].