diff --git a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-exp.rkt b/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-exp.rkt index 48976bd3bb..ab1eda5a1c 100644 --- a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-exp.rkt +++ b/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-exp.rkt @@ -35,8 +35,14 @@ (fl* x (fl/ (expm1-poly-numer x) (expm1-poly-denom x))))])) ;; Integer arguments for flexp2 -(: flexp2s (Vectorof Nonnegative-Flonum)) -(define flexp2s (build-vector (- 1024 -1074) (λ: ([n : Index]) (fl (expt 2 (- n 1074)))))) +(: flexp2s (U #f (Vectorof Nonnegative-Flonum))) +(define flexp2s #f) + +(define (build-flexp2s) + (: new-flexp2s (Vectorof Nonnegative-Flonum)) + (define new-flexp2s (build-vector (- 1024 -1074) (λ: ([n : Index]) (fl (expt 2 (- n 1074)))))) + (set! flexp2s new-flexp2s) + new-flexp2s) (begin-encourage-inline @@ -89,10 +95,14 @@ [else (fl* (flexp x) (flexp 1.0))])])) (: flexp2 (Flonum -> Nonnegative-Flonum)) + ;; Computes 2^x (define (flexp2 x) (cond [(fl<= x -1075.0) 0.0] [(fl>= x 1024.0) +inf.0] - [(fl= x (flround x)) (unsafe-vector-ref flexp2s (fl->exact-integer (fl+ x 1074.0)))] + [(fl= x (flround x)) + (let* ([flexp2s flexp2s] + [flexp2s (if flexp2s flexp2s (build-flexp2s))]) + (unsafe-vector-ref flexp2s (fl->exact-integer (fl+ x 1074.0))))] [else (flexpt 2.0 x)])) ) ; begin-encourage-inline diff --git a/pkgs/math-pkgs/math-lib/math/private/functions/incomplete-gamma/gamma-temme.rkt b/pkgs/math-pkgs/math-lib/math/private/functions/incomplete-gamma/gamma-temme.rkt index fe07d48637..63288268f4 100644 --- a/pkgs/math-pkgs/math-lib/math/private/functions/incomplete-gamma/gamma-temme.rkt +++ b/pkgs/math-pkgs/math-lib/math/private/functions/incomplete-gamma/gamma-temme.rkt @@ -15,36 +15,43 @@ (define temme-iters (make-parameter 32)) (define num-fs 100) -(define fs - (let () - (define: start-fs : (Vectorof Real) (vector 1 -1/3 1/12 -2/135)) - (define: fs : (Vectorof Real) (make-vector num-fs 0)) - (vector-copy! fs 0 start-fs) - ;; DP algorithm to compute f coefficients - (for ([m (in-range 4 num-fs)]) - (vector-set! - fs m - (* (- (/ (+ m 1) (+ m 2))) - (+ (* (/ (- m 1) (* 3 m)) (vector-ref fs (- m 1))) - (for/fold: ([sum : Real 0]) ([j (in-range 3 m)]) - (+ sum (/ (* (vector-ref fs (- j 1)) (vector-ref fs (+ m 1 (- j)))) - (+ m 2 (- j))))))))) - (vector->flvector fs))) +(define: fs : (U #f FlVector) #f) + +(: build-fs (-> FlVector)) +(define (build-fs) + (define: start-fs : (Vectorof Real) (vector 1 -1/3 1/12 -2/135)) + (define: vec-fs : (Vectorof Real) (make-vector num-fs 0)) + (vector-copy! vec-fs 0 start-fs) + ;; DP algorithm to compute f coefficients + (for ([m (in-range 4 num-fs)]) + (vector-set! + vec-fs m + (* (- (/ (+ m 1) (+ m 2))) + (+ (* (/ (- m 1) (* 3 m)) (vector-ref vec-fs (- m 1))) + (for/fold: ([sum : Real 0]) ([j (in-range 3 m)]) + (+ sum (/ (* (vector-ref vec-fs (- j 1)) + (vector-ref vec-fs (+ m 1 (- j)))) + (+ m 2 (- j))))))))) + (define new-fs (vector->flvector vec-fs)) + (set! fs new-fs) + new-fs) (: R-sum (Float Flonum -> Flonum)) (define (R-sum k n) (define num-fs (temme-iters)) - ;; This originally filled a vector of bs, because imperative programmers don't know how to do - ;; anything else besides bang an array full of values (sheesh) - (define-values (sum b2 b1) - (for/fold: ([sum : Flonum 0.0] - [b2 : Flonum (unsafe-flvector-ref fs (- num-fs 1))] - [b1 : Flonum (unsafe-flvector-ref fs (- num-fs 2))] - ) ([m (in-range (- num-fs 3) 0 -1)]) - (define c (unsafe-flvector-ref fs m)) - (define b0 (fl+ c (fl/ (fl* (fl+ (fl m) 1.0) b2) k))) - (values (fl+ (fl* n sum) b0) b1 b0))) - sum) + (let* ([fs fs] + [fs (if fs fs (build-fs))]) + ;; This originally filled a vector of bs, because imperative programmers don't know how to do + ;; anything else besides bang an array full of values (sheesh) + (define-values (sum b2 b1) + (for/fold: ([sum : Flonum 0.0] + [b2 : Flonum (unsafe-flvector-ref fs (- num-fs 1))] + [b1 : Flonum (unsafe-flvector-ref fs (- num-fs 2))] + ) ([m (in-range (- num-fs 3) 0 -1)]) + (define c (unsafe-flvector-ref fs m)) + (define b0 (fl+ c (fl/ (fl* (fl+ (fl m) 1.0) b2) k))) + (values (fl+ (fl* n sum) b0) b1 b0))) + sum)) (: R-log (Float Float -> (Values Float Float))) ;; Log-space version of `R' above diff --git a/pkgs/math-pkgs/math-lib/math/private/number-theory/small-primes.rkt b/pkgs/math-pkgs/math-lib/math/private/number-theory/small-primes.rkt index b95b5f82ce..8dc0caefa7 100644 --- a/pkgs/math-pkgs/math-lib/math/private/number-theory/small-primes.rkt +++ b/pkgs/math-pkgs/math-lib/math/private/number-theory/small-primes.rkt @@ -69,6 +69,8 @@ (clear-bit! y) (loop (cdr ds))])))) +(define: sieve-done? : Boolean #f) + (define: (sieve) : Void (define x 1) (let loop ([ds deltas]) @@ -85,9 +87,11 @@ (mark-composites x)) (loop (cdr ds))])))) -(sieve) - (define: (small-prime? [x : Integer]) : Boolean + (unless sieve-done? + (sieve) + (set! sieve-done? #t)) + (or (= x 2) (= x 3) (= x 5) (and (mod60->bit (modulo x 60)) (bit x))))