Complex erf and Fresnel functions

Author: bdeket

math-doc/math/scribblings/math-special-functions.scrbl

@@ -19,7 +19,7 @@
 The term ``special function'' has no formal definition. However, for the purposes of the
 @racketmodname[math] library, a @deftech{special function} is one that is not @tech{elementary}.
 
-The special functions are split into three groups: @secref{complex-functions}, @secref{real-functions} and
+The special functions are split into two groups: @secref{complex-real-functions} and
 @secref{flonum-functions}. Functions that accept real arguments are usually defined
 in terms of their flonum counterparts, but are different in two crucial ways:
 @itemlist[
@@ -56,7 +56,7 @@ The most general type @racket[Real -> (U Zero Flonum)] is used to generate
 @racket[lambert]'s contract when it is used in untyped code. Except for this discussion,
 this the only type documented for @racket[lambert].
 
-@section[#:tag "complex-functions"]{Complex Functions}
+@section[#:tag "complex-real-functions"]{Complex and real Functions}
 
 @defproc[(gamma [x Number]) Number]{
 Computes the @hyperlink["http://en.wikipedia.org/wiki/Gamma_function"]{gamma function},
@@ -132,8 +132,6 @@ error.
 If you need low relative error near negative roots, use @racket[bfpsi0].
 }
 
-@section[#:tag "real-functions"]{Real Functions}
-
 @defproc[(psi [m Integer] [x Real]) Flonum]{
 Computes a @hyperlink["http://en.wikipedia.org/wiki/Polygamma_function"]{polygamma function},
 or the @racket[m]th logarithmic derivative of the gamma function. The order @racket[m] must be
@@ -154,7 +152,7 @@ near negative roots. Near negative roots, relative error is apparently unbounded
 is low.
 }
 
-@deftogether[(@defproc[(erf [x Real]) Real]
+@deftogether[(@defproc[(erf [x Number]) Number]
               @defproc[(erfc [x Real]) Real])]{
 Compute the @hyperlink["http://en.wikipedia.org/wiki/Error_function"]{error function and
 complementary error function}, respectively. The only exact cases are @racket[(erf 0) = 0]
@@ -167,7 +165,8 @@ and @racket[(erfc 0) = 1].
                  (erf 1)
                  (- 1 (erfc 1))
                  (erf -1)
-                 (- (erfc 1) 1)]
+                 (- (erfc 1) 1)
+                 (erf 1+i)]
 
 Mathematically, erfc(@italic{x}) = 1 - erf(@italic{x}), but having separate implementations
 can help maintain accuracy. To compute an expression containing erf, use @racket[erf] for
@@ -185,8 +184,8 @@ manipulate @racket[(- 1.0 (erfc x))] and its surrounding expressions to avoid th
                     (flulp-error (fllog1p (- (erfc x))) log-erf-x)]
 For negative @racket[x] away from @racket[0.0], do the same with @racket[(- (erfc (- x)) 1.0)].
 
-For @racket[erf], error is no greater than 2 @tech{ulps} everywhere that has been tested, and
-is almost always no greater than 1. For @racket[erfc], observed error is no greater than 4 ulps,
+For @racket[erf], on the real line the error is no greater than 2 @tech{ulps} everywhere that has been tested, and
+is almost always no greater than 1. In the complex plane the relative error is smaller 1e-12. For @racket[erfc], observed error is no greater than 4 ulps,
 and is usually no greater than 2.
 }
 
@@ -394,10 +393,10 @@ have very dissimilar magnitudes (e.g. @racket[1e-16] and @racket[1e16]), it exhi
 @tech{catastrophic cancellation}. We are working on it.
 }
 
-@deftogether[(@defproc[(Fresnel-S [x Real]) Real]
-              @defproc[(Fresnel-C [x Real]) Real]
-              @defproc[(Fresnel-RS [x Real]) Real]
-              @defproc[(Fresnel-RC [x Real]) Real])]{
+@deftogether[(@defproc[(Fresnel-S [x Number]) Number]
+              @defproc[(Fresnel-C [x Number]) Number]
+              @defproc[(Fresnel-RS [x Number]) Number]
+              @defproc[(Fresnel-RC [x Number]) Number])]{
 Compute the @hyperlink["https://en.wikipedia.org/wiki/Fresnel_integral"]{Fresnel integrals}. Where
  @itemlist[
  @item{@racket[(Fresnel-S x)] calculates ∫sin(πt²/2) |0->x}
@@ -413,7 +412,9 @@ The first two are sometimes also referred to as the natural Fresnel integrals.
                  (Fresnel-RS 1)
                  (* (sqrt (/ pi 2)) (Fresnel-S (* (sqrt (/ 2 pi)) 1)))]
 
-Spot-checks within the region 0<=x<=150 sugest that the error is no greater than 1e-14 everywhere that has been tested, and usually is lower than 2e-15.
+Spot-checks on the real line within the region 0<=x<=150 sugest that the error is no greater than 1e-14
+everywhere that has been tested, and usually is lower than 2e-15. In the complex plane the relative error
+is smaller than 1e-12.
 }
 
 

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

@@ -14,7 +14,7 @@ IEEE Transactions on Communications, 2000, vol 48, pp 529--532
          "continued-fraction.rkt")
 
 (provide flerf flerfc*expsqr flerfc
-         erf erfc)
+         erf erfc complex-erf)
 
 ;; ===================================================================================================
 ;; erf
@@ -242,11 +242,14 @@ IEEE Transactions on Communications, 2000, vol 48, pp 529--532
 
 (: erf (case-> (Zero -> Zero)
                (Flonum -> Flonum)
-               (Real -> (U Zero Flonum))))
-(define (erf x)
-  (cond [(flonum? x)  (flerf x)]
+               (Real -> (U Zero Flonum))
+               (Number -> Number)))
+(define (erf z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
+  (cond [(flonum? x) (flerf x)]
         [(eqv? x 0)  x]
-        [else  (flerf (fl x))]))
+        [(real? x)  (flerf (fl x))]
+        [else       (complex-erf z)]))
 
 (: erfc (case-> (Zero -> One)
                 (Flonum -> Flonum)

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

@@ -20,10 +20,11 @@ would like to extend this together with erf to the complex plane
 |#
 
 (require "../../base.rkt"
-         "../../flonum.rkt")
+         "../../flonum.rkt"
+         "erf.rkt")
 
-(provide flFresnel-S Fresnel-S Fresnel-RS
-         flFresnel-C Fresnel-C Fresnel-RC)
+(provide flFresnel-S complex-Fresnel-S Fresnel-S Fresnel-RS
+         flFresnel-C complex-Fresnel-C Fresnel-C Fresnel-RC)
 
 ;------------------------------
 ;polinomials for the rational assymptotic aproximation for S & C
@@ -76,13 +77,15 @@ would like to extend this together with erf to the complex plane
 
 ;------------------------------
 ;polinomials for the rational powerseries aproximation for S
-(define sn (make-flpolyfun ( 3.18016297876567817986E11
+(define sn/d
+  (make-quotient-flpolyfun ( 3.18016297876567817986E11
                             -4.42979518059697779103E10
                              2.54890880573376359104E9
                             -6.29741486205862506537E7
                              7.08840045257738576863E5
-                            -2.99181919401019853726E3)))
-(define sd (make-flpolyfun ( 6.07366389490084639049E11
+                            -2.99181919401019853726E3
+                             0.0)
+                           ( 6.07366389490084639049E11
                              2.24411795645340920940E10
                              4.19320245898111231129E8
                              5.17343888770096400730E6
@@ -97,7 +100,7 @@ would like to extend this together with erf to the complex plane
     [(fl< x 1.5625)
      (define X2 (fl* x x))
      (define X4 (fl* X2 X2))
-     (fl* (fl* X2 x) (fl/ (sn X4)(sd X4)))]
+     (fl* (fl* X2 x) (sn/d X4))]
     [else
      (define t (fl* pi (fl* x x)))
      (define t/2 (fl/ t 2.0))
@@ -107,15 +110,23 @@ would like to extend this together with erf to the complex plane
      (fl- 0.5 (fl/ (fl+ (fl* f (flcos t/2))(fl* g (flsin t/2)))
                    (fl* pi x)))]))
 
+(define (complex-Fresnel-S [z : Number]) : Number
+  (define z* (* z (sqrt (/ pi 4))))
+  (* (/ 1+i 4)
+     (-       (complex-erf (* z* 1+i))
+        (* +i (complex-erf (* z* 1-i))))))
+
 ;------------------------------
 ;polinomials for the rational powerseries aproximation for C
-(define cn (make-flpolyfun ( 9.99999999999999998822E-1
+(define cn/d
+  (make-quotient-flpolyfun ( 9.99999999999999998822E-1
                             -2.05525900955013891793E-1
                              1.88843319396703850064E-2
                             -6.45191435683965050962E-4
                              9.50428062829859605134E-6
-                            -4.98843114573573548651E-8)))
-(define cd (make-flpolyfun ( 1.00000000000000000118E0
+                            -4.98843114573573548651E-8
+                             0.0)
+                           ( 1.00000000000000000118E0
                              4.12142090722199792936E-2
                              8.68029542941784300606E-4
                              1.22262789024179030997E-5
@@ -130,7 +141,7 @@ would like to extend this together with erf to the complex plane
     [(fl< x 1.5625)
      (define X2 (fl* x x))
      (define X4 (fl* X2 X2))
-     (fl* x (fl/ (cn X4)(cd X4)))]
+     (fl* x (cn/d X4))]
     [else
      (define t (fl* pi (fl* x x)))
      (define t/2 (fl/ t 2.0))
@@ -140,13 +151,57 @@ would like to extend this together with erf to the complex plane
      (fl+ 0.5 (fl/ (fl- (fl* f (flsin t/2))(fl* g (flcos t/2)))
                    (fl* pi x)))]))
 
+(define (complex-Fresnel-C [z : Number]) : Number
+  (define z* (* z (sqrt (/ pi 4))))
+  (* (/ 1-i 4)
+     (+       (complex-erf (* z* 1+i))
+        (* +i (complex-erf (* z* 1-i))))))
+
 ;------------------------------
-(define (Fresnel-S [x : Real]) : Real (flFresnel-S (fl x)))
-(define (Fresnel-C [x : Real]) : Real (flFresnel-C (fl x)))
-(define (Fresnel-RS [x : Real]) : Real
-  (* (flsqrt (fl/ pi 2.0)) (flFresnel-S (* (flsqrt (fl/ 2.0 pi)) (fl x)))))
-(define (Fresnel-RC [x : Real]) : Real
-  (* (flsqrt (fl/ pi 2.0)) (flFresnel-C (* (flsqrt (fl/ 2.0 pi)) (fl x)))))
+(: Fresnel-S (case-> (Zero -> Zero)
+                     (Flonum -> Flonum)
+                     (Real -> (U Zero Flonum))
+                     (Number -> Number)))
+(define (Fresnel-S z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
+  (cond
+    [(eqv? z 0) 0]
+    [(flonum? x)(flFresnel-S x)]
+    [(real? x)(flFresnel-S (fl x))]
+    [else (complex-Fresnel-S z)]))
+(: Fresnel-C (case-> (Zero -> Zero)
+                     (Flonum -> Flonum)
+                     (Real -> (U Zero Flonum))
+                     (Number -> Number)))
+(define (Fresnel-C z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
+  (cond
+    [(eqv? z 0) 0]
+    [(flonum? x)(flFresnel-C x)]
+    [(real? x)(flFresnel-C (fl x))]
+    [else (complex-Fresnel-C z)]))
+(: Fresnel-RS (case-> (Zero -> Zero)
+                      (Flonum -> Flonum)
+                      (Real -> (U Zero Flonum))
+                      (Number -> Number)))
+(define (Fresnel-RS z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
+  (cond
+    [(eqv? z 0) 0]
+    [(flonum? x)(* (flsqrt (fl/ pi 2.0))(flFresnel-S (* (flsqrt (fl/ 2.0 pi)) x)))]
+    [(real? x)  (* (flsqrt (fl/ pi 2.0))(flFresnel-S (* (flsqrt (fl/ 2.0 pi)) (fl x))))]
+    [else (* (sqrt (/ pi 2))(complex-Fresnel-S (* (sqrt (/ 2 pi)) z)))]))
+(: Fresnel-RC (case-> (Zero -> Zero)
+                      (Flonum -> Flonum)
+                      (Real -> (U Zero Flonum))
+                      (Number -> Number)))
+(define (Fresnel-RC z)
+  (define x (if (= (imag-part z) 0) (real-part z) z))
+  (cond
+    [(eqv? z 0) 0]
+    [(flonum? x)(* (flsqrt (fl/ pi 2.0))(flFresnel-C (* (flsqrt (fl/ 2.0 pi)) x)))]
+    [(real? x)  (* (flsqrt (fl/ pi 2.0))(flFresnel-C (* (flsqrt (fl/ 2.0 pi)) (fl x))))]
+    [else (* (sqrt (/ pi 2))(complex-Fresnel-C (* (sqrt (/ 2 pi)) z)))]))
 
 ;------------------------------
 (module* bfFresnel #f
@@ -209,18 +264,18 @@ 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 #f
     ;some functions to check the function.
     ;commented out, because (partly) put in function-tests from math-test
   (require math/bigfloat
-           rackunit)
+           typed/rackunit)
   (require (submod "..")
            (submod ".." bfFresnel))
 
-  (define (test a #:ε [ε 1e-15])
-    (check-= (Fresnel-S a)  (bigfloat->flonum (bfFresnel-S (bf a)))  ε (format "(Fresnel-S ~a)" a))
+  (define (test [a : Flonum] #:ε [ε : Flonum 1e-15])
+    (check-= (Fresnel-S  a) (bigfloat->flonum (bfFresnel-S  (bf a))) ε (format "(Fresnel-S ~a)" a))
     (check-= (Fresnel-RS a) (bigfloat->flonum (bfFresnel-RS (bf a))) ε (format "(Fresnel-RS ~a)" a))
-    (check-= (Fresnel-C a)  (bigfloat->flonum (bfFresnel-C (bf a)))  ε (format "(Fresnel-C ~a)" a))
+    (check-= (Fresnel-C  a) (bigfloat->flonum (bfFresnel-C  (bf a))) ε (format "(Fresnel-C ~a)" a))
     (check-= (Fresnel-RC a) (bigfloat->flonum (bfFresnel-RC (bf a))) ε (format "(Fresnel-RC ~a)" a)))
 
   (test 0.1)
@@ -242,6 +297,31 @@ would like to extend this together with erf to the complex plane
   (check-equal? (Fresnel-S -5)(-(Fresnel-S 5)))
   (check-equal? (Fresnel-C -5)(-(Fresnel-C 5)))
 
+  (check-= (magnitude
+            (/ (complex-Fresnel-S 1+i)
+               -2.0618882191948404680807165366857086008159083237378680520+2.0618882191948404680807165366857086008159083237378680520i))
+           1 1e-12)
+  (check-= (magnitude
+            (/ (complex-Fresnel-S 5+0.2i)
+               0.47365635370953447150430003290670950910437504109567910708+0.73462461062246762668741695076291749315532498862606992695i))
+           1 1e-12)
+  (check-= (magnitude
+            (/ (complex-Fresnel-S -8-25i)
+               4.33491319138289340482459027250885195089165476877024e270+1.38537242126661103439768955547584123133806956947632e270i))
+           1 1e-12)
+  (check-= (magnitude
+            (/ (complex-Fresnel-C 1+i)
+               2.55579377810243902463452238835219584215662360420358429635+2.55579377810243902463452238835219584215662360420358429635i))
+           1 1e-12)
+  (check-= (magnitude
+            (/ (complex-Fresnel-C 5+0.2i)
+               1.237351377588089955209810162536250644901272375752411527+0.02602758966318992966794130566838646516498797340046208293i))
+           1 1e-12)
+  (check-= (magnitude
+            (/ (complex-Fresnel-C -8-25i)
+               1.38537242126661103439768955547584123133806956947632e270-4.33491319138289340482459027250885195089165476877024e270i))
+           1 1e-12)
+
   #;(let ()
     (local-require plot)
     (define (mk)(* 100 (random)))

math-test/math/tests/function-tests.rkt

@@ -121,3 +121,21 @@
 (check-= (Fresnel-C 20) 0.4999873349723443881870062136976602164476 (ε))
 (check-= (Fresnel-C 50) 0.4999991894307279679558101639817919070024 (ε))
 
+  (let ([z 1+i]
+        [ans -2.0618882191948404680807165366857086008159083237378680520+2.0618882191948404680807165366857086008159083237378680520i])
+    (check-= (/ (Fresnel-S z) ans) 1 1e-12))
+  (let ([z 5+0.2i]
+        [ans 0.47365635370953447150430003290670950910437504109567910708+0.73462461062246762668741695076291749315532498862606992695i])
+    (check-= (/ (Fresnel-S z) ans) 1 1e-12))
+  (let ([z -8-25i]
+        [ans 4.33491319138289340482459027250885195089165476877024e270+1.38537242126661103439768955547584123133806956947632e270i])
+    (check-= (/ (Fresnel-S z) ans) 1 1e-12))
+  (let ([z 1+i]
+        [ans 2.55579377810243902463452238835219584215662360420358429635+2.55579377810243902463452238835219584215662360420358429635i])
+    (check-= (/ (Fresnel-C z) ans) 1 1e-12))
+  (let ([z 5+0.2i]
+        [ans 1.237351377588089955209810162536250644901272375752411527+0.02602758966318992966794130566838646516498797340046208293i])
+    (check-= (/ (Fresnel-C z) ans) 1 1e-12))
+  (let ([z -8-25i]
+        [ans 1.38537242126661103439768955547584123133806956947632e270-4.33491319138289340482459027250885195089165476877024e270i])
+    (check-= (/ (Fresnel-C z) ans) 1 1e-12))