Defer building tables; drops 600ms from (require math) on my machine

This commit is contained in:
Neil Toronto 2014-03-26 18:56:49 -06:00
parent cd2632c186
commit c9aba4de94
3 changed files with 52 additions and 31 deletions

View File

@ -35,8 +35,14 @@
(fl* x (fl/ (expm1-poly-numer x) (expm1-poly-denom x))))])) (fl* x (fl/ (expm1-poly-numer x) (expm1-poly-denom x))))]))
;; Integer arguments for flexp2 ;; Integer arguments for flexp2
(: flexp2s (Vectorof Nonnegative-Flonum)) (: flexp2s (U #f (Vectorof Nonnegative-Flonum)))
(define flexp2s (build-vector (- 1024 -1074) (λ: ([n : Index]) (fl (expt 2 (- n 1074)))))) (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 (begin-encourage-inline
@ -89,10 +95,14 @@
[else (fl* (flexp x) (flexp 1.0))])])) [else (fl* (flexp x) (flexp 1.0))])]))
(: flexp2 (Flonum -> Nonnegative-Flonum)) (: flexp2 (Flonum -> Nonnegative-Flonum))
;; Computes 2^x
(define (flexp2 x) (define (flexp2 x)
(cond [(fl<= x -1075.0) 0.0] (cond [(fl<= x -1075.0) 0.0]
[(fl>= x 1024.0) +inf.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)])) [else (flexpt 2.0 x)]))
) ; begin-encourage-inline ) ; begin-encourage-inline

View File

@ -15,36 +15,43 @@
(define temme-iters (make-parameter 32)) (define temme-iters (make-parameter 32))
(define num-fs 100) (define num-fs 100)
(define fs (define: fs : (U #f FlVector) #f)
(let ()
(define: start-fs : (Vectorof Real) (vector 1 -1/3 1/12 -2/135)) (: build-fs (-> FlVector))
(define: fs : (Vectorof Real) (make-vector num-fs 0)) (define (build-fs)
(vector-copy! fs 0 start-fs) (define: start-fs : (Vectorof Real) (vector 1 -1/3 1/12 -2/135))
;; DP algorithm to compute f coefficients (define: vec-fs : (Vectorof Real) (make-vector num-fs 0))
(for ([m (in-range 4 num-fs)]) (vector-copy! vec-fs 0 start-fs)
(vector-set! ;; DP algorithm to compute f coefficients
fs m (for ([m (in-range 4 num-fs)])
(* (- (/ (+ m 1) (+ m 2))) (vector-set!
(+ (* (/ (- m 1) (* 3 m)) (vector-ref fs (- m 1))) vec-fs m
(for/fold: ([sum : Real 0]) ([j (in-range 3 m)]) (* (- (/ (+ m 1) (+ m 2)))
(+ sum (/ (* (vector-ref fs (- j 1)) (vector-ref fs (+ m 1 (- j)))) (+ (* (/ (- m 1) (* 3 m)) (vector-ref vec-fs (- m 1)))
(+ m 2 (- j))))))))) (for/fold: ([sum : Real 0]) ([j (in-range 3 m)])
(vector->flvector fs))) (+ 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)) (: R-sum (Float Flonum -> Flonum))
(define (R-sum k n) (define (R-sum k n)
(define num-fs (temme-iters)) (define num-fs (temme-iters))
;; This originally filled a vector of bs, because imperative programmers don't know how to do (let* ([fs fs]
;; anything else besides bang an array full of values (sheesh) [fs (if fs fs (build-fs))])
(define-values (sum b2 b1) ;; This originally filled a vector of bs, because imperative programmers don't know how to do
(for/fold: ([sum : Flonum 0.0] ;; anything else besides bang an array full of values (sheesh)
[b2 : Flonum (unsafe-flvector-ref fs (- num-fs 1))] (define-values (sum b2 b1)
[b1 : Flonum (unsafe-flvector-ref fs (- num-fs 2))] (for/fold: ([sum : Flonum 0.0]
) ([m (in-range (- num-fs 3) 0 -1)]) [b2 : Flonum (unsafe-flvector-ref fs (- num-fs 1))]
(define c (unsafe-flvector-ref fs m)) [b1 : Flonum (unsafe-flvector-ref fs (- num-fs 2))]
(define b0 (fl+ c (fl/ (fl* (fl+ (fl m) 1.0) b2) k))) ) ([m (in-range (- num-fs 3) 0 -1)])
(values (fl+ (fl* n sum) b0) b1 b0))) (define c (unsafe-flvector-ref fs m))
sum) (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))) (: R-log (Float Float -> (Values Float Float)))
;; Log-space version of `R' above ;; Log-space version of `R' above

View File

@ -69,6 +69,8 @@
(clear-bit! y) (clear-bit! y)
(loop (cdr ds))])))) (loop (cdr ds))]))))
(define: sieve-done? : Boolean #f)
(define: (sieve) : Void (define: (sieve) : Void
(define x 1) (define x 1)
(let loop ([ds deltas]) (let loop ([ds deltas])
@ -85,9 +87,11 @@
(mark-composites x)) (mark-composites x))
(loop (cdr ds))])))) (loop (cdr ds))]))))
(sieve)
(define: (small-prime? [x : Integer]) : Boolean (define: (small-prime? [x : Integer]) : Boolean
(unless sieve-done?
(sieve)
(set! sieve-done? #t))
(or (= x 2) (= x 3) (= x 5) (or (= x 2) (= x 3) (= x 5)
(and (mod60->bit (modulo x 60)) (and (mod60->bit (modulo x 60))
(bit x)))) (bit x))))