Imported Upstream version 16.6

1 files changed, 47 insertions, 52 deletions

diff --git a/numerics.scm b/numerics.scm
index f0aaa71..7560354 100644
--- a/numerics.scm
+++ b/numerics.scm
@@ -76,8 +76,7 @@
 ;;; 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)
-  (if (or (< m 0) 
-	  (> m L) 
+  (if (or (not (<= 0 m L)) 
 	  (> (abs x) 1.0))
       (snd-error "invalid arguments to plgndr")
       (let ((pmm 1.0)
@@ -706,58 +705,55 @@
   (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)
+    (let ((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! 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)
+			  (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))))))))))))))
+	  (eps 1e-17)
 	  (s 0.0))
       (do ((k 0 (+ k 1)))
 	  ((= k id))
 	(let* ((ak (+ (* 8 k) m))
-	       (p (- id k))
-	       (t (expm p ak)))
+	       (t (expm (- id k) ak)))
 	  (set! s (+ s (/ t ak)))
 	  (set! s (- s (floor s)))))
       
@@ -765,8 +761,7 @@
       (let ((happy #f))
 	(do ((k id (+ k 1)))
 	    ((or (> k (+ id 100)) happy) s)
-	  (let* ((ak (+ (* 8 k) m))
-		 (t (/ (expt 16.0 (- id k)) ak)))
+	  (let ((t (/ (expt 16.0 (- id k)) (+ (* 8 k) m))))
 	    (set! happy (< t eps))
 	    (set! s (+ s t))
 	    (set! s (- s (floor s))))))))
@@ -778,7 +773,7 @@
 	 (s3 (series 5 id))
 	 (s4 (series 6 id))
 	 (pid (- (+ (* 4.0 s1) (* -2.0 s2)) s3 s4)))
-    (set! pid (+ 1.0 (- pid (floor pid))))
+    (set! pid (- (+ 1.0 pid) (floor pid)))
     (ihex pid 10 chx)
     (format #t " position = ~D~% fraction = ~,15F~% hex digits =  ~S~%" id pid chx)))