From 6009eed8d23d3c54a1cb5c85f422d3a69cbd1805 Mon Sep 17 00:00:00 2001 From: Neil Toronto Date: Wed, 28 Nov 2012 15:53:38 -0700 Subject: [PATCH] Moved flvector functions into math/flonum Sped up normal distribution sampling procedure (2x for large samples) --- collects/math/flonum.rkt | 6 +- .../math/private/array/fcarray-struct.rkt | 1 - .../math/private/array/flarray-pointwise.rkt | 2 +- .../math/private/array/flarray-struct.rkt | 1 - .../private/distributions/bernoulli-dist.rkt | 1 - .../math/private/distributions/beta-dist.rkt | 1 - .../private/distributions/binomial-dist.rkt | 1 - .../private/distributions/cauchy-dist.rkt | 1 - .../distributions/exponential-dist.rkt | 1 - .../math/private/distributions/gamma-dist.rkt | 1 - .../private/distributions/geometric-dist.rkt | 1 - .../distributions/impl/binomial-random.rkt | 1 - .../distributions/impl/gamma-random.rkt | 1 - .../distributions/impl/normal-random.rkt | 70 ++++++++++++++----- .../distributions/impl/poisson-random.rkt | 1 - .../private/distributions/logistic-dist.rkt | 1 - .../private/distributions/normal-dist.rkt | 1 - .../private/distributions/triangle-dist.rkt | 1 - .../private/distributions/uniform-dist.rkt | 1 - .../{vector => flonum}/flvector-syntax.rkt | 0 .../private/{vector => flonum}/flvector.rkt | 46 +----------- collects/math/private/functions/gamma.rkt | 1 - .../incomplete-gamma/gamma-temme.rkt | 1 - collects/math/private/functions/tan-diff.rkt | 3 +- collects/math/private/vector/vector-fft.rkt | 5 +- collects/math/vector.rkt | 4 +- 26 files changed, 64 insertions(+), 90 deletions(-) rename collects/math/private/{vector => flonum}/flvector-syntax.rkt (100%) rename collects/math/private/{vector => flonum}/flvector.rkt (84%) diff --git a/collects/math/flonum.rkt b/collects/math/flonum.rkt index eeece17536..0f2e831de7 100644 --- a/collects/math/flonum.rkt +++ b/collects/math/flonum.rkt @@ -15,7 +15,8 @@ "private/flonum/flonum-syntax.rkt" "private/flonum/expansion/expansion-base.rkt" "private/flonum/expansion/expansion-exp.rkt" - "private/flonum/expansion/expansion-log.rkt") + "private/flonum/expansion/expansion-log.rkt" + "private/flonum/flvector.rkt") (provide (all-from-out "private/flonum/flonum-bits.rkt" @@ -32,7 +33,8 @@ "private/flonum/flonum-syntax.rkt" "private/flonum/expansion/expansion-base.rkt" "private/flonum/expansion/expansion-exp.rkt" - "private/flonum/expansion/expansion-log.rkt") + "private/flonum/expansion/expansion-log.rkt" + "private/flonum/flvector.rkt") lg* lg/ lgprod) (define-syntax lg* (make-rename-transformer #'fl+)) diff --git a/collects/math/private/array/fcarray-struct.rkt b/collects/math/private/array/fcarray-struct.rkt index 68b985b336..c07cedc4a3 100644 --- a/collects/math/private/array/fcarray-struct.rkt +++ b/collects/math/private/array/fcarray-struct.rkt @@ -10,7 +10,6 @@ (module defs typed/racket/base (require "../../flonum.rkt" - "../vector/flvector.rkt" "../unsafe.rkt" "array-struct.rkt" "utils.rkt" diff --git a/collects/math/private/array/flarray-pointwise.rkt b/collects/math/private/array/flarray-pointwise.rkt index 5a23d38eb4..1fb27c154a 100644 --- a/collects/math/private/array/flarray-pointwise.rkt +++ b/collects/math/private/array/flarray-pointwise.rkt @@ -2,7 +2,7 @@ (require racket/flonum (for-syntax racket/base) - "../vector/flvector.rkt" + "../../flonum.rkt" "array-struct.rkt" "array-broadcast.rkt" "array-pointwise.rkt" diff --git a/collects/math/private/array/flarray-struct.rkt b/collects/math/private/array/flarray-struct.rkt index 648b6e96d9..75d5474913 100644 --- a/collects/math/private/array/flarray-struct.rkt +++ b/collects/math/private/array/flarray-struct.rkt @@ -10,7 +10,6 @@ (require "../../flonum.rkt" "../unsafe.rkt" - "../vector/flvector.rkt" "array-struct.rkt" "utils.rkt" "for-each.rkt") diff --git a/collects/math/private/distributions/bernoulli-dist.rkt b/collects/math/private/distributions/bernoulli-dist.rkt index 11c019cbbe..e331fc19e4 100644 --- a/collects/math/private/distributions/bernoulli-dist.rkt +++ b/collects/math/private/distributions/bernoulli-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") diff --git a/collects/math/private/distributions/beta-dist.rkt b/collects/math/private/distributions/beta-dist.rkt index 0af8ff42bd..ef35926287 100644 --- a/collects/math/private/distributions/beta-dist.rkt +++ b/collects/math/private/distributions/beta-dist.rkt @@ -4,7 +4,6 @@ racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "../functions/beta.rkt" "../functions/incomplete-beta.rkt" diff --git a/collects/math/private/distributions/binomial-dist.rkt b/collects/math/private/distributions/binomial-dist.rkt index f663b308cb..a41b72493e 100644 --- a/collects/math/private/distributions/binomial-dist.rkt +++ b/collects/math/private/distributions/binomial-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "../functions/incomplete-beta.rkt" (prefix-in impl: "impl/binomial-pdf.rkt") diff --git a/collects/math/private/distributions/cauchy-dist.rkt b/collects/math/private/distributions/cauchy-dist.rkt index 82294e32a9..e410ee0696 100644 --- a/collects/math/private/distributions/cauchy-dist.rkt +++ b/collects/math/private/distributions/cauchy-dist.rkt @@ -4,7 +4,6 @@ racket/promise "../../flonum.rkt" "../../base.rkt" - "../../vector.rkt" "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") diff --git a/collects/math/private/distributions/exponential-dist.rkt b/collects/math/private/distributions/exponential-dist.rkt index 272f48e9d1..f4e9c48383 100644 --- a/collects/math/private/distributions/exponential-dist.rkt +++ b/collects/math/private/distributions/exponential-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") diff --git a/collects/math/private/distributions/gamma-dist.rkt b/collects/math/private/distributions/gamma-dist.rkt index 22a72d9c03..147953c06e 100644 --- a/collects/math/private/distributions/gamma-dist.rkt +++ b/collects/math/private/distributions/gamma-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "../functions/incomplete-gamma.rkt" "impl/gamma-pdf.rkt" diff --git a/collects/math/private/distributions/geometric-dist.rkt b/collects/math/private/distributions/geometric-dist.rkt index 4b289b29f6..573bb07704 100644 --- a/collects/math/private/distributions/geometric-dist.rkt +++ b/collects/math/private/distributions/geometric-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") diff --git a/collects/math/private/distributions/impl/binomial-random.rkt b/collects/math/private/distributions/impl/binomial-random.rkt index c19dd58c56..e64d9fac80 100644 --- a/collects/math/private/distributions/impl/binomial-random.rkt +++ b/collects/math/private/distributions/impl/binomial-random.rkt @@ -6,7 +6,6 @@ Wolfgang Hormann. The Generation of Binomial Random Variates. (require "../../../base.rkt" "../../../flonum.rkt" - "../../../vector.rkt" "../../unsafe.rkt" "normal-random.rkt") diff --git a/collects/math/private/distributions/impl/gamma-random.rkt b/collects/math/private/distributions/impl/gamma-random.rkt index ddb4a421ec..adfda9b161 100644 --- a/collects/math/private/distributions/impl/gamma-random.rkt +++ b/collects/math/private/distributions/impl/gamma-random.rkt @@ -19,7 +19,6 @@ For others: sum of Gamma and Exponential variables, Normal approximation. |# (require "../../../flonum.rkt" - "../../../vector.rkt" "../../unsafe.rkt" "normal-random.rkt") diff --git a/collects/math/private/distributions/impl/normal-random.rkt b/collects/math/private/distributions/impl/normal-random.rkt index 9fa37f48b9..f7eea0d91d 100644 --- a/collects/math/private/distributions/impl/normal-random.rkt +++ b/collects/math/private/distributions/impl/normal-random.rkt @@ -1,30 +1,68 @@ #lang typed/racket/base -(require "../../../flonum.rkt" +;; Return random samples from a normal distribution using the Box-Muller transform + +(require racket/fixnum + "../../../flonum.rkt" "../../../base.rkt" - "../../../vector.rkt") + "../../unsafe.rkt") (provide flnormal-sample) -(: box-muller-transform (Float Float -> Float)) -(define (box-muller-transform x y) - (cond [(and (fl= x 0.0) (fl= y 0.0)) 0.0] - [else (fl* (flsqrt (fl* -2.0 (fllog x))) - (flsin (fl* (fl* 2.0 pi) y)))])) +;; Leaving these in, in case we discover in the future that it's actually important for them +;; to be accurate +#| +(: flsin2pix (Flonum -> Flonum)) +;; Computes sin(2*pi*x) accurately in the range [0,1] +(define (flsin2pix x) + (let*-values ([(x s) (if (x . fl> . 0.5) (values (fl- x 0.5) -1.0) (values x 1.0))] + [(x) (if (x . fl> . 0.25) (fl- 0.5 x) x)]) + (fl* s (flsin (fl* (fl* 2.0 pi) x))))) + +(: flcos2pix (Flonum -> Flonum)) +;; Computes cos(2*pi*x) accurately in the range [0,1] +(define (flcos2pix x) + (let*-values ([(x) (if (x . fl> . 0.5) (fl- 1.0 x) x)] + [(x s) (if (x . fl> . 0.25) (values (fl- 0.5 x) 1.0) (values x -1.0))]) + (fl* s (flsin (fl* (fl* 2.0 pi) (fl- x 0.25)))))) +|# + +(: nonzero-random (-> Flonum)) +(define (nonzero-random) + (let ([u (random)]) + (if (fl= u 0.0) (nonzero-random) u))) (: flnormal-sample (Flonum Flonum Integer -> FlVector)) ;; The Box-Muller method has an bad reputation, but undeservedly: ;; 1. There's nothing unstable about the floating-point implementation of the transform ;; 2. It has good tail behavior (i.e. it can return very unlikely numbers) -;; 3. With today's RNGs, there's no need to worry about generating two random numbers -;; 4. With today's FPUs, there's no need to worry about computing `log' and `sin' (sheesh) +;; 3. With today's FPUs, there's no need to worry about computing `log' and `sin' (sheesh) ;; Points in favor: it's simple and fast (define (flnormal-sample c s n) - (cond [(n . < . 0) (raise-argument-error 'flnormal-sample "Natural" 2 c s n)] + (cond [(not (index? n)) (raise-argument-error 'flnormal-sample "Natural" 2 c s n)] [else - (build-flvector - n (λ (_) - (let loop () - (define r (box-muller-transform (random) (random))) - (if (rational? r) (fl+ c (fl* s r)) (loop)))))])) - + (define xs (make-flvector n)) + (cond + [(fx= n 0) xs] + [else + (define n-1 (fx- n 1)) + (let loop ([#{i : Nonnegative-Fixnum} 0]) + (cond [(i . fx< . n-1) + (define u1 (nonzero-random)) + (define u2 (random)) + (define t (flsqrt (fl* -2.0 (fllog u1)))) + (define z (fl* (fl* 2.0 pi) u2)) + (define x (fl* t (flcos z))) + (define y (fl* t (flsin z))) + (unsafe-flvector-set! xs i (fl+ c (fl* s x))) + (unsafe-flvector-set! xs (fx+ i 1) (fl+ c (fl* s y))) + (loop (fx+ i 2))] + [(i . fx= . n-1) + (define u1 (nonzero-random)) + (define u2 (random)) + (define x (fl* (flsqrt (fl* -2.0 (fllog u1))) + (flsin (fl* (fl* 2.0 pi) u2)))) + (unsafe-flvector-set! xs i (fl+ c (fl* s x))) + xs] + [else + xs]))])])) diff --git a/collects/math/private/distributions/impl/poisson-random.rkt b/collects/math/private/distributions/impl/poisson-random.rkt index 3541393ca2..bc2ebe317b 100644 --- a/collects/math/private/distributions/impl/poisson-random.rkt +++ b/collects/math/private/distributions/impl/poisson-random.rkt @@ -3,7 +3,6 @@ (require racket/fixnum "../../../flonum.rkt" "../../../base.rkt" - "../../../vector.rkt" "../../functions/log-gamma.rkt") (provide flpoisson-sample) diff --git a/collects/math/private/distributions/logistic-dist.rkt b/collects/math/private/distributions/logistic-dist.rkt index 20df61931c..8ad7a14f83 100644 --- a/collects/math/private/distributions/logistic-dist.rkt +++ b/collects/math/private/distributions/logistic-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") diff --git a/collects/math/private/distributions/normal-dist.rkt b/collects/math/private/distributions/normal-dist.rkt index c6ae10e1d4..919c431139 100644 --- a/collects/math/private/distributions/normal-dist.rkt +++ b/collects/math/private/distributions/normal-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "impl/normal-pdf.rkt" "impl/normal-cdf.rkt" diff --git a/collects/math/private/distributions/triangle-dist.rkt b/collects/math/private/distributions/triangle-dist.rkt index 79a2d75152..a9aa0f6f0a 100644 --- a/collects/math/private/distributions/triangle-dist.rkt +++ b/collects/math/private/distributions/triangle-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "../inline-sort.rkt" "dist-struct.rkt" diff --git a/collects/math/private/distributions/uniform-dist.rkt b/collects/math/private/distributions/uniform-dist.rkt index 5992cafece..aa797545fa 100644 --- a/collects/math/private/distributions/uniform-dist.rkt +++ b/collects/math/private/distributions/uniform-dist.rkt @@ -3,7 +3,6 @@ (require racket/performance-hint racket/promise "../../flonum.rkt" - "../../vector.rkt" "../unsafe.rkt" "dist-struct.rkt" "utils.rkt") diff --git a/collects/math/private/vector/flvector-syntax.rkt b/collects/math/private/flonum/flvector-syntax.rkt similarity index 100% rename from collects/math/private/vector/flvector-syntax.rkt rename to collects/math/private/flonum/flvector-syntax.rkt diff --git a/collects/math/private/vector/flvector.rkt b/collects/math/private/flonum/flvector.rkt similarity index 84% rename from collects/math/private/vector/flvector.rkt rename to collects/math/private/flonum/flvector.rkt index 1c89c3c7dc..c07c87be37 100644 --- a/collects/math/private/vector/flvector.rkt +++ b/collects/math/private/flonum/flvector.rkt @@ -3,7 +3,7 @@ (require racket/fixnum racket/string (for-syntax racket/base syntax/parse) - "../../flonum.rkt" + "flonum-functions.rkt" "../unsafe.rkt" "flvector-syntax.rkt") @@ -12,9 +12,6 @@ ;; Construction unsafe-flvector-copy! flvector-copy! - ;; Loops - for/flvector: - for*/flvector: ;; Conversion list->flvector flvector->list @@ -90,47 +87,6 @@ [else (unsafe-flvector-copy! dest dest-start src src-start src-end)])])) -;; =================================================================================================== -;; Loops - -(define-syntax (base-for/flvector: stx) - (syntax-parse stx - [(_ for: #:length n-expr:expr (clauses ...) body ...+) - (syntax/loc stx - (let: ([n : Integer n-expr]) - (cond [(n . > . 0) - (define xs (make-flvector n)) - (define: i : Nonnegative-Fixnum 0) - (let/ec: break : Void - (for: (clauses ...) - (unsafe-flvector-set! xs i (let () body ...)) - (set! i (unsafe-fx+ i 1)) - (when (i . unsafe-fx>= . n) (break (void))))) - xs] - [else (flvector)])))] - [(_ for: (clauses ...) body ...+) - (syntax/loc stx - (let () - (define n 4) - (define xs (make-flvector 4)) - (define i 0) - (for: (clauses ...) - (let: ([x : Float (let () body ...)]) - (cond [(unsafe-fx= i n) (define new-n (unsafe-fx* 2 n)) - (define new-xs (make-flvector new-n x)) - (unsafe-flvector-copy! new-xs 0 xs 0 n) - (set! n new-n) - (set! xs new-xs)] - [else (unsafe-flvector-set! xs i x)])) - (set! i (unsafe-fx+ i 1))) - (flvector-copy xs 0 i)))])) - -(define-syntax-rule (for/flvector: e ...) - (base-for/flvector: for: e ...)) - -(define-syntax-rule (for*/flvector: e ...) - (base-for/flvector: for*: e ...)) - ;; =================================================================================================== ;; Conversion diff --git a/collects/math/private/functions/gamma.rkt b/collects/math/private/functions/gamma.rkt index 1d8b0111fb..864d259a98 100644 --- a/collects/math/private/functions/gamma.rkt +++ b/collects/math/private/functions/gamma.rkt @@ -25,7 +25,6 @@ Approximations: "../../flonum.rkt" "../../base.rkt" "../number-theory/factorial.rkt" - "../vector/flvector.rkt" "lanczos.rkt") (provide gamma flgamma) diff --git a/collects/math/private/functions/incomplete-gamma/gamma-temme.rkt b/collects/math/private/functions/incomplete-gamma/gamma-temme.rkt index a1327b174a..fe07d48637 100644 --- a/collects/math/private/functions/incomplete-gamma/gamma-temme.rkt +++ b/collects/math/private/functions/incomplete-gamma/gamma-temme.rkt @@ -2,7 +2,6 @@ (require "../../../flonum.rkt" "../../../base.rkt" - "../../../vector.rkt" "../../unsafe.rkt" "../../distributions/impl/normal-cdf.rkt" "../stirling-error.rkt" diff --git a/collects/math/private/functions/tan-diff.rkt b/collects/math/private/functions/tan-diff.rkt index 9d73707115..51c9fc23cb 100644 --- a/collects/math/private/functions/tan-diff.rkt +++ b/collects/math/private/functions/tan-diff.rkt @@ -2,8 +2,7 @@ (require racket/fixnum "../../flonum.rkt" - "../vector/vector.rkt" - "../vector/flvector.rkt") + "../vector/vector.rkt") (provide fltan-diff/y flcot-diff/y) diff --git a/collects/math/private/vector/vector-fft.rkt b/collects/math/private/vector/vector-fft.rkt index 00ee5a009f..9232c702da 100644 --- a/collects/math/private/vector/vector-fft.rkt +++ b/collects/math/private/vector/vector-fft.rkt @@ -5,10 +5,9 @@ racket/list racket/future "../../base.rkt" + "../../flonum.rkt" "../../parameters.rkt" - "../unsafe.rkt" - "flvector.rkt" - "vector.rkt") + "../unsafe.rkt") (provide vector-fft flvector-fft! vector-inverse-fft flvector-inverse-fft!) diff --git a/collects/math/vector.rkt b/collects/math/vector.rkt index 950db94fde..00daa8e72a 100644 --- a/collects/math/vector.rkt +++ b/collects/math/vector.rkt @@ -1,10 +1,8 @@ #lang typed/racket/base -(require "private/vector/flvector.rkt" - "private/vector/vector.rkt" +(require "private/vector/vector.rkt" "private/vector/vector-fft.rkt") (provide (all-from-out - "private/vector/flvector.rkt" "private/vector/vector.rkt" "private/vector/vector-fft.rkt"))