Complex erf, start fresnel

Author: bdeket

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

@@ -198,6 +198,48 @@ IEEE Transactions on Communications, 2000, vol 48, pp 529--532
 
 ;; ===================================================================================================
 
+(: complex-erf (Number -> Number))
+(define 2π (* 2 pi))
+(define (complex-erf z)
+  (cond
+    [(<= (magnitude z) 8)
+     (define nn 32)
+     (define x (real-part z))
+     (define y (imag-part z))
+     (define k1 (* (/ 2 pi) (exp (* -1 x x))))
+     (define k2 (exp (* -i 2 x y)))
+     (define f (+ (erf x)
+                  (if (= x 0)
+                      (* (/ +i pi) y)
+                      (* (/ k1 (* 4 x)) (- 1 k2)))))
+     (cond
+       [(= y 0) f]
+       [else
+        (define s5
+          (for/fold : Number
+            ([s5 : Number 0])
+            ([n (in-range 1 (+ nn 1))])
+            (define s3 (/ (exp (* -1 n n 1/4)) (+ (* n n) (* 4 x x))))
+            (define s4 (- (* 2 x) (* k2 (- (* 2 x (cosh (* n y))) (* +i n (sinh (* n y)))))))
+            (+ s5 (* s3 s4))))
+        (+ f (* k1 s5))])]
+    [else
+     (define neg? (< (real-part z) 0))
+     (define z+ (if neg? (- z) z))
+     (define nmax 193)
+     (define y (* 2 z z))
+     (define s (for/fold : Number
+                 ([s : Number 1])
+                 ([n (in-range nmax 0 -2)])
+                 (- 1 (* n (/ s y)))))
+     (define f (* (if neg? -1 1)
+                  (- 1 (* s (/ (exp (* -1 z z)) (* sqrtpi z))))))
+     (if (= (real-part z) 0)
+         (- f 1)
+         f)]))
+
+;; ===================================================================================================
+
 (: erf (case-> (Zero -> Zero)
                (Flonum -> Flonum)
                (Real -> (U Zero Flonum))))

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

@@ -27,50 +27,52 @@ would like to extend this together with erf to the complex plane
 
 ;------------------------------
 ;polinomials for the rational assymptotic aproximation for S & C
-(define fn (make-flpolyfun ( 3.76329711269987889006E-20
-                             1.34283276233062758925E-16
-                             1.72010743268161828879E-13
-                             1.02304514164907233465E-10
-                             3.05568983790257605827E-8
-                             4.63613749287867322088E-6
-                             3.45017939782574027900E-4
-                             1.15220955073585758835E-2
-                             1.43407919780758885261E-1
-                             4.21543555043677546506E-1)))
-(define fd (make-flpolyfun ( 1.25443237090011264384E-20
-                             4.52001434074129701496E-17
-                             5.88754533621578410010E-14
-                             3.60140029589371370404E-11
-                             1.12699224763999035261E-8
-                             1.84627567348930545870E-6
-                             1.55934409164153020873E-4
-                             6.44051526508858611005E-3
-                             1.16888925859191382142E-1
-                             7.51586398353378947175E-1
-                             1.0)))
-(define gn (make-flpolyfun ( 1.86958710162783235106E-22
-                             8.36354435630677421531E-19
-                             1.37555460633261799868E-15
-                             1.08268041139020870318E-12
-                             4.45344415861750144738E-10
-                             9.82852443688422223854E-8
-                             1.15138826111884280931E-5
-                             6.84079380915393090172E-4
-                             1.87648584092575249293E-2
-                             1.97102833525523411709E-1
-                             5.04442073643383265887E-1)))
-(define gd (make-flpolyfun ( 1.86958710162783236342E-22
-                             8.39158816283118707363E-19
-                             1.38796531259578871258E-15
-                             1.10273215066240270757E-12
-                             4.60680728146520428211E-10
-                             1.04314589657571990585E-7
-                             1.27545075667729118702E-5
-                             8.14679107184306179049E-4
-                             2.53603741420338795122E-2
-                             3.37748989120019970451E-1
-                             1.47495759925128324529E0
-                             1.0)))
+(define fn/d (make-quotient-flpolyfun
+              ( 3.76329711269987889006E-20
+                1.34283276233062758925E-16
+                1.72010743268161828879E-13
+                1.02304514164907233465E-10
+                3.05568983790257605827E-8
+                4.63613749287867322088E-6
+                3.45017939782574027900E-4
+                1.15220955073585758835E-2
+                1.43407919780758885261E-1
+                4.21543555043677546506E-1)
+              ( 1.25443237090011264384E-20
+                4.52001434074129701496E-17
+                5.88754533621578410010E-14
+                3.60140029589371370404E-11
+                1.12699224763999035261E-8
+                1.84627567348930545870E-6
+                1.55934409164153020873E-4
+                6.44051526508858611005E-3
+                1.16888925859191382142E-1
+                7.51586398353378947175E-1
+                1.0)))
+(define gn/d (make-quotient-flpolyfun
+              ( 1.86958710162783235106E-22
+                8.36354435630677421531E-19
+                1.37555460633261799868E-15
+                1.08268041139020870318E-12
+                4.45344415861750144738E-10
+                9.82852443688422223854E-8
+                1.15138826111884280931E-5
+                6.84079380915393090172E-4
+                1.87648584092575249293E-2
+                1.97102833525523411709E-1
+                5.04442073643383265887E-1)
+              ( 1.86958710162783236342E-22
+                8.39158816283118707363E-19
+                1.38796531259578871258E-15
+                1.10273215066240270757E-12
+                4.60680728146520428211E-10
+                1.04314589657571990585E-7
+                1.27545075667729118702E-5
+                8.14679107184306179049E-4
+                2.53603741420338795122E-2
+                3.37748989120019970451E-1
+                1.47495759925128324529E0
+                1.0)))
 
 ;------------------------------
 ;polinomials for the rational powerseries aproximation for S
@@ -100,8 +102,8 @@ would like to extend this together with erf to the complex plane
      (define t (fl* pi (fl* x x)))
      (define t/2 (fl/ t 2.0))
      (define U (fl/ 1.0 (fl* t t)))
-     (define f (fl- 1.0 (fl* U (fl/ (fn U)(fd U)))))
-     (define g (fl/ (fl/ (gn U)(gd U)) t))
+     (define f (fl- 1.0 (fl* U (fn/d U))))
+     (define g (fl/ (gn/d U) t))
      (fl- 0.5 (fl/ (fl+ (fl* f (flcos t/2))(fl* g (flsin t/2)))
                    (fl* pi x)))]))
 
@@ -133,8 +135,8 @@ would like to extend this together with erf to the complex plane
      (define t (fl* pi (fl* x x)))
      (define t/2 (fl/ t 2.0))
      (define U (fl/ 1.0 (fl* t t)))
-     (define f (fl- 1.0 (fl* U (fl/ (fn U)(fd U)))))
-     (define g (fl/ (fl/ (gn U)(gd U)) t))
+     (define f (fl- 1.0 (fl* U (fn/d U))))
+     (define g (fl/ (gn/d U) t))
      (fl+ 0.5 (fl/ (fl- (fl* f (flsin t/2))(fl* g (flcos t/2)))
                    (fl* pi x)))]))
 
@@ -207,7 +209,7 @@ would like to extend this together with erf to the complex plane
          (bfFresnel-RC (bf* (bfsqrt (bf/ pi.bf (bf 2))) x) maxp)))
     )
 
-#;(module* test racket/base
+(module* test racket/base
     ;some functions to check the function.
     ;commented out, because (partly) put in function-tests from math-test
   (require math/bigfloat