\ Regular spherical Bessel functions j_n(x), n=0-9 \ \ Uses Miller's method of downward recursion, as described \ in Abramowitz & Stegun, "Handbook of Mathematical Functions" \ 10.5 ff. The recursion is \ \ j(n-1) = (2n+1) j(n) / x - j(n+1) \ \ The downward recursion is started with j40 = 0, j39 = 1 . The \ resulting functions are normalized using \ \ Sum (n=0 to inf) { (2n+1) * jn(x)^2 } = 1 . \ \ Usage: % 3 SPHBES leaves jn(3), n=0-9, \ in the double-length (64-bit) array JBES{ \ \ Programmed by J.V. Noble 6/15/99 \ \ Assembly language version suitable for Win32Forth marker -jbes include arrays.f 40 long 1 dfloats 1array jbes{ FVARIABLE x HEX code ITERATE ( f: x -- ) \ initialization finit \ clear fpu stack mov ecx, FSP [edi] sub ecx, # B/FLOAT js L$2 \ error handler fld FSIZE FSTACK [ecx] [edi] ( 87: x) mov FSP [edi], ecx push ebx push # 4F \ 79 on data stack fldz ( 87: 0 x) fild dword 0 [esp] ( 87: 79 0 x) fldz fld1 ( 87: 1 0 79 0 x) pop ebx \ ebx = 79 ( 87: jn jn+1 2n+1 sum x) \ end of initialization L$1: dec ebx \ loop begins here fst double jbes{ 0 } [ebx*4] [edi] \ fwait \ may be needed fxch st(1) ( 87: jn+1 jn k=2n+1 sum x) fld st(1) ( 87: jn jn+1 jn k sum x) fmul st(0), st(3) ( 87: k*jn jn+1 jn k sum x) fld st(0) ( 87: k*jn k*jn jn+1 jn k sum x) fmul st(0), st(3) ( 87: k*jn^2 k*jn jn+1 jn k sum x) faddp st(5), st(0) ( 87: k*jn jn+1 jn k sum' x) fdiv st(0), st(5) ( 87: k*jn/x jn+1 jn k sum' x) fsubpr st(1), st(0) \ this is a sp. error in 486asm.f \ ^^^^^^ \ -- should be fsubrp fld1 fsub st(3), st(0) fsubp st(3), st(0) ( 87: jn-1 jn 2n-1 sum' x) dec ebx jns L$1 \ loop ends here ( 87: j0 j1 -1 sum x) fstp st(0) ( 87: j1 1 sum x) fstp st(0) ( 87: 1 sum x) fstp st(0) ( 87: sum x) mov ecx, FSP [edi] fstp FSIZE FSTACK [ecx] [edi] fwait add ecx, # B/FLOAT mov FSP [edi], ecx fstp st(0) ( 87: x -- ) pop ebx ( n' --) jmp L$3 L$2: mov esi, # ' FSTKUFLO >body add esi, edi L$3: next, end-code DECIMAL : DO_X=0 FDROP F1.0 JBES{ 0 } DF! 10 1 DO F0.0 JBES{ I } DF! LOOP ; : NORMALIZE ( f: sum --) FSQRT F1.0 FSWAP F/ 39 0 DO FDUP JBES{ I } DUP F@ F* F! LOOP FDROP ; : SPHBES ( F: x --) FDUP F0= IF DO_X=0 EXIT THEN ITERATE NORMALIZE ; FALSE [IF] To calculate j0-j39 say, e.g., 3.7e0 sphbes To access/display a value say, e.g., jbes{ 3 } F@ FS. [THEN]