Imported Upstream version 16.5~dfsg

1 files changed, 85 insertions, 92 deletions

diff --git a/numerics.scm b/numerics.scm
index e058304..f0aaa71 100644
--- a/numerics.scm
+++ b/numerics.scm
@@ -74,10 +74,10 @@
 
 ;;; --------------------------------------------------------------------------------
 ;;; from Numerical Recipes
-(define (plgndr l m x)			;Legendre polynomial P m/l (x), m and l integer
-					;0 <= m <= l and -1<= x <= 1 (x real)
+(define (plgndr L m x)			;Legendre polynomial P m/L (x), m and L integer
+					;0 <= m <= L and -1<= x <= 1 (x real)
   (if (or (< m 0) 
-	  (> m l) 
+	  (> m L) 
 	  (> (abs x) 1.0))
       (snd-error "invalid arguments to plgndr")
       (let ((pmm 1.0)
@@ -91,14 +91,14 @@
 		  ((> i m))
 		(set! pmm (* (- pmm) fact somx2))
 		(set! fact (+ fact 2.0)))))
-	(if (= l m) 
+	(if (= L m) 
 	    pmm
 	    (let ((pmmp1 (* x pmm (+ (* 2 m) 1))))
-	      (if (= l (+ m 1)) 
+	      (if (= L (+ m 1)) 
 		  pmmp1
 		  (let ((pk 0.0)) ; NR used "ll" which is unreadable
 		    (do ((k (+ m 2) (+ k 1)))
-			((> k l))
+			((> k L))
 		      (set! pk (/ (- (* x (- (* 2 k) 1) pmmp1) 
 				      (* (+ k m -1) pmm)) 
 				   (- k m)))
@@ -154,30 +154,23 @@
 
 
 (define* (gegenbauer n x (alpha 0.0))
-  (if (< alpha -0.5) (set! alpha -0.5))
-  (if (= n 0)
-      1.0
-      (if (= alpha 0.0) ; maxima and A&S 22.3.14 (gsl has bogus values here)
-	  (* (/ 2.0 n) 
-	     (cos (* n x)))
-	  (if (= n 1)   ; gsl splits out special cases 
-	      (* 2 alpha x)                             ; G&R 8.93(2)
-	      (if (= n 2)
-		  (- (* 2 alpha (+ alpha 1) x x) alpha) ; G&R 8.93(3)
-		  (let ((fn1 (* 2 x alpha))
-			(fn 0.0)
-			(fn2 1.0))
-		    (if (= n 1)
-			fn1
-			(do ((k 2 (+ k 1))
-			     (k0 2.0 (+ k0 1.0)))
-			    ((> k n) fn)
-			  (set! fn (/ (- (* 2 x fn1 (+ k alpha -1.0))
-					 (* fn2 (+ k (* 2 alpha) -2.0)))
-				      k0))
-			  (set! fn2 fn1)
-			  (set! fn1 fn)))))))))
-
+  (set! alpha (max alpha -0.5))
+  (cond ((= n 0)       1.0)
+	((= alpha 0.0) (* (/ 2.0 n) (cos (* n x))))           ; maxima and A&S 22.3.14 (gsl has bogus values here)
+	((= n 1)       (* 2 alpha x))                         ; G&R 8.93(2)
+	((= n 2)       (- (* 2 alpha (+ alpha 1) x x) alpha)) ; G&R 8.93(3)
+	(else
+	 (let ((fn1 (* 2 x alpha))
+	       (fn 0.0)
+	       (fn2 1.0))
+	   (if (= n 1)
+	       fn1
+	       (do ((k 2 (+ k 1))
+		    (k0 2.0 (+ k0 1.0)))
+		   ((> k n) fn)
+		 (set! fn (/ (- (* 2 x fn1 (+ k alpha -1.0)) (* fn2 (+ k (* 2 alpha) -2.0))) k0))
+		 (set! fn2 fn1)
+		 (set! fn1 fn)))))))
 
 ;;; (with-sound (:scaled-to 0.5) (do ((i 0 (+ i 1)) (x 0.0 (+ x .1))) ((= i 10000)) (outa i (gegenbauer 15 (cos x) 1.0))))
 
@@ -637,21 +630,20 @@
 	   
 	  ((or (= (modulo n 2) 0) (= (modulo n 3) 0) (= (modulo n 5) 0) (= (modulo n 7) 0) 
 	       (= (modulo n 17) 0) (= (modulo n 13) 0) (= (modulo n 257) 0) (= (modulo n 11) 0))
-	   (let ((divisor (if (= (modulo n 2) 0) 2
-			      (if (= (modulo n 3) 0) 3
-				  (if (= (modulo n 5) 0) 5
-				      (if (= (modulo n 7) 0) 7
-					  (if (= (modulo n 17) 0) 17
-					      (if (= (modulo n 13) 0) 13
-						  (if (= (modulo n 11) 0) 11
-						      257)))))))))
-	     (let ((val (sin-m*pi/n 1 (/ n divisor))))
-	       (and val
-		    `(let ((ex ,val))
-		       (/ (- (expt (+ (sqrt (- 1 (* ex ex))) (* 0+i ex)) (/ 1 ,divisor))
-			     (expt (- (sqrt (- 1 (* ex ex))) (* 0+i ex)) (/ 1 ,divisor)))
-			  0+2i))))))
-	   
+	   (let* ((divisor (cond ((= (modulo n 2) 0) 2)
+				 ((= (modulo n 3) 0) 3)
+				 ((= (modulo n 5) 0) 5)
+				 ((= (modulo n 7) 0) 7)
+				 ((= (modulo n 17) 0) 17)
+				 ((= (modulo n 13) 0) 13)
+				 ((= (modulo n 11) 0) 11)
+				 (else 257)))
+		  (val (sin-m*pi/n 1 (/ n divisor))))
+	     (and val
+		  `(let ((ex ,val))
+		     (/ (- (expt (+ (sqrt (- 1 (* ex ex))) (* 0+i ex)) (/ 1 ,divisor))
+			   (expt (- (sqrt (- 1 (* ex ex))) (* 0+i ex)) (/ 1 ,divisor)))
+			0+2i)))))
 	  (else #f))))
 
 #|
@@ -671,7 +663,7 @@
 		  (begin
 		    (set! maxerr err)
 		    (set! max-case (/ m n)))))))))
-  (format #t "sin-m*pi/n (~A cases) max err ~A at ~A~%" cases maxerr max-case))
+  (format () "sin-m*pi/n (~A cases) max err ~A at ~A~%" cases maxerr max-case))
 
 :(sin (/ pi (* 257 17)))
 0.00071906440440859
@@ -711,53 +703,54 @@
 	(set! (chx i) (hx (floor y))))
       chx))
   
-  (define expm
-    (let* ((ntp 25)
-	   (tp1 0)
-	   (tp (make-vector ntp)))
-      (lambda (p ak)
-	;; expm = 16^p mod ak.  This routine uses the left-to-right binary exponentiation scheme.
-	
-	;; If this is the first call to expm, fill the power of two table tp.
-	(if (= tp1 0)
-	    (begin
-	      (set! tp1 1)
-	      (set! (tp 0) 1.0)
-	      (do ((i 1 (+ i 1)))
-		  ((= i ntp))
-		(set! (tp i) (* 2.0 (tp (- i 1)))))))
-	
-	(if (= ak 1.0)
-	    0.0
-	    (let ((pl -1))
-	      ;;  Find the greatest power of two less than or equal to p.
-	      (do ((i 0 (+ i 1)))
-		  ((or (not (= pl -1)) 
-		       (= i ntp)))
-		(if (> (tp i) p)
-		    (set! pl i)))
-	      
-	      (if (= pl -1) (set! pl ntp))
-	      (let ((pt (tp (- pl 1)))
-		    (p1 p)
-		    (r 1.0))
-		;;  Perform binary exponentiation algorithm modulo ak.
-		
-		(do ((j 1 (+ j 1)))
-		    ((> j pl) r)
-		  (if (>= p1 pt)
-		      (begin
-			(set! r (* 16.0 r))
-			(set! r (- r (* ak (floor (/ r ak)))))
-			(set! p1 (- p1 pt))))
-		  (set! pt (* 0.5 pt))
-		  (if (>= pt 1.0)
-		      (begin
-			(set! r (* r r))
-			(set! r (- r (* ak (floor (/ r ak))))))))))))))
-  
   (define (series m id)
     ;; This routine evaluates the series  sum_k 16^(id-k)/(8*k+m) using the modular exponentiation technique.
+    
+    (define expm
+      (let* ((ntp 25)
+	     (tp1 0)
+	     (tp (make-vector ntp)))
+	(lambda (p ak)
+	  ;; expm = 16^p mod ak.  This routine uses the left-to-right binary exponentiation scheme.
+	  
+	  ;; If this is the first call to expm, fill the power of two table tp.
+	  (if (= tp1 0)
+	      (begin
+		(set! tp1 1)
+		(set! (tp 0) 1.0)
+		(do ((i 1 (+ i 1)))
+		    ((= i ntp))
+		  (set! (tp i) (* 2.0 (tp (- i 1)))))))
+	  
+	  (if (= ak 1.0)
+	      0.0
+	      (let ((pl -1))
+		;;  Find the greatest power of two less than or equal to p.
+		(do ((i 0 (+ i 1)))
+		    ((or (not (= pl -1)) 
+			 (= i ntp)))
+		  (if (> (tp i) p)
+		      (set! pl i)))
+		
+		(if (= pl -1) (set! pl ntp))
+		(let ((pt (tp (- pl 1)))
+		      (p1 p)
+		      (r 1.0))
+		  ;;  Perform binary exponentiation algorithm modulo ak.
+		  
+		  (do ((j 1 (+ j 1)))
+		      ((> j pl) r)
+		    (if (>= p1 pt)
+			(begin
+			  (set! r (* 16.0 r))
+			  (set! r (- r (* ak (floor (/ r ak)))))
+			  (set! p1 (- p1 pt))))
+		    (set! pt (* 0.5 pt))
+		    (if (>= pt 1.0)
+			(begin
+			  (set! r (* r r))
+			  (set! r (- r (* ak (floor (/ r ak))))))))))))))
+    
     (let ((eps 1e-17)
 	  (s 0.0))
       (do ((k 0 (+ k 1)))
@@ -784,7 +777,7 @@
 	 (s2 (series 4 id))
 	 (s3 (series 5 id))
 	 (s4 (series 6 id))
-	 (pid (+ (* 4.0 s1) (* -2.0 s2) (- s3) (- s4))))
+	 (pid (- (+ (* 4.0 s1) (* -2.0 s2)) s3 s4)))
     (set! pid (+ 1.0 (- pid (floor pid))))
     (ihex pid 10 chx)
     (format #t " position = ~D~% fraction = ~,15F~% hex digits =  ~S~%" id pid chx)))