diff --git a/math-doc/math/scribblings/math-special-functions.scrbl b/math-doc/math/scribblings/math-special-functions.scrbl
index cd997d1..ef1a385 100644
--- a/math-doc/math/scribblings/math-special-functions.scrbl
+++ b/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 two groups: @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[
@@ -28,13 +28,13 @@ in terms of their flonum counterparts, but are different in two crucial ways:
@racket[exn:fail:contract] instead of returning @racket[+nan.0].}
]
-Currently, @racketmodname[math/special-functions] does not export any functions that accept
-or return complex numbers. Mathematically, some of them could return complex numbers given
-real numbers, such @racket[hurwitz-zeta] when given a negative second argument. In
+Only functions in @secref{complex-real-functions} support accepting or returning complex arguments.
+Mathematically, some of the other functions could return complex numbers given
+real numbers, such as @racket[hurwitz-zeta] when given a negative second argument. In
these cases, they raise an @racket[exn:fail:contract] (for an exact argument) or return
@racket[+nan.0] (for an inexact argument).
-Most real functions have more than one type, but they are documented as having only
+Most real and complex functions have more than one type, but they are documented as having only
one. The documented type is the most general type, which is used to generate a contract for
uses in untyped code. Use @racket[:print-type] to see all of a function's types.
@@ -56,11 +56,11 @@ 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 "real-functions"]{Real Functions}
+@section[#:tag "complex-real-functions"]{Complex and Real Functions}
-@defproc[(gamma [x Real]) (U Positive-Integer Flonum)]{
+@defproc[(gamma [x Number]) Number]{
Computes the @hyperlink["http://en.wikipedia.org/wiki/Gamma_function"]{gamma function},
-a generalization of the factorial function to the entire real line, except nonpositive integers.
+a generalization of the factorial function to the entire complex plane, except nonpositive exact integers.
When @racket[x] is an exact integer, @racket[(gamma x)] is exact.
@examples[#:eval untyped-eval
@@ -76,16 +76,17 @@ When @racket[x] is an exact integer, @racket[(gamma x)] is exact.
(gamma -1.0)
(gamma 0.0)
(gamma -0.0)
+ (gamma 1+i)
(gamma 172.0)
(eval:alts
(bf (gamma 172))
(eval:result @racketresultfont{(bf "1.241018070217667823424840524103103992618e309")}))]
-Error is no more than 10 @tech{ulps} everywhere that has been tested, and is usually no more than 4
-ulps.
+On the real line the error is no more than 10 @tech{ulps} everywhere that has been tested, and is usually no more than 4
+ulps. In the rest of the complex plane the relative error is smaller than 1e-13 except in the very close neighbourhood of negative integers, where the error can increase signifiantly.
}
-@defproc[(log-gamma [x Real]) (U Zero Flonum)]{
+@defproc[(log-gamma [x Number]) Number]{
Like @racket[(log (abs (gamma x)))], but more accurate and without unnecessary overflow.
The only exact cases are @racket[(log-gamma 1) = 0] and @racket[(log-gamma 2) = 0].
@@ -99,14 +100,15 @@ The only exact cases are @racket[(log-gamma 1) = 0] and @racket[(log-gamma 2) =
(log-gamma -1)
(log-gamma -1.0)
(log-gamma 0.0)
+ (log-gamma 1+i)
(log (abs (gamma 172.0)))
(log-gamma 172.0)]
-Error is no more than 11 @tech{ulps} everywhere that has been tested, and is usually no more than 2
+On the real line error is no more than 11 @tech{ulps} everywhere that has been tested, and is usually no more than 2
ulps. Error reaches its maximum near negative roots.
}
-@defproc[(psi0 [x Real]) Flonum]{
+@defproc[(psi0 [x Number]) Number]{
Computes the @hyperlink["http://en.wikipedia.org/wiki/Digamma_function"]{digamma function},
the logarithmic derivative of the gamma function.
@@ -114,10 +116,11 @@ the logarithmic derivative of the gamma function.
(plot (function psi0 -2.5 4.5) #:y-min -5 #:y-max 5)
(psi0 0)
(psi0 1)
+ (psi0 1+i)
(- gamma.0)]
-Except near negative roots, maximum observed error is 2 @tech{ulps}, but is usually no more
-than 1.
+On the real line, except near negative roots, maximum observed error is 2 @tech{ulps}, but is usually no more
+than 1. In the complex plane the relative error is around 1e-13.
Near negative roots, which occur singly between each pair of negative integers, @racket[psi0]
exhibits @tech{catastrophic cancellation} from using the reflection formula, meaning that
@@ -149,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]
@@ -162,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
@@ -180,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 than 1e-12. For @racket[erfc], observed error is no greater than 4 ulps,
and is usually no greater than 2.
}
@@ -389,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}
@@ -408,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.
}
diff --git a/math-lib/math/private/bigfloat/bigfloat-fresnel.rkt b/math-lib/math/private/bigfloat/bigfloat-fresnel.rkt
new file mode 100644
index 0000000..2e6ed78
--- /dev/null
+++ b/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)))
diff --git a/math-lib/math/private/functions/erf.rkt b/math-lib/math/private/functions/erf.rkt
index 5ff597a..1c78e28 100644
--- a/math-lib/math/private/functions/erf.rkt
+++ b/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
@@ -198,13 +198,59 @@ 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 z² (* z+ z+))
+ (define nmax 193)
+ (define y (* 2 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²)) (* sqrtpi z+))))))
+ (if (= (real-part z) 0)
+ (- f 1)
+ f)]))
+
+;; ===================================================================================================
+
(: 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)
diff --git a/math-lib/math/private/functions/fresnel.rkt b/math-lib/math/private/functions/fresnel.rkt
index 743daaf..aaa4464 100644
--- a/math-lib/math/private/functions/fresnel.rkt
+++ b/math-lib/math/private/functions/fresnel.rkt
@@ -12,75 +12,76 @@ 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
|#
(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
-(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
-(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
@@ -95,25 +96,33 @@ 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))
(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)))]))
+(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
@@ -128,97 +137,81 @@ 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))
(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)))]))
-;------------------------------
-(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)))))
+(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))))))
;------------------------------
-(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)))
- )
+(: 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* 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)
+ (require typed/rackunit)
(require (submod "..")
- (submod ".." bfFresnel))
+ "../../bigfloat.rkt"
+ "../bigfloat/bigfloat-fresnel.rkt")
- (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)
@@ -240,6 +233,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)))
diff --git a/math-lib/math/private/functions/gamma.rkt b/math-lib/math/private/functions/gamma.rkt
index 864d259..2cb73c3 100644
--- a/math-lib/math/private/functions/gamma.rkt
+++ b/math-lib/math/private/functions/gamma.rkt
@@ -272,13 +272,65 @@ Approximations:
[(and (x . fl> . -4.5) (x . fl< . 4.5)) (flgamma-taylor x)]
[else (flgamma-lanczos x)])))]))
+(define √2π 2.5066282746310005024157652848110)
+(: complex-gamma (Float-Complex -> Float-Complex))
+(define (complex-gamma z+1)
+ (cond
+ [(or (= z+1 2)(= z+1 1)) 1.+0.i]
+ [else
+ (define neg? (< (real-part z+1)0))
+ (define z (- (if neg? (- z+1) z+1) 1))
+ (define zh (+ z 0.5))
+ (define zgh (+ zh lanczos-complex-g))
+ (define zp (expt zgh (* 0.5 zh)))
+ (define zzz (* zp (exp (- zgh)) zp))
+
+ (define m
+ (+ (car lanczos-complex-c)
+ (for/fold ([s : Float-Complex 0.0+0.0i])
+ ([a (in-list (cdr lanczos-complex-c))]
+ [i (in-naturals 1)])
+ (+ s (/ a (+ z i))))))
+
+ ;(println m)(println zp)(println zgh)(println zzz)
+ (define (mm [z : Float-Complex])(or (infinite? (real-part z))(infinite? (imag-part z))))
+ (cond
+ [(or (mm zh)(mm zgh)) +nan.0+nan.0i]
+ [(or (mm zp)(mm zzz))
+ (cond
+ [neg? 0.0+0.0i]
+ [else
+ (define ans (* (* √2π m)
+ (make-polar 1 (+ (* (magnitude zgh) (sin (angle zgh)))
+ (* 2 (angle zgh) (magnitude zh) (sin (+ (* .5 pi) (angle zh)))))) ))
+ (define α (angle ans))
+ (cond
+ [(< 0.0 α (* 0.5 pi)) +inf.0+inf.0i]
+ [(< (* 0.5 pi) α pi) -inf.0+inf.0i]
+ [(< (- pi) α (* -.5 pi)) -inf.0-inf.0i]
+ [(< (* -.5 pi) α 0.0) +inf.0-inf.0i]
+ [(= 0.0 α) +inf.0+0.0i]
+ [(= (* 0.5 pi) α) 0.0+inf.0i]
+ [(or (= pi α)(= (- pi) α)) -inf.0+0.0i]
+ [(= (* -.5 pi) α) 0.0-inf.0i]
+ [else (error "complex-gamma: angle outside of range [-π , π]")])])]
+ ;[(mm m) (error "todo")]
+ [else
+ (define ans (* (* √2π m) zzz))
+ (if neg?
+ (/ (- pi) ans z+1 (sin (* pi z+1)))
+ ans)])]))
+
(: gamma (case-> (One -> One)
(Integer -> Positive-Integer)
(Float -> Float)
- (Real -> (U Positive-Integer Flonum))))
-(define (gamma x)
+ (Real -> (U Positive-Integer Flonum))
+ (Number -> Number)))
+(define (gamma z)
+ (define x (if (= (imag-part z) 0) (real-part z) z))
(cond [(double-flonum? x) (flgamma x)]
[(exact-integer? x)
(cond [(x . > . 0) (factorial (- x 1))]
[else (raise-argument-error 'gamma "Real, not Zero or Negative-Integer" x)])]
- [else (flgamma (fl x))]))
+ [(real? x) (flgamma (fl x))]
+ [else (complex-gamma (number->float-complex z))]))
diff --git a/math-lib/math/private/functions/lanczos.rkt b/math-lib/math/private/functions/lanczos.rkt
index b822185..fc79189 100644
--- a/math-lib/math/private/functions/lanczos.rkt
+++ b/math-lib/math/private/functions/lanczos.rkt
@@ -2,7 +2,7 @@
(require "../../flonum.rkt")
-(provide lanczos-sum lanczos-g)
+(provide (all-defined-out))
;; Lanczos polynomial for N=13 G=6.024680040776729583740234375
;; Max experimental error (with arbitary precision arithmetic) is 1.196214e-17
@@ -36,3 +36,25 @@
1.0)))
(define lanczos-g 6.024680040776729583740234375)
+
+
+;; Lanczos series approximation for the complex plane
+;; ref Paul Godfrey's MATLAB implementations | https://my.fit.edu/~gabdo/paulbio.html
+;; ->relative error < 1e-13
+(define lanczos-complex-c
+ (list 0.99999999999999709182
+ 57.156235665862923517
+ -59.597960355475491248
+ 14.136097974741747174
+ -0.49191381609762019978
+ 0.33994649984811888699e-4
+ 0.46523628927048575665e-4
+ -0.98374475304879564677e-4
+ 0.15808870322491248884e-3
+ -0.21026444172410488319e-3
+ 0.21743961811521264320e-3
+ -0.16431810653676389022e-3
+ 0.84418223983852743293e-4
+ -0.26190838401581408670e-4
+ 0.36899182659531622704e-5))
+(define lanczos-complex-g 607/128)
diff --git a/math-lib/math/private/functions/log-gamma.rkt b/math-lib/math/private/functions/log-gamma.rkt
index 37b2d8e..7a1f31c 100644
--- a/math-lib/math/private/functions/log-gamma.rkt
+++ b/math-lib/math/private/functions/log-gamma.rkt
@@ -4,7 +4,8 @@
"../../flonum.rkt"
"../../base.rkt"
"log-gamma-zeros.rkt"
- "gamma.rkt")
+ "gamma.rkt"
+ "lanczos.rkt")
(provide fllog-gamma log-gamma)
@@ -398,10 +399,37 @@
[(x . fl< . 4.5) (fllog-gamma-small-positive x)]
[else (fllog-gamma-large-positive x)])))]))
+(define log-√2π 0.9189385332046727417803297)
+(define log-π 1.14472988584940017414342735)
+(: complex-log-gamma (Number -> Number))
+(define (complex-log-gamma z)
+ (define neg? (< (real-part z)0))
+ (define z* (if neg? (- z) z))
+ (cond
+ [(or (= z* 1)(= z* 2))0]
+ [else
+ (define z- (- z* 0.5))
+ (define zg (+ z- lanczos-complex-g))
+
+ (define ans
+ (+ (- zg) (* z- (log zg))
+ log-√2π (log (+ (car lanczos-complex-c)
+ (for/sum : Number ([a (in-list (cdr lanczos-complex-c))]
+ [i (in-naturals 0)])
+ (/ a (+ z* i)))))))
+ (cond
+ [neg?
+ (define lpi (+ log-π (* (if (< (imag-part z) 0) +i -i) pi)))
+ (- lpi (log z) ans (log (sin (* pi z))))]
+ [else ans])]))
+
+
(: log-gamma (case-> (One -> Zero)
(Flonum -> Flonum)
- (Real -> (U Zero Flonum))))
-(define (log-gamma x)
+ (Real -> (U Zero Flonum))
+ (Number -> Number)))
+(define (log-gamma z)
+ (define x (if (= (imag-part z) 0) (real-part z) z))
(cond [(flonum? x) (fllog-gamma x)]
[(single-flonum? x) (fllog-gamma (fl x))]
[(integer? x)
@@ -409,4 +437,5 @@
(raise-argument-error 'log-gamma "Real, not Zero or Negative-Integer" x)]
[(eqv? x 1) 0]
[else (fllog-factorial (fl (- x 1)))])]
- [else (fllog-gamma (fl x))]))
+ [(real? x) (fllog-gamma (fl x))]
+ [else (complex-log-gamma z)]))
diff --git a/math-lib/math/private/functions/psi.rkt b/math-lib/math/private/functions/psi.rkt
index d57d2c5..8ccdc9c 100644
--- a/math-lib/math/private/functions/psi.rkt
+++ b/math-lib/math/private/functions/psi.rkt
@@ -6,7 +6,8 @@
"../number-theory/bernoulli.rkt"
"gamma.rkt"
"hurwitz-zeta.rkt"
- "tan-diff.rkt")
+ "tan-diff.rkt"
+ "lanczos.rkt")
(provide flpsi0 flpsi psi0 psi)
@@ -208,6 +209,29 @@
[(x . = . +inf.0) +inf.0]
[else +nan.0]))
+(: complex-psi0 (Number -> Number))
+(define (complex-psi0 z*)
+ (define neg? (< (real-part z*) 0.5))
+ (define z (if neg? (- 1 z*) z*))
+
+ (define-values (n d)
+ (for/fold : (Values Number Number)
+ ([n : Number 0][d : Number 0])
+ ([c (in-list (reverse (cdr lanczos-complex-c)))]
+ [k (in-range (length lanczos-complex-c) -1 -1)])
+ (define dz (/ 1 (+ z k -2)))
+ (define dd (* c dz))
+ (values (- n (* dd dz)) (+ d dd))))
+
+ (define gg (+ z lanczos-complex-g -0.5))
+
+ (define ans (+ (log gg) (- (/ n (+ d (car lanczos-complex-c)))
+ (/ lanczos-complex-g gg))))
+
+ (if neg?
+ (- ans (/ pi (tan (* pi z*))))
+ ans))
+
(define pi.128 267257146016241686964920093290467695825/85070591730234615865843651857942052864)
(: flexppi (Flonum -> Flonum))
@@ -232,13 +256,16 @@
(fl- (fl* sgn (flpsi m (fl- 1.0 x))) t)]
[else +nan.0]))
-(: psi0 (Real -> Flonum))
-(define (psi0 x)
+(: psi0 (case-> (Real -> Flonum)
+ (Number -> Number)))
+(define (psi0 z)
+ (define x (if (= (imag-part z) 0) (real-part z) z))
(cond [(flonum? x) (flpsi0 x)]
[(single-flonum? x) (flpsi0 (fl x))]
[(and (integer? x) (x . <= . 0))
(raise-argument-error 'psi0 "Real, not Zero or Negative-Integer" x)]
- [else (flpsi0 (fl x))]))
+ [(real? x) (flpsi0 (fl x))]
+ [else (complex-psi0 z)]))
(: psi (Integer Real -> Flonum))
(define (psi m x)
diff --git a/math-lib/math/special-functions.rkt b/math-lib/math/special-functions.rkt
index c620527..3ab9a43 100644
--- a/math-lib/math/special-functions.rkt
+++ b/math-lib/math/special-functions.rkt
@@ -9,19 +9,23 @@
(except-in "private/functions/lambert.rkt" lambert)
(except-in "private/functions/zeta.rkt" eta zeta)
(except-in "private/functions/hurwitz-zeta.rkt" hurwitz-zeta)
- "private/functions/psi.rkt"
+ (except-in "private/functions/psi.rkt" psi0)
"private/functions/incomplete-gamma.rkt"
"private/functions/incomplete-beta.rkt"
"private/functions/stirling-error.rkt"
- "private/functions/fresnel.rkt")
+ (except-in "private/functions/fresnel.rkt" Fresnel-S Fresnel-RS Fresnel-C Fresnel-RC))
(require/untyped-contract
"private/functions/gamma.rkt"
- [gamma (Real -> Real)])
+ [gamma (Number -> Number)])
(require/untyped-contract
"private/functions/log-gamma.rkt"
- [log-gamma (Real -> Real)])
+ [log-gamma (Number -> Number)])
+
+(require/untyped-contract
+ "private/functions/psi.rkt"
+ [psi0 (Number -> Number)])
(require/untyped-contract
"private/functions/beta.rkt"
@@ -30,8 +34,15 @@
(require/untyped-contract
"private/functions/erf.rkt"
- [erf (Real -> Real)]
+ [erf (Number -> Number)]
[erfc (Real -> Real)])
+(require/untyped-contract
+ "private/functions/fresnel.rkt"
+ [Fresnel-S (Number -> Number)]
+ [Fresnel-RS (Number -> Number)]
+ [Fresnel-C (Number -> Number)]
+ [Fresnel-RC (Number -> Number)])
+
(require/untyped-contract
"private/functions/lambert.rkt"
@@ -61,9 +72,11 @@
"private/functions/fresnel.rkt")
gamma
log-gamma
+ psi0
beta log-beta
erf erfc
lambert
eta zeta
hurwitz-zeta
- )
+ Fresnel-S Fresnel-RS
+ Fresnel-C Fresnel-RC)
diff --git a/math-test/math/tests/function-tests.rkt b/math-test/math/tests/function-tests.rkt
index 4a1fdf7..0d9d026 100644
--- a/math-test/math/tests/function-tests.rkt
+++ b/math-test/math/tests/function-tests.rkt
@@ -46,7 +46,31 @@
(check-exn exn:fail:contract? (λ () (gamma 0)))
(check-equal? (gamma 4) 6)
(check-equal? (gamma 4.0) (flgamma 4.0))
-
+;checks calculated at https://wolframalpha.com/input/?i=gamma(z) / log-gamma(z) / digamma(z)
+(let ([z 1.5+0.2i]
+ [ans 0.869862133805271779271655350309163858030873370542544143+0.00730170831118402614972641679459165642341675697339983422i])
+ (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 1.5+4i]
+ [ans -0.018686297879039387689174933326088586248898218845351961+0.0026241318932335593901663964582546948097937382268327333i])
+ (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 196+285i]
+ [ans 5.72038483529989178989490686567823916571020691258673e291-2.00024207656771597752271221062051820963168450311374e291i])
+ (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 28.45-0.05i]
+ [ans 4.788665109060112126493799529515601713890011469557001e28-8.048791086601393564384851295197107978228121888100461e27i])
+ (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z 0.5-0.00333i]
+ [ans 1.772367471310459385811039514905244407455502656029765525+0.01158858570803772503302975596977297632272041615128277137i])
+ (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z -0.5-0.00333i]
+ [ans -3.544732073864574326606517758100934969124507449000061+0.0004307441958626149491398963294062742489283873077748653i])
+ (check-= (/ (gamma z) ans) 1 1e-13))
+(let ([z -1+i]
+ [ans -0.9974967895818289935938328922330359491032588243198502190-4.228631137454774745965835886896275145906361523162647513i])
+ (check-= (/ (log-gamma z) ans) 1 1e-13))
+(let ([z -1+i]
+ [ans 0.59465032062247697727187848272191072247626297176354162323+2.5766740474685811741340507947500004904456562664038166655i])
+ (check-= (/ (psi0 z) ans) 1 1e-13))
;; ---------------------------------------------------------------------------------------------------
;; hypot
@@ -97,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))
\ No newline at end of file