From d0bde91a127a45f9d39b7b1d3122a99547b96a64 Mon Sep 17 00:00:00 2001 From: Bert De Ketelaere Date: Fri, 11 Oct 2019 14:28:07 +0200 Subject: [PATCH 1/8] Complex special functions Algorithms for calculating gamma, log-gamma and digamma in the complex field --- .../scribblings/math-special-functions.scrbl | 33 ++++++++++------- math-lib/math/private/functions/gamma.rkt | 31 ++++++++++++++-- math-lib/math/private/functions/lanczos.rkt | 24 +++++++++++- math-lib/math/private/functions/log-gamma.rkt | 37 +++++++++++++++++-- math-lib/math/private/functions/psi.rkt | 35 ++++++++++++++++-- math-lib/math/special-functions.rkt | 11 ++++-- math-test/math/tests/function-tests.rkt | 26 ++++++++++++- 7 files changed, 167 insertions(+), 30 deletions(-) diff --git a/math-doc/math/scribblings/math-special-functions.scrbl b/math-doc/math/scribblings/math-special-functions.scrbl index cd997d1..3eb64fd 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 three groups: @secref{complex-functions}, @secref{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 +Only functions in @secref{complex-functions} support accepting or returning complex arguments. +Mathematically, some of the other functions could return complex numbers given real numbers, such @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-functions"]{Complex 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 (@tech{ulps}=450). } -@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 @@ -129,6 +132,8 @@ 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 diff --git a/math-lib/math/private/functions/gamma.rkt b/math-lib/math/private/functions/gamma.rkt index 864d259..bf8dd8f 100644 --- a/math-lib/math/private/functions/gamma.rkt +++ b/math-lib/math/private/functions/gamma.rkt @@ -272,13 +272,38 @@ Approximations: [(and (x . fl> . -4.5) (x . fl< . 4.5)) (flgamma-taylor x)] [else (flgamma-lanczos x)])))])) +(define √2π 2.5066282746310005024157652848110) +(: complex-gamma (Number -> Number)) +(define (complex-gamma z+1) + (define neg? (< (real-part z+1)0)) + (define z (- (if neg? (- z+1) z+1) 1)) + (cond + [(or (= z 1)(= z 0))1] + [else + (define zh (+ z 0.5)) + (define zgh (+ zh lanczos-complex-g)) + (define zp (expt zgh (* 0.5 zh))) + + (define ans + (* zp (exp (- zgh)) zp + √2π (+ (car lanczos-complex-c) + (for/sum : Number ([a (in-list (cdr lanczos-complex-c))] + [i (in-naturals 1)]) + (/ a (+ z i)))))) + (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 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..61155cc 100644 --- a/math-lib/math/special-functions.rkt +++ b/math-lib/math/special-functions.rkt @@ -9,7 +9,7 @@ (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" @@ -17,11 +17,15 @@ (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" @@ -61,6 +65,7 @@ "private/functions/fresnel.rkt") gamma log-gamma + psi0 beta log-beta erf erfc lambert diff --git a/math-test/math/tests/function-tests.rkt b/math-test/math/tests/function-tests.rkt index 4a1fdf7..2e1de32 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 From 3f3cd78c118269dd945bac8864a7ba0742c1459e Mon Sep 17 00:00:00 2001 From: Bert De Ketelaere Date: Fri, 11 Oct 2019 18:03:43 +0200 Subject: [PATCH 2/8] Complex erf, start fresnel --- math-lib/math/private/functions/erf.rkt | 42 ++++++++ math-lib/math/private/functions/fresnel.rkt | 100 ++++++++++---------- 2 files changed, 93 insertions(+), 49 deletions(-) diff --git a/math-lib/math/private/functions/erf.rkt b/math-lib/math/private/functions/erf.rkt index 5ff597a..c700698 100644 --- a/math-lib/math/private/functions/erf.rkt +++ b/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)))) diff --git a/math-lib/math/private/functions/fresnel.rkt b/math-lib/math/private/functions/fresnel.rkt index 743daaf..d52a816 100644 --- a/math-lib/math/private/functions/fresnel.rkt +++ b/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 From 6d959a9aaddbafecf3f00c1ba70d662b266eddb1 Mon Sep 17 00:00:00 2001 From: Bert De Ketelaere Date: Tue, 15 Oct 2019 12:45:13 +0200 Subject: [PATCH 3/8] Complex erf and Fresnel functions --- .../scribblings/math-special-functions.scrbl | 27 ++-- math-lib/math/private/functions/erf.rkt | 13 +- math-lib/math/private/functions/fresnel.rkt | 124 ++++++++++++++---- math-test/math/tests/function-tests.rkt | 18 +++ 4 files changed, 142 insertions(+), 40 deletions(-) diff --git a/math-doc/math/scribblings/math-special-functions.scrbl b/math-doc/math/scribblings/math-special-functions.scrbl index 3eb64fd..90529e8 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 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. } diff --git a/math-lib/math/private/functions/erf.rkt b/math-lib/math/private/functions/erf.rkt index c700698..0b03b32 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 @@ -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) diff --git a/math-lib/math/private/functions/fresnel.rkt b/math-lib/math/private/functions/fresnel.rkt index d52a816..dbdb282 100644 --- a/math-lib/math/private/functions/fresnel.rkt +++ b/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))) diff --git a/math-test/math/tests/function-tests.rkt b/math-test/math/tests/function-tests.rkt index 2e1de32..0d9d026 100644 --- a/math-test/math/tests/function-tests.rkt +++ b/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)) \ No newline at end of file From fe463f53b75f3d4b27e617a4af89c9e243cb3d4a Mon Sep 17 00:00:00 2001 From: Bert De Ketelaere Date: Tue, 15 Oct 2019 12:55:25 +0200 Subject: [PATCH 4/8] separate bigfloat-fresnel into it's own file --- .../private/bigfloat/bigfloat-fresnel.rkt | 65 +++++++++++++++++ math-lib/math/private/functions/fresnel.rkt | 72 ++----------------- 2 files changed, 69 insertions(+), 68 deletions(-) create mode 100644 math-lib/math/private/bigfloat/bigfloat-fresnel.rkt 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/fresnel.rkt b/math-lib/math/private/functions/fresnel.rkt index dbdb282..aaa4464 100644 --- a/math-lib/math/private/functions/fresnel.rkt +++ b/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)) From 9363d58c45a368d5b286490a638f2fcdb00f2b09 Mon Sep 17 00:00:00 2001 From: Bert De Ketelaere Date: Tue, 15 Oct 2019 14:31:35 +0200 Subject: [PATCH 5/8] Update 'type' for untyped code --- math-lib/math/special-functions.rkt | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/math-lib/math/special-functions.rkt b/math-lib/math/special-functions.rkt index 61155cc..3ab9a43 100644 --- a/math-lib/math/special-functions.rkt +++ b/math-lib/math/special-functions.rkt @@ -13,7 +13,7 @@ "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" @@ -34,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" @@ -71,4 +78,5 @@ lambert eta zeta hurwitz-zeta - ) + Fresnel-S Fresnel-RS + Fresnel-C Fresnel-RC) From 91856a87e067b651df0105336d320c8a094e57cb Mon Sep 17 00:00:00 2001 From: Bert De Ketelaere Date: Sat, 25 Jan 2020 22:00:58 +0100 Subject: [PATCH 6/8] Update according to comments --- math-doc/math/scribblings/math-special-functions.scrbl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/math-doc/math/scribblings/math-special-functions.scrbl b/math-doc/math/scribblings/math-special-functions.scrbl index 90529e8..ef1a385 100644 --- a/math-doc/math/scribblings/math-special-functions.scrbl +++ b/math-doc/math/scribblings/math-special-functions.scrbl @@ -28,9 +28,9 @@ in terms of their flonum counterparts, but are different in two crucial ways: @racket[exn:fail:contract] instead of returning @racket[+nan.0].} ] -Only functions in @secref{complex-functions} support accepting or returning complex arguments. +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 @racket[hurwitz-zeta] when given a negative second argument. In +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). @@ -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-real-functions"]{Complex and real 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}, @@ -83,7 +83,7 @@ When @racket[x] is an exact integer, @racket[(gamma x)] is exact. (eval:result @racketresultfont{(bf "1.241018070217667823424840524103103992618e309")}))] 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 (@tech{ulps}=450). +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 Number]) Number]{ @@ -185,7 +185,7 @@ manipulate @racket[(- 1.0 (erfc x))] and its surrounding expressions to avoid th For negative @racket[x] away from @racket[0.0], do the same with @racket[(- (erfc (- x)) 1.0)]. 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, +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. } From 77c855e07a2ce63df79067434eabafd03aca2c86 Mon Sep 17 00:00:00 2001 From: bdeket Date: Sun, 5 Jul 2020 23:50:20 +0200 Subject: [PATCH 7/8] right sign for values with negative real-part --- math-lib/math/private/functions/erf.rkt | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/math-lib/math/private/functions/erf.rkt b/math-lib/math/private/functions/erf.rkt index 0b03b32..1c78e28 100644 --- a/math-lib/math/private/functions/erf.rkt +++ b/math-lib/math/private/functions/erf.rkt @@ -226,14 +226,15 @@ IEEE Transactions on Communications, 2000, vol 48, pp 529--532 [else (define neg? (< (real-part z) 0)) (define z+ (if neg? (- z) z)) + (define z² (* z+ z+)) (define nmax 193) - (define y (* 2 z z)) + (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 z)) (* sqrtpi z)))))) + (- 1 (* s (/ (exp (* -1 z²)) (* sqrtpi z+)))))) (if (= (real-part z) 0) (- f 1) f)])) From db737eb3f1cf91e498d4fb17519d4148dbda5c06 Mon Sep 17 00:00:00 2001 From: bdeket Date: Tue, 7 Jul 2020 00:53:27 +0200 Subject: [PATCH 8/8] complex-gamma: improve accuracy and sign for large (magnitude of) z --- math-lib/math/private/functions/gamma.rkt | 55 +++++++++++++++++------ 1 file changed, 41 insertions(+), 14 deletions(-) diff --git a/math-lib/math/private/functions/gamma.rkt b/math-lib/math/private/functions/gamma.rkt index bf8dd8f..2cb73c3 100644 --- a/math-lib/math/private/functions/gamma.rkt +++ b/math-lib/math/private/functions/gamma.rkt @@ -273,26 +273,53 @@ Approximations: [else (flgamma-lanczos x)])))])) (define √2π 2.5066282746310005024157652848110) -(: complex-gamma (Number -> Number)) +(: complex-gamma (Float-Complex -> Float-Complex)) (define (complex-gamma z+1) - (define neg? (< (real-part z+1)0)) - (define z (- (if neg? (- z+1) z+1) 1)) (cond - [(or (= z 1)(= z 0))1] + [(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 ans - (* zp (exp (- zgh)) zp - √2π (+ (car lanczos-complex-c) - (for/sum : Number ([a (in-list (cdr lanczos-complex-c))] - [i (in-naturals 1)]) - (/ a (+ z i)))))) - (if neg? - (- (/ pi (* ans z+1 (sin (* pi z+1))))) - ans)])) + (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) @@ -306,4 +333,4 @@ Approximations: (cond [(x . > . 0) (factorial (- x 1))] [else (raise-argument-error 'gamma "Real, not Zero or Negative-Integer" x)])] [(real? x) (flgamma (fl x))] - [else (complex-gamma z)])) + [else (complex-gamma (number->float-complex z))]))