separate bigfloat-fresnel into it's own file

Author: bdeket

math-lib/math/private/bigfloat/bigfloat-fresnel.rkt

@@ -0,0 +1,65 @@
+#lang typed/racket/base
+
+;Bigfloat implementation:
+;  adaptation of powerseries as shown on wikipedia
+;  this implementation is potentially exact
+;  but for large x (> 5!) it is really slow and needs a lot of bits in bf-precision (~2x²)!
+
+
+(require "bigfloat-struct.rkt")
+
+(provide bfFresnel-S bfFresnel-RS bfFresnel-C bfFresnel-RC)
+
+(define (precision-check [a : Bigfloat][maxp : Integer]) : Integer
+  (define p (bf-precision))
+  (define a2 (expt (bigfloat->real a) 2))
+  (define min-precision-needed (* 2 a2))
+  (define expected-precision-loss (* 4/5 a2))
+  (define mp (+ 5 (round (inexact->exact (max min-precision-needed (+ p expected-precision-loss))))))
+  (cond
+    [(<= mp maxp) mp]
+    [else
+     (error (format "bfFresnel: calculation aborted
+ Minimum precision needed for calculating ~a... is ~a
+ This is more than the maximum allowed calculating precision ~a->~a."
+                    (bigfloat->flonum a) mp p maxp))]))
+(define (bfFresnel-RS [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
+  (define p (bf-precision))
+  (bfcopy
+   (parameterize ([bf-precision (precision-check x maxp)])
+     (define X3 (bfexpt x (bf 3)))
+     (define X4 (bfexpt x (bf 4)))
+     (define prsn (bfexpt (bf 1/2) (bf p)))
+     (define-values (s l)
+       (for/fold : (Values Bigfloat Bigfloat)
+         ([s : Bigfloat (bf/ X3 (bf 3))]
+          [l : Bigfloat (bf/ X3 (bf 3))])
+         ([n (in-naturals 1)]
+          #:break (bf< (bfabs (bf/ l s)) prsn))
+         (define l+ (bf* (bf -1) X4 (bf* l (bf (/ (- (* 4 n) 1)
+                                                  (* 2 n)(+ (* 2 n) 1)(+ (* 4 n) 3))))))
+         (values (bf+ s l+) l+)))
+     s)))
+(define (bfFresnel-S [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
+  (bf* (bfsqrt (bf/ (bf 2) pi.bf))
+       (bfFresnel-RS (bf* (bfsqrt (bf/ pi.bf (bf 2))) x) maxp)))
+  
+(define (bfFresnel-RC [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
+  (define p (bf-precision))
+  (bfcopy
+   (parameterize ([bf-precision (precision-check x maxp)])
+     (define X4 (bfexpt x (bf 4)))
+     (define prsn (bfexpt (bf 1/2) (bf (bf-precision))))
+     (define-values (s l)
+       (for/fold : (Values Bigfloat Bigfloat)
+         ([s : Bigfloat x]
+          [l : Bigfloat x])
+         ([n (in-naturals 1)]
+          #:break (bf< (bfabs (bf/ l s)) prsn))
+         (define l+ (bf* (bf -1) X4 (bf* l (bf (/ (- (* 4 n) 3)
+                                                  (* 2 n)(- (* 2 n) 1)(+ (* 4 n) 1))))))
+         (values (bf+ s l+) l+)))
+     s)))
+(define (bfFresnel-C [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
+  (bf* (bfsqrt (bf/ (bf 2) pi.bf))
+       (bfFresnel-RC (bf* (bfsqrt (bf/ pi.bf (bf 2))) x) maxp)))

math-lib/math/private/functions/fresnel.rkt

@@ -12,10 +12,6 @@ floating point implementation:
   adaptation of the algorithm in the fresnl.c file from Cephes,
   but with the limit for the rational powerseries lowered to x<1.5625 (coming from 2.5625)
   above x>1.6 the error becomes too big (~1e-8 x~2.5)
-Bigfloat implementation:
-  adaptation of powerseries as shown on wikipedia
-
-would like to extend this together with erf to the complex plane
 
 |#
 
@@ -204,73 +200,13 @@ would like to extend this together with erf to the complex plane
     [else (* (sqrt (/ pi 2))(complex-Fresnel-C (* (sqrt (/ 2 pi)) z)))]))
 
 ;------------------------------
-(module* bfFresnel #f
-  (require math/bigfloat)
-  (provide bfFresnel-S bfFresnel-RS bfFresnel-C bfFresnel-RC)
-  ;this implementation is potentially exact
-  ;but for large x (> 5!) it is really slow and needs a lot of bits in bf-precision (~2x²)!
-  (define (precision-check [a : Bigfloat][maxp : Integer]) : Integer
-    (define p (bf-precision))
-    (define a2 (expt (bigfloat->real a) 2))
-    (define min-precision-needed (* 2 a2))
-    (define expected-precision-loss (* 4/5 a2))
-    (define mp (+ 5 (round (inexact->exact (max min-precision-needed (+ p expected-precision-loss))))))
-    (cond
-      [(<= mp maxp) mp]
-      [else
-       (error (format "bfFresnel: calculation aborted
- Minimum precision needed for calculating ~a... is ~a
- This is more than the maximum allowed calculating precision ~a->~a."
-                      (bigfloat->flonum a) mp p maxp))]))
-  (define (bfFresnel-RS [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
-    (define p (bf-precision))
-    (bfcopy
-     (parameterize ([bf-precision (precision-check x maxp)])
-       (define X3 (bfexpt x (bf 3)))
-       (define X4 (bfexpt x (bf 4)))
-       (define prsn (bfexpt (bf 1/2) (bf p)))
-       (define-values (s l)
-         (for/fold : (Values Bigfloat Bigfloat)
-           ([s : Bigfloat (bf/ X3 (bf 3))]
-            [l : Bigfloat (bf/ X3 (bf 3))])
-           ([n (in-naturals 1)]
-            #:break (bf< (bfabs (bf/ l s)) prsn))
-           (define l+ (bf* (bf -1) X4 (bf* l (bf (/ (- (* 4 n) 1)
-                                                    (* 2 n)(+ (* 2 n) 1)(+ (* 4 n) 3))))))
-           (values (bf+ s l+) l+)))
-       s)))
-  (define (bfFresnel-S [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
-    (bf* (bfsqrt (bf/ (bf 2) pi.bf))
-         (bfFresnel-RS (bf* (bfsqrt (bf/ pi.bf (bf 2))) x) maxp)))
-  
-  (define (bfFresnel-RC [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
-    (define p (bf-precision))
-    (bfcopy
-     (parameterize ([bf-precision (precision-check x maxp)])
-       (define X4 (bfexpt x (bf 4)))
-       (define prsn (bfexpt (bf 1/2) (bf (bf-precision))))
-       (define-values (s l)
-         (for/fold : (Values Bigfloat Bigfloat)
-           ([s : Bigfloat x]
-            [l : Bigfloat x])
-           ([n (in-naturals 1)]
-            #:break (bf< (bfabs (bf/ l s)) prsn))
-           (define l+ (bf* (bf -1) X4 (bf* l (bf (/ (- (* 4 n) 3)
-                                                    (* 2 n)(- (* 2 n) 1)(+ (* 4 n) 1))))))
-           (values (bf+ s l+) l+)))
-       s)))
-  (define (bfFresnel-C [x : Bigfloat][maxp : Integer (* 2 (bf-precision))])
-    (bf* (bfsqrt (bf/ (bf 2) pi.bf))
-         (bfFresnel-RC (bf* (bfsqrt (bf/ pi.bf (bf 2))) x) maxp)))
-    )
-
-#;(module* test #f
+(module* test #f
     ;some functions to check the function.
     ;commented out, because (partly) put in function-tests from math-test
-  (require math/bigfloat
-           typed/rackunit)
+  (require typed/rackunit)
   (require (submod "..")
-           (submod ".." bfFresnel))
+           "../../bigfloat.rkt"
+           "../bigfloat/bigfloat-fresnel.rkt")
 
   (define (test [a : Flonum] #:ε [ε : Flonum 1e-15])
     (check-= (Fresnel-S  a) (bigfloat->flonum (bfFresnel-S  (bf a))) ε (format "(Fresnel-S ~a)" a))