From 7ad153c190d20c4fbbb3d0a0e2f63027245d3c4a Mon Sep 17 00:00:00 2001 From: Tim Daly Date: Sun, 27 Dec 2015 17:53:38 -0500 Subject: [PATCH] src/interp/sfsfun.lisp removed, merged with bookvol5 Goal: literate interpreter Move the lisp support code for DoubleFloatSpecialFunctions into the general algebra support section of the interpreter. --- books/bookvol5.pamphlet | 67 ++++- changelog | 4 + patch | 2 +- src/axiom-website/patches.html | 2 + src/interp/Makefile.pamphlet | 5 +- src/interp/sfsfun.lisp.pamphlet | 771 --------------------------------------- 6 files changed, 77 insertions(+), 774 deletions(-) delete mode 100644 src/interp/sfsfun.lisp.pamphlet diff --git a/books/bookvol5.pamphlet b/books/bookvol5.pamphlet index d21434f..8de89ee 100644 --- a/books/bookvol5.pamphlet +++ b/books/bookvol5.pamphlet @@ -48306,6 +48306,70 @@ AXIOM recurrence for $u_{k}$ \end{chunk} +\defun{chebf01}{chebf01} +Where: + +$c$ parameter to 0F1, possibly complex\\ +$z$ argument to 0F1\\ +$w$ scale factor so that $0 < \frac{z}{w} < 1$\\ +$n$, $n+2$ coefficients will be produced stored in an array +indexed from $0$ to $n+1$. + +Program transcribed from Fortran,, p. 80 \cite{Luke77} + +\begin{chunk}{defun chebf01} +(defun chebf01 (c z) + (let (cc temp b p sum rho divfac c1 x1 arr + ncount z1 a1 a2 a3 n2 n1 start four w n) + (setq n 75) ; ad hoc decision + (setq w (* 2.0 z)) + ; arr will be used to store the Cheb. series coefficients + (setq four 4.0) + (setq start (expt 10.0 -200)) + (setq n1 (+ n 1)) + (setq n2 (+ n 2)) + (setq a3 0.0) + (setq a2 0.0) + (setq a1 start) ; arbitrary starting value + (setq z1 (/ four w)) + (setq ncount n1) + (setq arr (make-array n2)) + (setf (aref arr ncount) start) ; start off + (setq x1 n2) + (setq c1 (- 1.0 c)) + (loop for ncount from n downto 0 by 1 do + (setq divfac (/ 1.0 x1)) + (setq x1 (- x1 1.0)) + (setf (aref arr ncount) + (* x1 (- (+ (* (+ divfac (* z1 (- x1 c1))) a1) + (* (+ (/ 1.0 x1) (* z1 (+ (+ x1 c1) 1.0))) a2)) + (* divfac a3)))) + (setq a3 a2) + (setq a2 a1) + (setq a1 (aref arr ncount))) + (setf (aref arr 0) (/ (aref arr 0) 2.0)) + ; compute scale factor + (setq rho (aref arr 0)) + (setq sum rho) + (setq p 1.0) + (loop for i from 1 to n1 do + (setq rho (- rho (* p (aref arr i)))) + (setq sum (+ sum (aref arr i))) + (setq p (- p))) + (loop for m from 0 to n1 do + (setf (aref arr m) (/ (aref arr m) rho))) + (setq sum (/ sum rho)) + ; Now evaluate array at argument + (setq b 0.0) + (setq temp 0.0) + (loop for i from (+ n 1) downto 0 by 1 do + (setq cc b) + (setq b temp) + (setq temp (+ (- cc) (aref arr i)))) + temp)) + +\end{chunk} + \chapter{Common Lisp Algebra Support} These functions are called directly from the algebra source code. They fall into two basic categories, one are the functions that are @@ -48728,7 +48792,7 @@ up to overflow. See Hart and Cheney. \calls{chyper0f1}{s-to-c} \begin{chunk}{defun chyper0f1} (defun chyper0f1 (a z) - (c-to-s (|chebf01| (s-to-c a) (s-to-c z)) )) + (c-to-s (chebf01 (s-to-c a) (s-to-c z)) )) \end{chunk} @@ -63244,6 +63308,7 @@ Note that ACL2 does not support FLOATP. \getchunk{defun changeHistListLen} \getchunk{defun changeToNamedInterpreterFrame} \getchunk{defun charDigitVal} +\getchunk{defun chebf01} \getchunk{defun checkCondition} \getchunk{defun checkFilter} \getchunk{defun chkAllNonNegativeInteger} diff --git a/changelog b/changelog index f44b21c..8ce58c3 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,7 @@ +20151227 tpd src/axiom-website/patches.html 20151227.01.tpd.patch +20151227 tpd books/bookvol5 move more functions from sfsfun +20151227 tpd src/interp/sfsfun.lisp removed, merged with bookvol5 +20151227 tpd src/interp/Makefile remove sfsfun.lisp 20151226 tpd src/axiom-website/patches.html 20151226.01.tpd.patch 20151226 tpd books/bookvol5 move more functions from sfsfun 20151226 tpd src/interp/sfsfun.lisp move more functions to bookvol5 diff --git a/patch b/patch index a1517df..95bff10 100644 --- a/patch +++ b/patch @@ -1,4 +1,4 @@ -books/bookvol5 move more functions from sfsfun +src/interp/sfsfun.lisp removed, merged with bookvol5 Goal: literate interpreter diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index f554be8..e663fc0 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -5172,6 +5172,8 @@ books/bookvol5 move more functions from sfsfun
books/bookvol5 move more functions from sfsfun
20151226.01.tpd.patch books/bookvol5 move more functions from sfsfun
+20151227.01.tpd.patch +src/interp/sfsfun.lisp removed, merged with bookvol5
diff --git a/src/interp/Makefile.pamphlet b/src/interp/Makefile.pamphlet index bf6e36d..5516201 100644 --- a/src/interp/Makefile.pamphlet +++ b/src/interp/Makefile.pamphlet @@ -177,7 +177,6 @@ OBJS= ${OUT}/vmlisp.${O} \ ${OUT}/record.${O} \ ${OUT}/rulesets.${O} \ ${OUT}/server.${O} \ - ${OUT}/sfsfun.${O} \ ${OUT}/simpbool.${O} ${OUT}/slam.${O} \ ${OUT}/sockio.${O} \ ${OUT}/template.${O} ${OUT}/termrw.${O} \ @@ -191,6 +190,10 @@ OBJS= ${OUT}/vmlisp.${O} \ \end{chunk} +\begin{verbatim} + sfsfun.lisp.pamphlet contains 0F1 code +\end{verbatim} + Before we save the {\bf SAVESYS} image we need to run some initialization code. These files perform initialization for various parts of the system. The {\bf patches.lisp} \cite{5} diff --git a/src/interp/sfsfun.lisp.pamphlet b/src/interp/sfsfun.lisp.pamphlet deleted file mode 100644 index a060edd..0000000 --- a/src/interp/sfsfun.lisp.pamphlet +++ /dev/null @@ -1,771 +0,0 @@ -\documentclass{article} -\usepackage{axiom} -\begin{document} -\title{\$SPAD/src/interp sfsfun.lisp} -\author{The Axiom Team} -\maketitle -\begin{abstract} -\end{abstract} -\eject -\tableofcontents -\eject -\begin{verbatim} -NOTEfrom TTT: at least BesselJAsymptOrder needs work - -1. This file contains the contents of BWC's original files: - floaterrors.boot - floatutils.boot - rgamma.boot - cgamma.boot - rpsi.boot - cpsi.boot - f01.boot - chebf01cmake.boot - chebevalsf.boot - besselIJ.boot - -2. All declarations have been commented out with "--@@" - since the boot translator is generating bad lisp code from them. - -4. BesselIJ was not compiling. I have modified the code from that - file to make it compile. It should be checked for correctness. - -SMW June 25, 1991 - -"Fixes" to BesselJ, B. Char June 14, 1992. Needs extensive testing and - further fixes to BesselI and BesselJ. -More fixes to BesselJ, T. Tsikas 24 Feb, 1995. - -\end{verbatim} -\begin{chunk}{*} - -(IN-PACKAGE "BOOT") - -;negintp(x) == -; if ZEROP IMAGPART(x) and x<0.0 and ZEROP fracpart(x) -; then -; true -; else -; false - -(defun |negintp| (x) - (cond - ((and (zerop (imagpart x)) (< x 0.0) - (zerop (fracpart x))) - t) - (t nil))) - -;-- Lisp PI is a long float and causes type errors, here we give -;-- enough digits to have double accuracy even after conversion -;-- to binary -;DEFCONSTANT(dfPi,3.14159265358979323846264338328) - -;--- Small float implementation of Gamma function. Valid for -;--- real arguments. Maximum error (relative) seems to be -;--- 2-4 ulps for real x 2l -; then -; t2 := 0 -; else -; t2 := -AREF(a,k+1)*(k+1) -; SETF(AREF(a,k), t1+t2) -; --- evaluate d^N/dX^N cot(z) via Horner-like rule -; v := COT(z) -; sq := v**2 -; s := AREF(a,n+1) -; for i in (n-1)..0 by -2 repeat -; s := s*sq + AREF(a,i) -; m := MOD(n,2) -; if m=0 -; then -; s := s*v -; if skipit=1 -; then -; if m=0 -; then -; return 0 -; else -; return AREF(a,0) -; else -; return s - -;cpsireflect(n,z) == -; m := MOD(n,2) -; sign := (-1)**m -; sign*dfPi**(n+1)*cotdiffeval(n,dfPi*z,0) + cPsi(n,1.0-z) - -(DEFUN |cpsireflect| (|n| |z|) - (PROG (|sign| |m|) - (RETURN - (PROGN - (SETQ |m| (MOD |n| 2)) - (SETQ |sign| (EXPT (- 1) |m|)) - (+ (* (* |sign| (EXPT Pi (+ |n| 1))) - (cotdiffeval |n| (* Pi |z|) 0)) - (cPsiImpl |n| (- 1.0 |z|))))))) - -;--- c parameter to 0F1, possibly complex -;--- z argument to 0F1 -;--- Depends on files floaterror, floatutils - -;--- Program transcribed from Fortran,, p. 80 Luke 1977 - -;chebf01 (c,z) == -;--- w scale factor so that 0 200.0 and ABS(z)>10000.0 -;--- then -;--- FloatError('"cheb0F1 not implemented for ~S < 1",[c,z]) -; w := 2.0*z -;--- arr will be used to store the Cheb. series coefficients -; four:= 4.0 -; start := EXPT(10.0, -200) -; n1 := n+1 -; n2 := n+2 -; a3 := 0.0 -; a2 := 0.0 -; a1 := start -- arbitrary starting value -; z1 := four/w -; ncount := n1 -; arr := MAKE_-ARRAY(n2) -; SETF(AREF(arr,ncount) , start) -- start off -; x1 := n2 -; c1 := 1.0 - c -; for ncount in n..0 by -1 repeat -; divfac := 1.0/x1 -; x1 := x1 -1.0 -; SETF(AREF(arr,ncount) ,_ -; x1*((divfac+z1*(x1-c1))*a1 +_ -; (1.0/x1 + z1*(x1+c1+1.0))*a2-divfac*a3)) -; a3 := a2 -; a2 := a1 -; a1 := AREF(arr,ncount) -; SETF(AREF(arr,0),AREF(arr,0)/2.0) -;-- compute scale factor -; rho := AREF(arr,0) -; sum := rho -; p := 1.0 -; for i in 1..n1 repeat -; rho := rho - p*AREF(arr,i) -; sum := sum+AREF(arr,i) -; p := -p -; for l in 0..n1 repeat -; SETF(AREF(arr,l), AREF(arr,l)/rho) -; sum := sum/rho -;--- Now evaluate array at argument -; b := 0.0 -; temp := 0.0 -; for i in (n+1)..0 by -1 repeat -; cc := b -; b := temp -; temp := -cc + AREF(arr,i) -; temp - -(DEFUN |chebf01| (|c| |z|) - (PROG (|cc| |temp| |b| |p| |sum| |rho| |divfac| |c1| |x1| |arr| - |ncount| |z1| |a1| |a2| |a3| |n2| |n1| |start| |four| |w| - |n|) - (RETURN - (PROGN - (SETQ |n| 75) - (SETQ |w| (* 2.0 |z|)) - (SETQ |four| 4.0) - (SETQ |start| (EXPT 10.0 (- 200))) - (SETQ |n1| (+ |n| 1)) - (SETQ |n2| (+ |n| 2)) - (SETQ |a3| 0.0) - (SETQ |a2| 0.0) - (SETQ |a1| |start|) - (SETQ |z1| (/ |four| |w|)) - (SETQ |ncount| |n1|) - (SETQ |arr| (MAKE-ARRAY |n2|)) - (SETF (AREF |arr| |ncount|) |start|) - (SETQ |x1| |n2|) - (SETQ |c1| (- 1.0 |c|)) - ((LAMBDA (|bfVar#12| |ncount|) - (LOOP - (COND - ((COND - ((MINUSP |bfVar#12|) (< |ncount| 0)) - (T (> |ncount| 0))) - (RETURN NIL)) - ('T - (PROGN - (SETQ |divfac| (/ 1.0 |x1|)) - (SETQ |x1| (- |x1| 1.0)) - (SETF (AREF |arr| |ncount|) - (* |x1| - (- (+ (* (+ |divfac| (* |z1| (- |x1| |c1|))) - |a1|) - (* (+ (/ 1.0 |x1|) - (* |z1| (+ (+ |x1| |c1|) 1.0))) - |a2|)) - (* |divfac| |a3|)))) - (SETQ |a3| |a2|) - (SETQ |a2| |a1|) - (SETQ |a1| (AREF |arr| |ncount|))))) - (SETQ |ncount| (+ |ncount| |bfVar#12|)))) - (- 1) |n|) - (SETF (AREF |arr| 0) (/ (AREF |arr| 0) 2.0)) - (SETQ |rho| (AREF |arr| 0)) - (SETQ |sum| |rho|) - (SETQ |p| 1.0) - ((LAMBDA (|i|) - (LOOP - (COND - ((> |i| |n1|) (RETURN NIL)) - ('T - (PROGN - (SETQ |rho| (- |rho| (* |p| (AREF |arr| |i|)))) - (SETQ |sum| (+ |sum| (AREF |arr| |i|))) - (SETQ |p| (- |p|))))) - (SETQ |i| (+ |i| 1)))) - 1) - ((LAMBDA (|l|) - (LOOP - (COND - ((> |l| |n1|) (RETURN NIL)) - ('T (SETF (AREF |arr| |l|) (/ (AREF |arr| |l|) |rho|)))) - (SETQ |l| (+ |l| 1)))) - 0) - (SETQ |sum| (/ |sum| |rho|)) - (SETQ |b| 0.0) - (SETQ |temp| 0.0) - ((LAMBDA (|bfVar#13| |i|) - (LOOP - (COND - ((COND ((MINUSP |bfVar#13|) (< |i| 0)) (T (> |i| 0))) - (RETURN NIL)) - ('T - (PROGN - (SETQ |cc| |b|) - (SETQ |b| |temp|) - (SETQ |temp| (+ (- |cc|) (AREF |arr| |i|)))))) - (SETQ |i| (+ |i| |bfVar#13|)))) - (- 1) (+ |n| 1)) - |temp|)))) - -;brutef01(c,z) == -;-- Use series definition. Won't work well if cancellation occurs -; term := 1.0 -; sum := term -; n := 0.0 -; oldsum := 0.0 -; maxnterms := 10000 -; for i in 0..maxnterms until oldsum=sum repeat -; oldsum := sum -; term := term*z/(c+n)/(n+1.0) -; sum := sum + term -; n := n+1.0 -; sum - -(DEFUN |brutef01| (|c| |z|) - (PROG (|maxnterms| |oldsum| |n| |sum| |term|) - (RETURN - (PROGN - (SETQ |term| 1.0) - (SETQ |sum| |term|) - (SETQ |n| 0.0) - (SETQ |oldsum| 0.0) - (SETQ |maxnterms| 10000) - ((LAMBDA (|i| |bfVar#14|) - (LOOP - (COND - ((OR (> |i| |maxnterms|) |bfVar#14|) (RETURN NIL)) - ('T - (PROGN - (SETQ |oldsum| |sum|) - (SETQ |term| - (/ (/ (* |term| |z|) (+ |c| |n|)) (+ |n| 1.0))) - (SETQ |sum| (+ |sum| |term|)) - (SETQ |n| (+ |n| 1.0))))) - (SETQ |i| (+ |i| 1)) - (SETQ |bfVar#14| (EQUAL |oldsum| |sum|)))) - 0 NIL) - |sum|)))) - -;f01(c,z)== -; if negintp(c) -; then -; FloatError('"0F1 not defined for negative integer parameter value ~S",c) -;-- conditions when we'll permit the computation -; else if ABS(c)<1000.0 and ABS(z)<1000.0 -; then -; brutef01(c,z) -; else if ZEROP IMAGPART(z) and ZEROP IMAGPART(c) and z>=0.0 and c>=0.0 -; then -; brutef01(c,z) -;--- else -;--- t := SQRT(-z) -;--- c1 := c-1.0 -;--- p := PHASE(c) -;--- if ABS(c)>10.0*ABS(t) and p>=0.0 and PHASE(c)<.90*dfPi -;--- then BesselJAsymptOrder(c1,2*t)*cgamma(c/(t**(c1))) -;--- else if ABS(t)>10.0*ABS(c) and ABS(t)>50.0 -;--- then BesselJAsympt(c1,2*t)*cgamma(c/(t**(c1))) -;--- else -;--- FloatError('"0F1 not implemented for ~S",[c,z]) -; else if (10.0*ABS(c)>ABS(z)) and ABS(c)<1.0E4 and ABS(z)<1.0E4 -; then -; brutef01(c,z) -; else -; FloatError('"0F1 not implemented for ~S",[c,z]) - -(DEFUN |f01| (|c| |z|) - (PROG () - (RETURN - (COND - ((|negintp| |c|) - (FloatError - "0F1 not defined for negative integer parameter value ~S" - |c|)) - ((AND (< (ABS |c|) 1000.0) (< (ABS |z|) 1000.0)) - (|brutef01| |c| |z|)) - ((AND (ZEROP (IMAGPART |z|)) (ZEROP (IMAGPART |c|)) - (NOT (< |z| 0.0)) (NOT (< |c| 0.0))) - (|brutef01| |c| |z|)) - ((AND (< (ABS |z|) (* 10.0 (ABS |c|))) (< (ABS |c|) 10000.0) - (< (ABS |z|) 10000.0)) - (|brutef01| |c| |z|)) - ('T (FloatError "0F1 not implemented for ~S" (LIST |c| |z|))))))) - -;--- c parameter to 0F1 -;--- w scale factor so that 0 |i| |n|) (RETURN NIL)) - ('T - (PROGN - (SETQ |c| |b|) - (SETQ |b| |temp|) - (SETQ |temp| - (+ (- (* |y| |b|) |c|) (ELT |coefarr| |i|)))))) - (SETQ |i| (+ |i| 1)))) - 1) - (/ (- |temp| |c|) 2))))) - -;--- If plist is a list of coefficients for the Chebychev approximation -;--- of a function f(x), then chebderiveval computes the Chebychev approximation -;--- of f'(x). See Luke, "Special Functions and their approximations, vol. 1 -;--- Academic Press 1969., p. 329 (from Clenshaw and Cooper) - -;--- < definition to be supplied> - -;--- chebstareval(plist,n) computes a Chebychev approximation from a -;--- coefficient list, using shifted Chebychev polynomials of the first kind -;--- The defining relation is that T*(n,x) = T(n,2*x-1). Thus the interval -;--- [0,1] of T*n is the interval [-1,1] of Tn. - -;chebstareval(coeflist,x) == -- evaluation of T*(n,x) -; b := 0 -; temp := 0 -; y := 2*(2*x-1) -; -; for el in coeflist repeat -; c := b; -; b := temp -; temp := y*b -c + el -; temp - y*b/2 - -(defun |chebstareval| (coeflist x) - (let (c y (temp 0) (b 0)) - (setq y (* 2 (- (* 2 x) 1))) - (loop for el in coeflist do - (setq c b) - (setq b temp) - (setq temp (+ (- (* y b) c) el))) - (- temp (/ (* y b) 2)))) - -;--Float definitions for Bessel functions I and J. -;--External references: cgamma, rgamma, chebf01coefmake, chebevalstarsf -;-- floatutils - -;BesselIAsympt(v,z,n) == -; i := COMPLEX(0.0, 1.0) -; if (REALPART(z) = 0.0) -; then return EXPT(i,v)*BesselJ(v,-IMAGPART(z)) -; sum1 := 0.0 -; sum2 := 0.0 -; fourvsq := 4.0*v**2 -; two := 2.0 -; eight := 8.0 -; term1 := 1.0 -;--- sum1, sum2, fourvsq,two,i,eight,term1]) -; for r in 1..n repeat -; term1 := -term1 *(fourvsq-(two*float(r)-1.0)**2)/_ -; (float(r)*eight*z) -; sum1 := sum1 + term1 -; sum2 := sum2 + ABS(term1) -; sqrttwopiz := SQRT(two*dfPi*z) -; EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_ -; EXP(-(float(n)+.5)*dfPi*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2) - - -;defun{BesselIAsympt}{Asymptotic series for BesselJ} -;Asymptotic series for I. See Whittaker, p. 373. -;This is valid for $\frac{-3}{2}\pi < arg z < \frac{1}{2}\pi$ -(DEFUN |BesselIAsympt| (|v| |z| |n|) - (PROG (|sqrttwopiz| |term1| |eight| |two| |fourvsq| |sum2| |sum1| - |i|) - (RETURN - (PROGN - (SETQ |i| (COMPLEX 0.0 1.0)) - (COND - ((EQUAL (REALPART |z|) 0.0) - (RETURN - (* (EXPT |i| |v|) (BesselJ |v| (- (IMAGPART |z|))))))) - (SETQ |sum1| 0.0) - (SETQ |sum2| 0.0) - (SETQ |fourvsq| (* 4.0 (EXPT |v| 2))) - (SETQ |two| 2.0) - (SETQ |eight| 8.0) - (SETQ |term1| 1.0) - ((LAMBDA (|r|) - (LOOP - (COND - ((> |r| |n|) (RETURN NIL)) - ('T - (PROGN - (SETQ |term1| - (- (/ (* |term1| - (- |fourvsq| - (EXPT - (- (* |two| (float |r|)) 1.0) 2))) - (* (* (float |r|) |eight|) |z|)))) - (SETQ |sum1| (+ |sum1| |term1|)) - (SETQ |sum2| (+ |sum2| (ABS |term1|)))))) - (SETQ |r| (+ |r| 1)))) - 1) - (SETQ |sqrttwopiz| (SQRT (* (* |two| Pi) |z|))) - (+ (* (/ (EXP |z|) |sqrttwopiz|) (+ 1.0 |sum1|)) - (* (/ (* (EXP (- (* (* (+ (float |n|) 0.5) Pi) |i|))) - (EXP (- |z|))) - |sqrttwopiz|) - (+ 1.0 |sum2|))))))) - -;--- See Olver, p. 376-382. -;BesselIAsymptOrder(v,vz) == -; z := vz/v -; Pi := dfPi -;--- Use reflection formula (Atlas, p. 492) if v not in right half plane; Is this always accurate? -; if REALPART(v)<0.0 -; then return BesselIAsymptOrder(-v,vz) + 2.0/Pi*SIN(-v*Pi)*BesselKAsymptOrder(-v,vz) -;--- Use the reflection formula (Atlas, p. 496) if z not in right half plane; -; if REALPART(vz) < 0.0 -; then return EXPT(-1.0,v)*BesselIAsymptOrder(v,-vz) -; vinv := 1.0/v -; opzsqroh := SQRT(1.0+z*z) -; eta := opzsqroh + LOG(z/(1.0+opzsqroh)) -; p := 1.0/opzsqroh -; p2 := p*p -; p4 := p2*p2 -; u0p := 1. -; u1p := 1.0/8.0*p-5.0/24.0*p*p2 -; u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2 -; u3p := (75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_ -; *p2)*p2)*p2*p -; u4p := (3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4 -; u5p := (59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p -; hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv) -; EXP(v*eta)/(SQRT(2.0*Pi*v)*SQRT(opzsqroh))*hornerresult - -(DEFUN |BesselIAsymptOrder| (|v| |vz|) - (PROG (|hornerresult| |u5p| |u4p| |u3p| |u2p| |u1p| |u0p| |p4| |p2| - |p| |eta| |opzsqroh| |vinv| |z|) - (RETURN - (PROGN - (SETQ |z| (/ |vz| |v|)) - (COND - ((< (REALPART |v|) 0.0) - (RETURN - (+ (|BesselIAsymptOrder| (- |v|) |vz|) - (* (* (/ 2.0 Pi) (SIN (- (* |v| Pi)))) - (|BesselKAsymptOrder| (- |v|) |vz|)))))) - (COND - ((< (REALPART |vz|) 0.0) - (RETURN - (* (EXPT (- 1.0) |v|) (|BesselIAsymptOrder| |v| (- |vz|)))))) - (SETQ |vinv| (/ 1.0 |v|)) - (SETQ |opzsqroh| (SQRT (+ 1.0 (* |z| |z|)))) - (SETQ |eta| (+ |opzsqroh| (LOG (/ |z| (+ 1.0 |opzsqroh|))))) - (SETQ |p| (/ 1.0 |opzsqroh|)) - (SETQ |p2| (* |p| |p|)) - (SETQ |p4| (* |p2| |p2|)) - (SETQ |u0p| 1.0) - (SETQ |u1p| - (- (* (/ 1.0 8.0) |p|) (* (* (/ 5.0 24.0) |p|) |p2|))) - (SETQ |u2p| - (* (+ (/ 9.0 128.0) - (* (+ (- (/ 77.0 192.0)) (* (/ 385.0 1152.0) |p2|)) - |p2|)) - |p2|)) - (SETQ |u3p| - (* (* (+ (/ 75.0 1024.0) - (* (+ (- (/ 4563.0 5120.0)) - (* (- (/ 17017.0 9216.0) - (* (/ 85085.0 82944.0) |p2|)) - |p2|)) - |p2|)) - |p2|) - |p|)) - (SETQ |u4p| - (* (+ (/ 3675.0 32768.0) - (* (+ (- (/ 96833.0 40960.0)) - (* (+ (/ 144001.0 16384.0) - (* (+ (- (/ 7436429.0 663552.0)) - (* (/ 3.7182145E7 7962624.0) |p2|)) - |p2|)) - |p2|)) - |p2|)) - |p4|)) - (SETQ |u5p| - (* (* (+ (/ 59535.0 262144.0) - (* (+ (- (/ 6.7608983E7 9175040.0)) - (* (+ (/ 2.50881631E8 5898240.0) - (* - (+ (- (/ 1.08313205E8 1179648.0)) - (* - (- (/ 5.391411025E9 6.3700992E7) - (* - (/ 5.391411025E9 1.91102976E8) - |p2|)) - |p2|)) - |p2|)) - |p2|)) - |p2|)) - |p4|) - |p|)) - (SETQ |hornerresult| - (horner (LIST |u5p| |u4p| |u3p| |u2p| |u1p| |u0p|) - |vinv|)) - (* (/ (EXP (* |v| |eta|)) - (* (SQRT (* (* 2.0 Pi) |v|)) (SQRT |opzsqroh|))) - |hornerresult|))))) - -;---See also Olver, pp. 376-382 -;BesselKAsymptOrder (v,vz) == -; z := vz/v -; vinv := 1.0/v -; opzsqroh := SQRT(1.0+z*z) -; eta := opzsqroh + LOG(z/(1.0+opzsqroh)) -; p := 1.0/opzsqroh -; p2 := p**2 -; p4 := p2**2 -; u0p := 1. -; u1p := (1.0/8.0*p-5.0/24.0*p**3)*(-1.0) -; u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2 -; u3p := ((75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_ -; *p2)*p2)*p2*p)*(-1.0) -; u4p := (3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4 -; u5p := ((59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p)*(-1.0) -; hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv) -; SQRT(dfPi/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult - -(DEFUN |BesselKAsymptOrder| (|v| |vz|) - (PROG (|hornerresult| |u5p| |u4p| |u3p| |u2p| |u1p| |u0p| |p4| |p2| - |p| |eta| |opzsqroh| |vinv| |z|) - (RETURN - (PROGN - (SETQ |z| (/ |vz| |v|)) - (SETQ |vinv| (/ 1.0 |v|)) - (SETQ |opzsqroh| (SQRT (+ 1.0 (* |z| |z|)))) - (SETQ |eta| (+ |opzsqroh| (LOG (/ |z| (+ 1.0 |opzsqroh|))))) - (SETQ |p| (/ 1.0 |opzsqroh|)) - (SETQ |p2| (EXPT |p| 2)) - (SETQ |p4| (EXPT |p2| 2)) - (SETQ |u0p| 1.0) - (SETQ |u1p| - (* (- (* (/ 1.0 8.0) |p|) (* (/ 5.0 24.0) (EXPT |p| 3))) - (- 1.0))) - (SETQ |u2p| - (* (+ (/ 9.0 128.0) - (* (+ (- (/ 77.0 192.0)) (* (/ 385.0 1152.0) |p2|)) - |p2|)) - |p2|)) - (SETQ |u3p| - (* (* (* (+ (/ 75.0 1024.0) - (* (+ (- (/ 4563.0 5120.0)) - (* (- (/ 17017.0 9216.0) - (* (/ 85085.0 82944.0) |p2|)) - |p2|)) - |p2|)) - |p2|) - |p|) - (- 1.0))) - (SETQ |u4p| - (* (+ (/ 3675.0 32768.0) - (* (+ (- (/ 96833.0 40960.0)) - (* (+ (/ 144001.0 16384.0) - (* (+ (- (/ 7436429.0 663552.0)) - (* (/ 3.7182145E7 7962624.0) |p2|)) - |p2|)) - |p2|)) - |p2|)) - |p4|)) - (SETQ |u5p| - (* (* (* (+ (/ 59535.0 262144.0) - (* (+ (- (/ 6.7608983E7 9175040.0)) - (* (+ (/ 2.50881631E8 5898240.0) - (* - (+ (- (/ 1.08313205E8 1179648.0)) - (* - (- (/ 5.391411025E9 6.3700992E7) - (* - (/ 5.391411025E9 1.91102976E8) - |p2|)) - |p2|)) - |p2|)) - |p2|)) - |p2|)) - |p4|) - |p|) - (- 1.0))) - (SETQ |hornerresult| - (horner (LIST |u5p| |u4p| |u3p| |u2p| |u1p| |u0p|) - |vinv|)) - (* (/ (* (SQRT (/ Pi (* 2.0 |v|))) (EXP (- (* |v| |eta|)))) - (SQRT |opzsqroh|)) - |hornerresult|))))) - -\end{chunk} -\eject -\begin{thebibliography}{99} -\bibitem{1} nothing -\end{thebibliography} -\end{document} -- 1.7.5.4