From 5a43f2c6bc3b6cb846f3ee0dbaafeefd76dc688a Mon Sep 17 00:00:00 2001 From: Neil Toronto Date: Thu, 29 Nov 2012 13:34:09 -0700 Subject: [PATCH] Finished array documentation! Cleaned up other docs in preparation for alpha-testing announcement Created `math/utils' module for stuff that doesn't go anywhere else (e.g. FFT scaling convention, max-math-threads parameters) Reduced the number of macros that expand to applications of `array-map' Added `flvector-sum', defined `flsum' in terms of it Reduced the number of pointwise `flvector', `flarray' and `fcarray' operations Reworked `inline-build-flvector' and `inline-flvector-map' to be faster and expand to less code in both typed and untyped Racket Redefined conversions like `list->flvector' in terms of for loops (can do it now that TR has working `for/flvector:', etc.) --- collects/math/flonum.rkt | 2 - collects/math/main.rkt | 4 +- collects/math/private/array/array-fft.rkt | 5 +- .../math/private/array/array-parallel.rkt | 2 +- .../math/private/array/array-pointwise.rkt | 56 +-- .../math/private/array/fcarray-pointwise.rkt | 91 ++-- .../math/private/array/flarray-pointwise.rkt | 111 +---- collects/math/private/base/base-functions.rkt | 2 +- collects/math/private/flonum/flonum-log.rkt | 4 +- .../private/flonum/flonum-more-functions.rkt | 3 +- collects/math/private/flonum/flonum-sum.rkt | 38 -- .../math/private/flonum/flvector-syntax.rkt | 168 ++++---- collects/math/private/flonum/flvector.rkt | 256 ++++++----- collects/math/{ => private}/parameters.rkt | 0 collects/math/private/utils.rkt | 14 +- collects/math/private/vector/vector-fft.rkt | 2 +- collects/math/scribblings/math-array.scrbl | 405 ++++++++++-------- collects/math/scribblings/math-base.scrbl | 2 +- collects/math/scribblings/math-bigfloat.scrbl | 2 + .../math/scribblings/math-distributions.scrbl | 34 +- collects/math/scribblings/math-flonum.scrbl | 96 ++++- .../scribblings/math-special-functions.scrbl | 2 +- .../math/scribblings/math-statistics.scrbl | 7 +- collects/math/scribblings/math-utils.scrbl | 45 ++ collects/math/scribblings/math.scrbl | 5 +- collects/math/utils.rkt | 6 + 26 files changed, 699 insertions(+), 663 deletions(-) delete mode 100644 collects/math/private/flonum/flonum-sum.rkt rename collects/math/{ => private}/parameters.rkt (100%) create mode 100644 collects/math/scribblings/math-utils.scrbl create mode 100644 collects/math/utils.rkt diff --git a/collects/math/flonum.rkt b/collects/math/flonum.rkt index 0f2e831de7..80d0db1d75 100644 --- a/collects/math/flonum.rkt +++ b/collects/math/flonum.rkt @@ -7,7 +7,6 @@ "private/flonum/flonum-search.rkt" "private/flonum/flonum-exp.rkt" "private/flonum/flonum-log.rkt" - "private/flonum/flonum-sum.rkt" "private/flonum/flonum-more-functions.rkt" "private/flonum/flonum-factorial.rkt" "private/flonum/flonum-log1pmx.rkt" @@ -25,7 +24,6 @@ "private/flonum/flonum-search.rkt" "private/flonum/flonum-exp.rkt" "private/flonum/flonum-log.rkt" - "private/flonum/flonum-sum.rkt" "private/flonum/flonum-more-functions.rkt" "private/flonum/flonum-factorial.rkt" "private/flonum/flonum-log1pmx.rkt" diff --git a/collects/math/main.rkt b/collects/math/main.rkt index 76b5d84d5b..f038ae03dd 100644 --- a/collects/math/main.rkt +++ b/collects/math/main.rkt @@ -10,7 +10,7 @@ "vector.rkt" "array.rkt" ;"matrix.rkt" - "parameters.rkt") + "utils.rkt") (provide (all-from-out "base.rkt" @@ -23,4 +23,4 @@ "vector.rkt" "array.rkt" ;"matrix.rkt" - "parameters.rkt")) + "utils.rkt")) diff --git a/collects/math/private/array/array-fft.rkt b/collects/math/private/array/array-fft.rkt index 88b7d8144b..711958c98f 100644 --- a/collects/math/private/array/array-fft.rkt +++ b/collects/math/private/array/array-fft.rkt @@ -1,8 +1,8 @@ #lang typed/racket/base -(require "../../parameters.rkt" - "../../base.rkt" +(require "../../base.rkt" "../../flonum.rkt" + "../parameters.rkt" "../unsafe.rkt" "../vector/vector-fft.rkt" "fcarray-struct.rkt" @@ -13,7 +13,6 @@ (provide array-axis-fft array-fft - fcarray-fft array-axis-inverse-fft array-inverse-fft) diff --git a/collects/math/private/array/array-parallel.rkt b/collects/math/private/array/array-parallel.rkt index f24157250b..fffbf53e56 100644 --- a/collects/math/private/array/array-parallel.rkt +++ b/collects/math/private/array/array-parallel.rkt @@ -3,7 +3,7 @@ (require racket/future racket/list "../unsafe.rkt" - "../../parameters.rkt" + "../parameters.rkt" "array-struct.rkt" "mutable-array.rkt" "utils.rkt") diff --git a/collects/math/private/array/array-pointwise.rkt b/collects/math/private/array/array-pointwise.rkt index 01d55a3091..f624f50890 100644 --- a/collects/math/private/array/array-pointwise.rkt +++ b/collects/math/private/array/array-pointwise.rkt @@ -28,38 +28,23 @@ (define-syntax-rule (array-scale arr x) (inline-array-map (λ (y) (* x y)) arr)) +(define-array-op1 array-sqr sqr) +(define-array-op1 array-sqrt sqrt) (define-array-op1 array-abs abs) -(define-array-op1 array-round round) -(define-array-op1 array-floor floor) -(define-array-op1 array-ceiling ceiling) -(define-array-op1 array-truncate truncate) -(define-array-op1 array-conjugate conjugate) (define-array-op1 array-magnitude magnitude) (define-array-op1 array-angle angle) -(define-array-op1 array-sqrt sqrt) -(define-array-op1 array-log log) -(define-array-op1 array-sqr sqr) -(define-array-op1 array-exp exp) -(define-array-op1 array-sin sin) -(define-array-op1 array-cos cos) -(define-array-op1 array-tan tan) -(define-array-op1 array-asin asin) -(define-array-op1 array-acos acos) -(define-array-op1 array-atan atan) -(define-array-op1 array-inexact->exact inexact->exact) -(define-array-op1 array-exact->inexact exact->inexact) -(define-array-op1 array-fl real->double-flonum) -(define-array-op1 array-fc number->float-complex) +(define-array-op1 array-conjugate conjugate) (define-array-op1 array-real-part real-part) (define-array-op1 array-imag-part imag-part) + (define-array-op2 array-make-rectangular make-rectangular) +(define-array-op2 array-make-polar make-polar) (define-array-op array+ +) (define-array-op array* *) (define-array-op1+ array- -) (define-array-op1+ array/ /) -(define-array-op2 array-expt expt) (define-array-op1+ array-min min) (define-array-op1+ array-max max) @@ -81,29 +66,20 @@ array-map ;; Lifted operators array-scale - array-abs - array-round - array-floor - array-ceiling - array-truncate array-sqr array-sqrt - array-conjugate + array-abs array-magnitude array-angle - array-log - array-exp - array-sin - array-cos - array-tan - array-asin - array-acos - array-atan + array-conjugate + array-real-part + array-imag-part + array-make-rectangular + array-make-polar array+ array- array* array/ - array-expt array-min array-max array= @@ -114,12 +90,4 @@ array-not array-and array-or - array-if - ;; Number conversions - array-inexact->exact - array-exact->inexact - array-fl - array-fc - array-real-part - array-imag-part - array-make-rectangular) + array-if) diff --git a/collects/math/private/array/fcarray-pointwise.rkt b/collects/math/private/array/fcarray-pointwise.rkt index 602098b5fa..76cd0fed31 100644 --- a/collects/math/private/array/fcarray-pointwise.rkt +++ b/collects/math/private/array/fcarray-pointwise.rkt @@ -20,25 +20,17 @@ fcarray-scale fcarray-sqr fcarray-sqrt - fcarray-conjugate fcarray-magnitude fcarray-angle - fcarray-log - fcarray-exp - fcarray-sin - fcarray-cos - fcarray-tan - fcarray-asin - fcarray-acos - fcarray-atan fcarray+ fcarray* fcarray- fcarray/ - fcarray-expt + fcarray-conjugate fcarray-real-part fcarray-imag-part fcarray-make-rectangular + fcarray-make-polar ) ;; =================================================================================================== @@ -143,9 +135,6 @@ ;; =================================================================================================== ;; Pointwise operations -(define-syntax-rule (lift1 f) - (λ (arr) (inline-fcarray-map f arr))) - (define-syntax-rule (lift1->fl f) (λ (arr) (define ds (array-shape arr)) @@ -162,73 +151,63 @@ [else (unsafe-flarray ds new-xs)])))) +(define-syntax-rule (lift1 f) + (λ: ([arr : FCArray]) (inline-fcarray-map f arr))) + (define-syntax-rule (lift2 f) - (λ (arr1 arr2) (inline-fcarray-map f arr1 arr2))) + (λ: ([arr1 : FCArray] [arr2 : FCArray]) (inline-fcarray-map f arr1 arr2))) (: fcarray-scale (FCArray (U Float Float-Complex) -> FCArray)) - -(: fcarray-sqr (FCArray -> FCArray)) -(: fcarray-sqrt (FCArray -> FCArray)) -(: fcarray-conjugate (FCArray -> FCArray)) -(: fcarray-magnitude (FCArray -> FlArray)) -(: fcarray-angle (FCArray -> FlArray)) -(: fcarray-log (FCArray -> FCArray)) -(: fcarray-exp (FCArray -> FCArray)) -(: fcarray-sin (FCArray -> FCArray)) -(: fcarray-cos (FCArray -> FCArray)) -(: fcarray-tan (FCArray -> FCArray)) -(: fcarray-asin (FCArray -> FCArray)) -(: fcarray-acos (FCArray -> FCArray)) -(: fcarray-atan (FCArray -> FCArray)) - -(: fcarray+ (FCArray FCArray -> FCArray)) -(: fcarray* (FCArray FCArray -> FCArray)) -(: fcarray- (case-> (FCArray -> FCArray) - (FCArray FCArray -> FCArray))) -(: fcarray/ (case-> (FCArray -> FCArray) - (FCArray FCArray -> FCArray))) -(: fcarray-expt (FCArray FCArray -> FCArray)) - -(: fcarray-real-part (FCArray -> FlArray)) -(: fcarray-imag-part (FCArray -> FlArray)) -(: fcarray-make-rectangular (FlArray FlArray -> FCArray)) - (define (fcarray-scale arr y) (if (flonum? y) (inline-fcarray-map (λ (z) (* z y)) arr) (inline-fcarray-map (λ (z) (* z y)) arr))) -(define fcarray-sqr (lift1 (λ (x) (* x x)))) -(define fcarray-sqrt (lift1 sqrt)) -(define fcarray-conjugate (lift1 conjugate)) -(define fcarray-magnitude (lift1->fl magnitude)) -(define fcarray-angle (lift1->fl angle)) -(define fcarray-log (lift1 log)) -(define fcarray-exp (lift1 exp)) -(define fcarray-sin (lift1 sin)) -(define fcarray-cos (lift1 cos)) -(define fcarray-tan (lift1 tan)) -(define fcarray-asin (lift1 asin)) -(define fcarray-acos (lift1 acos)) -(define fcarray-atan (lift1 atan)) +(: fcarray-sqr (FCArray -> FCArray)) +(define fcarray-sqr (lift1 (λ: ([z : Float-Complex]) (* z z)))) +(: fcarray-sqrt (FCArray -> FCArray)) +(define fcarray-sqrt (lift1 sqrt)) + +(: fcarray-magnitude (FCArray -> FlArray)) +(define fcarray-magnitude (lift1->fl magnitude)) + +(: fcarray-angle (FCArray -> FlArray)) +(define fcarray-angle (lift1->fl angle)) + +(: fcarray-conjugate (FCArray -> FCArray)) +(define fcarray-conjugate (lift1 conjugate)) + +(: fcarray+ (FCArray FCArray -> FCArray)) (define fcarray+ (lift2 +)) + +(: fcarray* (FCArray FCArray -> FCArray)) (define fcarray* (lift2 *)) +(: fcarray- (case-> (FCArray -> FCArray) + (FCArray FCArray -> FCArray))) (define fcarray- (case-lambda [(arr) (inline-fcarray-map (λ (z) (- 0.0 z)) arr)] [(arr1 arr2) (inline-fcarray-map - arr1 arr2)])) +(: fcarray/ (case-> (FCArray -> FCArray) + (FCArray FCArray -> FCArray))) (define fcarray/ (case-lambda [(arr) (inline-fcarray-map (λ (z) (/ 1.0 z)) arr)] [(arr1 arr2) (inline-fcarray-map / arr1 arr2)])) -(define fcarray-expt (lift2 expt)) - +(: fcarray-real-part (FCArray -> FlArray)) (define fcarray-real-part (lift1->fl real-part)) + +(: fcarray-imag-part (FCArray -> FlArray)) (define fcarray-imag-part (lift1->fl imag-part)) +(: fcarray-make-rectangular (FlArray FlArray -> FCArray)) (define (fcarray-make-rectangular arr1 arr2) (array->fcarray (array-make-rectangular arr1 arr2))) + +(: fcarray-make-polar (FlArray FlArray -> FCArray)) +(define (fcarray-make-polar arr1 arr2) + (array->fcarray (array-make-polar arr1 arr2))) diff --git a/collects/math/private/array/flarray-pointwise.rkt b/collects/math/private/array/flarray-pointwise.rkt index 1fb27c154a..e0e0945dcc 100644 --- a/collects/math/private/array/flarray-pointwise.rkt +++ b/collects/math/private/array/flarray-pointwise.rkt @@ -16,33 +16,15 @@ flarray-map ;; Pointwise operations flarray-scale - flarray-round - flarray-floor - flarray-ceiling - flarray-truncate - flarray-abs flarray-sqr flarray-sqrt - flarray-log - flarray-exp - flarray-sin - flarray-cos - flarray-tan - flarray-asin - flarray-acos - flarray-atan + flarray-abs flarray+ flarray* flarray- flarray/ - flarray-expt flarray-min - flarray-max - flarray= - flarray< - flarray<= - flarray> - flarray>=) + flarray-max) ;; =================================================================================================== ;; Mapping @@ -53,7 +35,7 @@ [(_ f arr-expr) (syntax/loc stx (let: ([arr : FlArray arr-expr]) - (unsafe-flarray (array-shape arr) (flvector-map f (flarray-data arr)))))] + (unsafe-flarray (array-shape arr) (inline-flvector-map f (flarray-data arr)))))] [(_ f arr-expr arr-exprs ...) (with-syntax ([(arrs ...) (generate-temporaries #'(arr-exprs ...))] [(dss ...) (generate-temporaries #'(arr-exprs ...))] @@ -65,7 +47,7 @@ (define dss (array-shape arrs)) ... (cond [(and (equal? ds dss) ...) (unsafe-flarray - ds (flvector-map f (flarray-data arr) (flarray-data arrs) ...))] + ds (inline-flvector-map f (flarray-data arr) (flarray-data arrs) ...))] [else (define new-ds (array-shape-broadcast (list ds dss ...))) (define proc (unsafe-array-proc (array-broadcast arr new-ds))) @@ -115,91 +97,42 @@ (cond [(equal? ds1 ds2) (unsafe-flarray ds1 (f (flarray-data arr1) (flarray-data arr2)))] [else (array->flarray (array-f arr1 arr2))]))) -(define-syntax (lift2 stx) - (syntax-case stx () - [(_ f) (syntax/loc stx (λ (arr1 arr2) (inline-flarray-map f arr1 arr2)))])) - -(define-syntax-rule (lift-flvector-compare flvector-comp array-comp) - (λ (arr1 arr2) - (define ds1 (array-shape arr1)) - (define ds2 (array-shape arr2)) - (cond [(equal? ds1 ds2) - (unsafe-mutable-array ds1 (flvector-comp (flarray-data arr1) (flarray-data arr2)))] - [else (array-comp arr1 arr2)]))) - (: flarray-scale (FlArray Float -> FlArray)) - -(: flarray-round (FlArray -> FlArray)) -(: flarray-floor (FlArray -> FlArray)) -(: flarray-ceiling (FlArray -> FlArray)) -(: flarray-truncate (FlArray -> FlArray)) -(: flarray-abs (FlArray -> FlArray)) -(: flarray-sqr (FlArray -> FlArray)) -(: flarray-sqrt (FlArray -> FlArray)) -(: flarray-log (FlArray -> FlArray)) -(: flarray-exp (FlArray -> FlArray)) -(: flarray-sin (FlArray -> FlArray)) -(: flarray-cos (FlArray -> FlArray)) -(: flarray-tan (FlArray -> FlArray)) -(: flarray-asin (FlArray -> FlArray)) -(: flarray-acos (FlArray -> FlArray)) -(: flarray-atan (FlArray -> FlArray)) - -(: flarray+ (FlArray FlArray -> FlArray)) -(: flarray* (FlArray FlArray -> FlArray)) -(: flarray- (case-> (FlArray -> FlArray) - (FlArray FlArray -> FlArray))) -(: flarray/ (case-> (FlArray -> FlArray) - (FlArray FlArray -> FlArray))) -(: flarray-expt (FlArray FlArray -> FlArray)) -(: flarray-min (FlArray FlArray -> FlArray)) -(: flarray-max (FlArray FlArray -> FlArray)) - -(: flarray= (FlArray FlArray -> (Array Boolean))) -(: flarray< (FlArray FlArray -> (Array Boolean))) -(: flarray<= (FlArray FlArray -> (Array Boolean))) -(: flarray> (FlArray FlArray -> (Array Boolean))) -(: flarray>= (FlArray FlArray -> (Array Boolean))) - (define (flarray-scale arr y) (define-syntax-rule (fun xs) (flvector-scale xs y)) ((lift-flvector1 fun) arr)) -(define flarray-round (lift-flvector1 flvector-round)) -(define flarray-floor (lift-flvector1 flvector-floor)) -(define flarray-ceiling (lift-flvector1 flvector-ceiling)) -(define flarray-truncate (lift-flvector1 flvector-truncate)) -(define flarray-abs (lift-flvector1 flvector-abs)) -(define flarray-sqr (lift-flvector1 flvector-sqr)) -(define flarray-sqrt (lift-flvector1 flvector-sqrt)) -(define flarray-log (lift-flvector1 flvector-log)) -(define flarray-exp (lift-flvector1 flvector-exp)) -(define flarray-sin (lift-flvector1 flvector-sin)) -(define flarray-cos (lift-flvector1 flvector-cos)) -(define flarray-tan (lift-flvector1 flvector-tan)) -(define flarray-asin (lift-flvector1 flvector-asin)) -(define flarray-acos (lift-flvector1 flvector-acos)) -(define flarray-atan (lift-flvector1 flvector-atan)) +(: flarray-sqr (FlArray -> FlArray)) +(define flarray-sqr (lift-flvector1 flvector-sqr)) +(: flarray-sqrt (FlArray -> FlArray)) +(define flarray-sqrt (lift-flvector1 flvector-sqrt)) + +(: flarray-abs (FlArray -> FlArray)) +(define flarray-abs (lift-flvector1 flvector-abs)) + +(: flarray+ (FlArray FlArray -> FlArray)) (define flarray+ (lift-flvector2 flvector+ array+)) + +(: flarray* (FlArray FlArray -> FlArray)) (define flarray* (lift-flvector2 flvector* array*)) +(: flarray- (case-> (FlArray -> FlArray) + (FlArray FlArray -> FlArray))) (define flarray- (case-lambda [(arr) ((lift-flvector1 flvector-) arr)] [(arr1 arr2) ((lift-flvector2 flvector- array-) arr1 arr2)])) +(: flarray/ (case-> (FlArray -> FlArray) + (FlArray FlArray -> FlArray))) (define flarray/ (case-lambda [(arr) ((lift-flvector1 flvector/) arr)] [(arr1 arr2) ((lift-flvector2 flvector/ array/) arr1 arr2)])) -(define flarray-expt (lift2 flexpt)) +(: flarray-min (FlArray FlArray -> FlArray)) (define flarray-min (lift-flvector2 flvector-min array-min)) -(define flarray-max (lift-flvector2 flvector-max array-max)) -(define flarray= (lift-flvector-compare flvector= array=)) -(define flarray< (lift-flvector-compare flvector< array<)) -(define flarray<= (lift-flvector-compare flvector<= array<=)) -(define flarray> (lift-flvector-compare flvector> array>)) -(define flarray>= (lift-flvector-compare flvector>= array>=)) +(: flarray-max (FlArray FlArray -> FlArray)) +(define flarray-max (lift-flvector2 flvector-max array-max)) diff --git a/collects/math/private/base/base-functions.rkt b/collects/math/private/base/base-functions.rkt index 9535dd4b75..55d53ed2f6 100644 --- a/collects/math/private/base/base-functions.rkt +++ b/collects/math/private/base/base-functions.rkt @@ -1,7 +1,7 @@ #lang typed/racket/base (require racket/flonum - "../flonum/flonum-sum.rkt" + "../flonum/flvector.rkt" "../flonum/flonum-functions.rkt" "../flonum/flonum-more-functions.rkt") diff --git a/collects/math/private/flonum/flonum-log.rkt b/collects/math/private/flonum/flonum-log.rkt index 7ef2abac7d..1ca75ccdc9 100644 --- a/collects/math/private/flonum/flonum-log.rkt +++ b/collects/math/private/flonum/flonum-log.rkt @@ -5,8 +5,8 @@ "flonum-functions.rkt" "flonum-constants.rkt" "flonum-exp.rkt" - "flonum-sum.rkt" - "flonum-syntax.rkt") + "flonum-syntax.rkt" + "flvector.rkt") (provide fllog1p fllog+ lg1+ lg+ lg1- lg- lgsum diff --git a/collects/math/private/flonum/flonum-more-functions.rkt b/collects/math/private/flonum/flonum-more-functions.rkt index f0b159b317..1a4e133fda 100644 --- a/collects/math/private/flonum/flonum-more-functions.rkt +++ b/collects/math/private/flonum/flonum-more-functions.rkt @@ -5,7 +5,8 @@ "flonum-constants.rkt" "flonum-exp.rkt" "flonum-log.rkt" - "flonum-syntax.rkt") + "flonum-syntax.rkt" + "flvector.rkt") (provide flsqrt1pm1 flsinh flcosh fltanh diff --git a/collects/math/private/flonum/flonum-sum.rkt b/collects/math/private/flonum/flonum-sum.rkt deleted file mode 100644 index 2cadd06107..0000000000 --- a/collects/math/private/flonum/flonum-sum.rkt +++ /dev/null @@ -1,38 +0,0 @@ -#lang typed/racket/base - -;; Computes sums without incurring rounding error more than once -;; See (flsum '(1.0 1e-16 1e-16)) vs. (+ 1.0 1e-16 1e-16) - -#| -Algorithm adapted from: - -J R Shewchuk. Adaptive Precision Floating-Point Arithmetic and Fast Geometric Predicates. -Discrete & Computational Geometry, 1996, vol 18, pp 305--363. -|# - -(require "../unsafe.rkt" - "flonum-functions.rkt") - -(provide flsum) - -(: flsum ((Listof Flonum) -> Flonum)) -(define (flsum xs) - (define ys (make-flvector (length xs))) - (define num-ys - (for/fold: ([num-ys : Nonnegative-Fixnum 0]) ([x (in-list xs)]) - (define-values (new-x i) - (for/fold: ([x : Flonum x] [i : Nonnegative-Fixnum 0]) ([p (in-range num-ys)]) - (define y (unsafe-flvector-ref ys p)) - (let-values ([(x y) (if ((flabs x) . fl< . (flabs y)) (values y x) (values x y))]) - (define z (fl+ x y)) - (define-values (hi lo) - (cond [(rational? z) (values z (fl- y (fl- z x)))] - [else (values x y)])) - (cond [(fl= lo 0.0) (values hi i)] - [else (unsafe-flvector-set! ys i lo) - (values hi (unsafe-fx+ i 1))])))) - (unsafe-flvector-set! ys i new-x) - (unsafe-fx+ i 1))) - - (for/fold: ([sum : Flonum 0.0]) ([p (in-range num-ys)]) - (fl+ sum (unsafe-flvector-ref ys p)))) diff --git a/collects/math/private/flonum/flvector-syntax.rkt b/collects/math/private/flonum/flvector-syntax.rkt index ecd3b6e4c1..f7d2b1a544 100644 --- a/collects/math/private/flonum/flvector-syntax.rkt +++ b/collects/math/private/flonum/flvector-syntax.rkt @@ -1,79 +1,119 @@ #lang racket/base -(module syntax-defs* racket/base +(provide inline-build-flvector + build-flvector + inline-flvector-map + flvector-map) + +(module syntax-defs racket/base - (require (for-syntax racket/base) + (require (for-syntax racket/base + syntax/parse + typed/untyped-utils) typed/racket/base + racket/unsafe/ops racket/flonum racket/fixnum - racket/unsafe/ops "../syntax-utils.rkt" - "../exception.rkt") + "../exception.rkt" + "../utils.rkt") (provide (all-defined-out)) - (define-syntax-rule (ensure-flvector/length name xs-expr n) - (let*: ([xs : FlVector (ensure-flvector name xs-expr)]) - (if (fx= (unsafe-flvector-length xs) n) xs (raise-length-error name "FlVector" xs n)))) + (define-syntax (unsafe-flvector-fill! stx) + (syntax-parse stx + [(_ xs:id n:id f-expr:expr) + (syntax/loc stx + (let: loop : FlVector ([i : Nonnegative-Fixnum 0]) + (if (i . unsafe-fx< . n) + (begin (unsafe-flvector-set! xs i (f-expr i)) + (loop (unsafe-fx+ i 1))) + xs)))])) - (define-syntax-rule (unsafe-build-flvector n f) - (let ([xs (make-flvector n)]) - (let: loop : FlVector ([i : Nonnegative-Fixnum 0]) - (if (i . unsafe-fx< . n) - (begin (unsafe-flvector-set! xs i (f i)) - (loop (unsafe-fx+ i 1))) - xs)))) + (define-syntax (inline-build-flvector stx) + (syntax-case stx () + [(_ n-expr f-expr) + (cond + [(syntax-local-typed-context?) + (syntax/loc stx + (let*: ([xs : FlVector (make-flvector (ann n-expr Integer))] + [n : Index (flvector-length xs)]) + (unsafe-flvector-fill! xs n (ann f-expr (Index -> Flonum)))))] + [else + (syntax/loc stx + (let* ([xs (make-flvector n-expr)] + [n (flvector-length xs)]) + (define-syntax-rule (new-f i) + (let ([x (f-expr i)]) + (if (flonum? x) x (raise-result-error 'build-flvector "Flonum" x)))) + (unsafe-flvector-fill! xs n new-f)))])])) - (define-syntax-rule (inline-build-flvector* n-expr f-expr) - (let ([n (ensure-index 'build-flvector n-expr)] - [f (ensure-procedure 'build-flvector f-expr (Index -> Flonum))]) - (define-syntax-rule (new-f i) - (let ([x (f i)]) - (if (flonum? x) x (raise-result-error 'build-flvector "Flonum" x)))) - (unsafe-build-flvector n new-f))) - - (define-syntax (inline-flvector-map* stx) + (define-syntax (inline-flvector-map stx) (syntax-case stx () [(_ f-expr xs-expr) - (syntax/loc stx - (let ([f (ensure-procedure 'flvector-map f-expr (Flonum -> Flonum))] - [xs (ensure-flvector 'flvector-map xs-expr)]) - (define n (unsafe-flvector-length xs)) - (define-syntax-rule (new-f i) - (let ([y (f (unsafe-flvector-ref xs i))]) - (if (flonum? y) y (raise-result-error 'flvector-map "Flonum" y)))) - (unsafe-build-flvector n new-f)))] + (cond + [(syntax-local-typed-context?) + (syntax/loc stx + (let*: ([xs : FlVector xs-expr] + [n : Index (flvector-length xs)]) + (define-syntax-rule (new-f i) + ((ann f-expr (Flonum -> Flonum)) (unsafe-flvector-ref xs i))) + (define ys (make-flvector n)) + (unsafe-flvector-fill! ys n new-f)))] + [else + (syntax/loc stx + (let* ([xs (ensure-flvector 'flvector-map xs-expr)] + [n (flvector-length xs)]) + (define-syntax-rule (new-f i) + (let ([y (f-expr (unsafe-flvector-ref xs i))]) + (if (flonum? y) y (raise-result-error 'flvector-map "Flonum" y)))) + (define ys (make-flvector n)) + (unsafe-flvector-fill! ys n new-f)))])] [(_ f-expr xs-expr xss-expr ...) (with-syntax ([(f) (generate-temporaries #'(f-expr))] [(xs xss ...) (generate-temporaries #'(xs-expr xss-expr ...))] [(n ns ...) (generate-temporaries #'(xs-expr xss-expr ...))] [(Flonums ...) (build-list (length (syntax->list #'(xs-expr xss-expr ...))) (λ _ #'Flonum))]) - (syntax/loc stx - (let* ([f (ensure-procedure 'flvector-map f-expr (Flonums ... -> Flonum))] - [xs (ensure-flvector 'flvector-map xs-expr)] - [n (unsafe-flvector-length xs)] - [xss (ensure-flvector/length 'flvector-map xss-expr n)] ...) - (define-syntax-rule (new-f i) - (let ([y (f (unsafe-flvector-ref xs i) (unsafe-flvector-ref xss i) ...)]) - (cond [(flonum? y) y] - [else (raise-result-error 'flvector-map "Flonum" y)]))) - (unsafe-build-flvector n new-f))))])) + (cond + [(syntax-local-typed-context?) + (syntax/loc stx + (let*: ([xs : FlVector xs-expr] + [n : Index (flvector-length xs)] + [xss : FlVector xss-expr] ...) + (check-flvector-lengths! 'flvector-map n xss ...) + (define-syntax-rule (new-f i) + ((ann f-expr (Flonums ... -> Flonum)) (unsafe-flvector-ref xs i) + (unsafe-flvector-ref xss i) ...)) + (define ys (make-flvector n)) + (unsafe-flvector-fill! ys n new-f)))] + [else + (syntax/loc stx + (let* ([xs (ensure-flvector 'flvector-map xs-expr)] + [n (unsafe-flvector-length xs)] + [xss (ensure-flvector 'flvector-map xss-expr)] ...) + (check-flvector-lengths! 'flvector-map n xss ...) + (define-syntax-rule (new-f i) + (let ([y (f-expr (unsafe-flvector-ref xs i) (unsafe-flvector-ref xss i) ...)]) + (cond [(flonum? y) y] + [else (raise-result-error 'flvector-map "Flonum" y)]))) + (define ys (make-flvector n)) + (unsafe-flvector-fill! ys n new-f)))]))])) - ) + ) ; module (module defs typed/racket/base (require racket/flonum racket/fixnum racket/unsafe/ops - (submod ".." syntax-defs*) + (submod ".." syntax-defs) + "../utils.rkt" "../exception.rkt") (provide (all-defined-out)) (: build-flvector (Integer (Index -> Flonum) -> FlVector)) - (define (build-flvector n f) - (inline-build-flvector* n f)) + (define (build-flvector n f) (inline-build-flvector n f)) (: flvector-map (case-> ((Flonum -> Flonum) FlVector -> FlVector) ((Flonum Flonum Flonum * -> Flonum) FlVector FlVector FlVector * @@ -81,46 +121,20 @@ (define flvector-map (case-lambda: [([f : (Flonum -> Flonum)] [xs : FlVector]) - (inline-flvector-map* f xs)] + (inline-flvector-map f xs)] [([f : (Flonum Flonum -> Flonum)] [xs : FlVector] [ys : FlVector]) - (inline-flvector-map* f xs ys)] + (inline-flvector-map f xs ys)] [([f : (Flonum Flonum Flonum * -> Flonum)] [xs : FlVector] [ys : FlVector] . [yss : FlVector *]) (define n (flvector-length xs)) - (for: ([ys (in-list (cons ys yss))]) - (unless (fx= n (flvector-length ys)) (raise-length-error 'flvector-map "FlVector" ys n))) - (inline-build-flvector* + (apply check-flvector-lengths! 'flvector-map n ys yss) + (inline-build-flvector n (λ: ([i : Index]) (apply f (unsafe-flvector-ref xs i) (unsafe-flvector-ref ys i) (map (λ: ([ys : FlVector]) (unsafe-flvector-ref ys i)) yss))))])) - ) + ) ; module -(module syntax-defs racket/base - (require (for-syntax racket/base) - (submod ".." syntax-defs*) - (submod ".." defs)) - - (provide (all-defined-out)) - - (define-syntax (inline-build-flvector stx) - (syntax-case stx () - [(_ n-expr f-expr) (syntax/loc stx (inline-build-flvector* n-expr f-expr))] - [(_ . args) (syntax/loc stx (build-flvector . args))] - [_ (syntax/loc stx build-flvector)])) - - (define-syntax (inline-flvector-map stx) - (syntax-case stx () - [(_ f-expr xs-expr xss-expr ...) - (syntax/loc stx (inline-flvector-map* f-expr xs-expr xss-expr ...))] - [(_ . args) (syntax/loc stx (flvector-map . args))] - [_ (syntax/loc stx flvector-map)])) - - ) +(require 'syntax-defs 'defs) -(require 'defs - 'syntax-defs) - -(provide (rename-out [inline-build-flvector build-flvector] - [inline-flvector-map flvector-map])) diff --git a/collects/math/private/flonum/flvector.rkt b/collects/math/private/flonum/flvector.rkt index c07c87be37..dcf5f47720 100644 --- a/collects/math/private/flonum/flvector.rkt +++ b/collects/math/private/flonum/flvector.rkt @@ -1,10 +1,8 @@ #lang typed/racket/base (require racket/fixnum - racket/string - (for-syntax racket/base syntax/parse) - "flonum-functions.rkt" "../unsafe.rkt" + "flonum-functions.rkt" "flvector-syntax.rkt") (provide @@ -19,38 +17,22 @@ flvector->vector ;; Pointwise operations flvector-scale - flvector-round - flvector-floor - flvector-ceiling - flvector-truncate - flvector-abs flvector-sqr flvector-sqrt - flvector-log - flvector-exp - flvector-sin - flvector-cos - flvector-tan - flvector-asin - flvector-acos - flvector-atan + flvector-abs flvector+ flvector* flvector- flvector/ - flvector-expt flvector-min flvector-max - flvector= - flvector< - flvector<= - flvector> - flvector>= - ;; - flvector-sums) + ;; Sum + flvector-sum + flvector-sums + flsum) ;; =================================================================================================== -;; flvector-copy +;; flvector-copy! (: unsafe-flvector-copy! (FlVector Integer FlVector Integer Integer -> Void)) (define (unsafe-flvector-copy! dest dest-start src src-start src-end) @@ -92,12 +74,8 @@ (: list->flvector ((Listof Real) -> FlVector)) (define (list->flvector vs) - (define n (length vs)) - (define xs (make-flvector n)) - (let loop ([#{i : Nonnegative-Fixnum} 0] [vs vs]) - (cond [(i . < . n) (unsafe-flvector-set! xs i (real->double-flonum (unsafe-car vs))) - (loop (+ i 1) (unsafe-cdr vs))] - [else xs]))) + (for/flvector: #:length (length vs) ([v (in-list vs)]) + (fl v))) (: flvector->list (FlVector -> (Listof Float))) (define (flvector->list xs) @@ -105,148 +83,156 @@ (: vector->flvector ((Vectorof Real) -> FlVector)) (define (vector->flvector vs) - (define n (vector-length vs)) - (define xs (make-flvector n)) - (let loop ([#{i : Nonnegative-Fixnum} 0]) - (cond [(i . < . n) (unsafe-flvector-set! xs i (real->double-flonum (unsafe-vector-ref vs i))) - (loop (+ i 1))] - [else xs]))) + (for/flvector: #:length (vector-length vs) ([v (in-vector vs)]) + (fl v))) (: flvector->vector (FlVector -> (Vectorof Float))) (define (flvector->vector xs) - (define n (flvector-length xs)) - (define vs (make-vector n 0.0)) - (let loop ([#{i : Nonnegative-Fixnum} 0]) - (cond [(i . < . n) (unsafe-vector-set! vs i (unsafe-flvector-ref xs i)) - (loop (+ i 1))] - [else vs]))) + (for/vector: #:length (flvector-length xs) ([x (in-flvector xs)]) : Flonum + x)) ;; =================================================================================================== ;; Pointwise operations -(define-syntax (lift1 stx) - (syntax-case stx () - [(_ f) (syntax/loc stx (λ (arr) (flvector-map f arr)))])) +(define-syntax-rule (lift1 f) + (λ: ([arr : FlVector]) + (inline-flvector-map f arr))) -(define-syntax (lift2 stx) - (syntax-case stx () - [(_ f) (syntax/loc stx (λ (arr1 arr2) (flvector-map f arr1 arr2)))])) - -(define-syntax-rule (lift-comparison name comp) - (λ (xs1 xs2) - (define n1 (flvector-length xs1)) - (define n2 (flvector-length xs2)) - (unless (= n1 n2) (error name "flvectors must be the same length; given lengths ~e and ~e" n1 n2)) - (build-vector - n1 (λ: ([j : Index]) - (comp (unsafe-flvector-ref xs1 j) - (unsafe-flvector-ref xs2 j)))))) +(define-syntax-rule (lift2 f) + (λ: ([arr0 : FlVector] [arr1 : FlVector]) + (inline-flvector-map f arr0 arr1))) (: flvector-scale (FlVector Float -> FlVector)) -(define (flvector-scale arr y) (flvector-map (λ (x) (fl* x y)) arr)) +(define (flvector-scale arr y) (inline-flvector-map (λ: ([x : Flonum]) (fl* x y)) arr)) + +(: flvector-sqr (FlVector -> FlVector)) +(define flvector-sqr (lift1 (λ: ([x : Flonum]) (fl* x x)))) -(: flvector-round (FlVector -> FlVector)) -(: flvector-floor (FlVector -> FlVector)) -(: flvector-ceiling (FlVector -> FlVector)) -(: flvector-truncate (FlVector -> FlVector)) -(: flvector-abs (FlVector -> FlVector)) -(: flvector-sqr (FlVector -> FlVector)) (: flvector-sqrt (FlVector -> FlVector)) -(: flvector-log (FlVector -> FlVector)) -(: flvector-exp (FlVector -> FlVector)) -(: flvector-sin (FlVector -> FlVector)) -(: flvector-cos (FlVector -> FlVector)) -(: flvector-tan (FlVector -> FlVector)) -(: flvector-asin (FlVector -> FlVector)) -(: flvector-acos (FlVector -> FlVector)) -(: flvector-atan (FlVector -> FlVector)) +(define flvector-sqrt (lift1 flsqrt)) + +(: flvector-abs (FlVector -> FlVector)) +(define flvector-abs (lift1 flabs)) (: flvector+ (FlVector FlVector -> FlVector)) -(: flvector* (FlVector FlVector -> FlVector)) -(: flvector- (case-> (FlVector -> FlVector) - (FlVector FlVector -> FlVector))) -(: flvector/ (case-> (FlVector -> FlVector) - (FlVector FlVector -> FlVector))) -(: flvector-expt (FlVector FlVector -> FlVector)) -(: flvector-min (FlVector FlVector -> FlVector)) -(: flvector-max (FlVector FlVector -> FlVector)) - -(: flvector= (FlVector FlVector -> (Vectorof Boolean))) -(: flvector< (FlVector FlVector -> (Vectorof Boolean))) -(: flvector<= (FlVector FlVector -> (Vectorof Boolean))) -(: flvector> (FlVector FlVector -> (Vectorof Boolean))) -(: flvector>= (FlVector FlVector -> (Vectorof Boolean))) - -(define flvector-round (lift1 flround)) -(define flvector-floor (lift1 flfloor)) -(define flvector-ceiling (lift1 flceiling)) -(define flvector-truncate (lift1 fltruncate)) -(define flvector-abs (lift1 flabs)) -(define flvector-sqr (lift1 (λ: ([x : Float]) (fl* x x)))) -(define flvector-sqrt (lift1 flsqrt)) -(define flvector-log (lift1 fllog)) -(define flvector-exp (lift1 flexp)) -(define flvector-sin (lift1 flsin)) -(define flvector-cos (lift1 flcos)) -(define flvector-tan (lift1 fltan)) -(define flvector-asin (lift1 flasin)) -(define flvector-acos (lift1 flacos)) -(define flvector-atan (lift1 flatan)) - (define flvector+ (lift2 fl+)) + +(: flvector* (FlVector FlVector -> FlVector)) (define flvector* (lift2 fl*)) +(: flvector- (case-> (FlVector -> FlVector) + (FlVector FlVector -> FlVector))) (define flvector- - (case-lambda - [(arr) (flvector-map (λ: ([x : Float]) (fl- 0.0 x)) arr)] - [(arr1 arr2) (flvector-map fl- arr1 arr2)])) + (case-lambda: + [([arr0 : FlVector]) + (inline-flvector-map (λ: ([x : Float]) (fl- 0.0 x)) arr0)] + [([arr0 : FlVector] [arr1 : FlVector]) + (inline-flvector-map fl- arr0 arr1)])) +(: flvector/ (case-> (FlVector -> FlVector) + (FlVector FlVector -> FlVector))) (define flvector/ - (case-lambda - [(arr) (flvector-map (λ: ([x : Float]) (fl/ 1.0 x)) arr)] - [(arr1 arr2) (flvector-map fl/ arr1 arr2)])) + (case-lambda: + [([arr0 : FlVector]) + (inline-flvector-map (λ: ([x : Float]) (fl/ 1.0 x)) arr0)] + [([arr0 : FlVector] [arr1 : FlVector]) + (inline-flvector-map fl/ arr0 arr1)])) -(define flvector-expt (lift2 flexpt)) -(define flvector-min (lift2 flmin)) -(define flvector-max (lift2 flmax)) +(: flvector-min (FlVector FlVector -> FlVector)) +(define flvector-min (lift2 flmin)) -(define flvector= (lift-comparison 'flvector= fl=)) -(define flvector< (lift-comparison 'flvector< fl<)) -(define flvector<= (lift-comparison 'flvector<= fl<=)) -(define flvector> (lift-comparison 'flvector> fl>)) -(define flvector>= (lift-comparison 'flvector>= fl>=)) +(: flvector-max (FlVector FlVector -> FlVector)) +(define flvector-max (lift2 flmax)) ;; =================================================================================================== +;; Summation + +#| +Algorithm adapted from: + +J R Shewchuk. Adaptive Precision Floating-Point Arithmetic and Fast Geometric Predicates. +Discrete & Computational Geometry, 1996, vol 18, pp 305--363. +|# + +(: flvector-sum (FlVector -> Flonum)) +;; Returns the sum of the elements in xs in a way that incurs rounding error only once +(define (flvector-sum xs) + (define n (flvector-length xs)) + ;; Vector of remainders + (define rs (make-flvector n)) + ;; Loop over `xs' + (let i-loop ([#{i : Nonnegative-Fixnum} 0] + ;; p = Number of valid remainders in `rs' + [#{p : Nonnegative-Fixnum} 0]) + (cond + [(i . fx< . n) + ;; Add x consecutively to each remainder, storing the remainder of *those* additions in `rs' + (let j-loop ([#{j : Nonnegative-Fixnum} 0] + ;; q = Number of remainders generated by this j-loop: + [#{q : Nonnegative-Fixnum} 0] + [x (unsafe-flvector-ref xs i)]) + (cond + [(j . fx< . p) + (define r (unsafe-flvector-ref rs j)) + ;; Get the largest of x and r, or x if it's not rational + (let-values ([(x r) (if ((flabs x) . fl< . (flabs r)) (values r x) (values x r))]) + ;; Add with remainder + (define z (fl+ x r)) + (define-values (hi lo) + (cond [(flrational? z) (values z (fl- r (fl- z x)))] + [else (values x r)])) + (cond [(fl= lo 0.0) + ;; No remainder: don't store (makes average case O(n*log(n))) + (j-loop (unsafe-fx+ j 1) q hi)] + [else + ;; Store the remainder, increment the counter + (unsafe-flvector-set! rs q lo) + (j-loop (unsafe-fx+ j 1) (unsafe-fx+ q 1) hi)]))] + [else + ;; Store the sum so far as the last remainder + (unsafe-flvector-set! rs q x) + (i-loop (fx+ i 1) (unsafe-fx+ q 1))]))] + [else + ;; Add all the remainders + (let j-loop ([#{j : Nonnegative-Fixnum} 0] [acc 0.0]) + (cond [(j . fx< . p) (j-loop (unsafe-fx+ j 1) (fl+ acc (unsafe-flvector-ref rs j)))] + [else acc]))]))) (: flvector-sums (FlVector -> FlVector)) +;; Returns the partial sums of the elements in xs in a way that incurs rounding error only once +;; for each +;; This function works just like `flvector-sum', but keeps track of partial sums instead of +;; summing all the remainders at the end (define (flvector-sums xs) (define n (flvector-length xs)) (define rs (make-flvector n)) (define ss (make-flvector n)) - (let j-loop ([#{j : Nonnegative-Fixnum} 0] - [#{m : Nonnegative-Fixnum} 0]) + (let i-loop ([#{i : Nonnegative-Fixnum} 0] + [#{p : Nonnegative-Fixnum} 0]) (cond - [(j . fx< . n) - (define x (unsafe-flvector-ref xs j)) - (let p-loop ([#{p : Nonnegative-Fixnum} 0] - [#{x : Flonum} x] - [#{s : Flonum} 0.0] - [#{i : Nonnegative-Fixnum} 0]) + [(i . fx< . n) + (let j-loop ([#{j : Nonnegative-Fixnum} 0] + [#{q : Nonnegative-Fixnum} 0] + [x (unsafe-flvector-ref xs i)] + [s 0.0]) (cond - [(p . fx< . m) - (define r (unsafe-flvector-ref rs p)) + [(j . fx< . p) + (define r (unsafe-flvector-ref rs j)) (let-values ([(x r) (if ((flabs x) . fl< . (flabs r)) (values r x) (values x r))]) (define z (fl+ x r)) (define-values (hi lo) (cond [(flrational? z) (values z (fl- r (fl- z x)))] [else (values x r)])) (cond [(fl= lo 0.0) - (p-loop (unsafe-fx+ p 1) hi s i)] + (j-loop (unsafe-fx+ j 1) q hi s)] [else - (unsafe-flvector-set! rs i lo) - (p-loop (unsafe-fx+ p 1) hi (fl+ s lo) (unsafe-fx+ i 1))]))] + (unsafe-flvector-set! rs q lo) + (j-loop (unsafe-fx+ j 1) (unsafe-fx+ q 1) hi (fl+ s lo))]))] [else - (unsafe-flvector-set! rs i x) - (unsafe-flvector-set! ss j (fl+ s x)) - (j-loop (fx+ j 1) (unsafe-fx+ i 1))]))] + (unsafe-flvector-set! rs q x) + (unsafe-flvector-set! ss i (fl+ s x)) + (i-loop (fx+ i 1) (unsafe-fx+ q 1))]))] [else ss]))) + +(: flsum ((Listof Flonum) -> Flonum)) +(define (flsum xs) (flvector-sum (list->flvector xs))) diff --git a/collects/math/parameters.rkt b/collects/math/private/parameters.rkt similarity index 100% rename from collects/math/parameters.rkt rename to collects/math/private/parameters.rkt diff --git a/collects/math/private/utils.rkt b/collects/math/private/utils.rkt index 29223b75d3..12ce6fb503 100644 --- a/collects/math/private/utils.rkt +++ b/collects/math/private/utils.rkt @@ -1,8 +1,12 @@ #lang typed/racket/base -(require racket/pretty) +(require racket/pretty + racket/fixnum + racket/flonum + "exception.rkt") -(provide pretty-print-constructor) +(provide pretty-print-constructor + check-flvector-lengths!) (: port-next-column (Output-Port -> Natural)) ;; Helper to avoid the annoying #f column value @@ -59,3 +63,9 @@ [else ;; No pretty printer, or printing to infinite-width lines, so print on one line (print-all port 'one-line)])) + +(: check-flvector-lengths! (Symbol Index FlVector * -> Void)) +(define (check-flvector-lengths! name n . xss) + (for: ([xs (in-list xss)]) + (unless (fx= n (flvector-length xs)) + (raise-length-error name "FlVector" xs n)))) diff --git a/collects/math/private/vector/vector-fft.rkt b/collects/math/private/vector/vector-fft.rkt index 9232c702da..bbb38b299c 100644 --- a/collects/math/private/vector/vector-fft.rkt +++ b/collects/math/private/vector/vector-fft.rkt @@ -6,7 +6,7 @@ racket/future "../../base.rkt" "../../flonum.rkt" - "../../parameters.rkt" + "../parameters.rkt" "../unsafe.rkt") (provide vector-fft flvector-fft! diff --git a/collects/math/scribblings/math-array.scrbl b/collects/math/scribblings/math-array.scrbl index 194c86cbba..95a2990b03 100644 --- a/collects/math/scribblings/math-array.scrbl +++ b/collects/math/scribblings/math-array.scrbl @@ -7,8 +7,9 @@ (only-in typed/racket/base ann inst : λ: define: make-predicate -> Flonum Real Boolean Any Integer Index Natural Exact-Positive-Integer - Nonnegative-Real Sequenceof Fixnum Values Number - All U List Vector Listof Vectorof Struct)) + Nonnegative-Real Sequenceof Fixnum Values Number Float-Complex + All U List Vector Listof Vectorof Struct FlVector + Symbol Output-Port)) "utils.rkt") @(define typed-eval (make-math-eval)) @@ -30,15 +31,20 @@ For now, if you need speed, use the @racketmodname[typed/racket] language. @defmodule[math/array] -One of the most common ways to structure data is with an array: a grid of homogeneous, -independent elements, usually consisting of rows and columns. But an array data type +One of the most common ways to structure data is with an array: a rectangular grid of +homogeneous, independent elements. But an array data type is usually absent from functional languages' libraries. This is probably because arrays are perceived as requiring users to operate on them using destructive updates, write loops that micromanage array elements, and in general, stray far from the declarative ideal. -@margin-note{TODO: Cite Haskell array paper} -Normally, they do. However, experience in Python, and more recently Data-Parallel Haskell, +@(define array-pdf + "http://research.microsoft.com/en-us/um/people/simonpj/papers/ndp/RArrays.pdf") + +@margin-note*{* @bold{Regular, shape-polymorphic, parallel arrays in Haskell}, + Gabriele Keller, Manuel Chakravarty, Roman Leshchinskiy, Simon Peyton Jones, + and Ben Lippmeier. ICFP 2010. @hyperlink[array-pdf]{(PDF)}} +Normally, they do. However, experience in Python, and more recently Data-Parallel Haskell*, has shown that providing the right data types and a rich collection of whole-array operations allows working effectively with arrays in a functional, declarative style. As a bonus, doing so opens the possibility of parallelizing nearly every operation. @@ -277,6 +283,14 @@ because @racket[nums]'s second axis was stretched during broadcasting. Also, the @racket[#[#["aa"] #["ba"] #["ca"]]] in @racket[lets] is repeated because @racket[lets]'s first axis was stretched. +Even more concretely: +@interaction[#:eval typed-eval + (define ds + (array-shape-broadcast (list (array-shape nums) (array-shape lets)))) + ds + (array-broadcast nums ds) + (array-broadcast lets ds)] + @subsubsection[#:tag "array:broadcasting:control"]{Broadcasting Control} The parameter @racket[array-broadcasting] controls how pointwise operations @tech{broadcast} @@ -514,17 +528,17 @@ Returns an array with @tech{shape} @racket[ds] and @tech{procedure} @racket[proc Analogous to @racket[build-vector], but the result is @tech{non-strict}. @examples[#:eval typed-eval (eval:alts - (define: fibs : (Array Exact-Positive-Integer) + (define: fibs : (Array Natural) (build-array #(10) (λ: ([js : Indexes]) (define j (vector-ref js 0)) - (cond [(j . < . 2) 1] + (cond [(j . < . 2) j] [else (+ (array-ref fibs (vector (- j 1))) (array-ref fibs (vector (- j 2))))])))) (void)) (eval:alts fibs - (ann (array #[1 1 2 3 5 8 13 21 34 55]) (Array Exact-Positive-Integer)))] + (ann (array #[0 1 1 2 3 5 8 13 21 34]) (Array Natural)))] Because @racket[build-array] returns a non-strict array, @racket[fibs] may refer to itself within its definition. Of course, this naïve implementation computes its elements in time exponential in the size of @racket[fibs]. A quick, widely applicable fix is given in @@ -550,18 +564,18 @@ The example in @racket[build-array]'s documentation computes the Fibonacci numbe time. Speeding it up to linear time only requires wrapping its definition with @racket[array-lazy]: @interaction[#:eval typed-eval (eval:alts - (define: fibs : (Array Exact-Positive-Integer) + (define: fibs : (Array Natural) (array-lazy (build-array #(10) (λ: ([js : Indexes]) (define j (vector-ref js 0)) - (cond [(j . < . 2) 1] + (cond [(j . < . 2) j] [else (+ (array-ref fibs (vector (- j 1))) (array-ref fibs (vector (- j 2))))]))))) (void)) (eval:alts fibs - (ann (array #[1 1 2 3 5 8 13 21 34 55]) (Array Exact-Positive-Integer)))] + (ann (array #[0 1 1 2 3 5 8 13 21 34]) (Array Natural)))] Printing a lazy array computes and caches all of its elements, as does applying @racket[array-strict] to it. } @@ -633,6 +647,16 @@ have the value @racket[on-value]; the rest have the value @racket[off-value]. @section{Conversion} +@defform[(Listof* A)]{ +Equivalent to @racket[(U A (Listof A) (Listof (Listof A)) ...)] if infinite unions were allowed. +This is used as an argument type to @racket[list*->array] and as the return type of +@racket[array->list*]. +} + +@defform[(Vectorof* A)]{ +Like @racket[(Listof* A)], but for vectors. See @racket[vector*->array] and @racket[array->vector*]. +} + @deftogether[(@defproc[(list->array [lst (Listof A)]) (Array A)] @defproc[(array->list [arr (Array A)]) (Listof A)])]{ Convert lists to single-axis arrays and back. If @racket[arr] has no axes or more than one axis, @@ -652,12 +676,6 @@ For conversion between nested lists and multidimensional arrays, see @racket[lis Like @racket[list->array] and @racket[array->list], but for vectors. } -@defform[(Listof* A)]{ -Equivalent to @racket[(U A (Listof A) (Listof (Listof A)) ...)] if infinite unions were allowed. -This is used as an argument type to @racket[list*->array] and as the return type of -@racket[array->list*]. -} - @defproc[(list*->array [lsts (Listof* A)] [pred? ((Listof* A) -> Any : A)]) (Array A)]{ Converts a nested list of elements of type @racket[A] to an array. The predicate @racket[pred?] identifies elements of type @racket[A]. The shape of @racket[lsts] must be rectangular. @@ -676,10 +694,6 @@ without it, there is no way to distinguish between rows and elements. The inverse of @racket[list*->array]. } -@defform[(Vectorof* A)]{ -Like @racket[(Listof* A)], but for vectors. See @racket[vector*->array] and @racket[array->vector*]. -} - @defproc[(vector*->array [vecs (Vectorof* A)] [pred? ((Vectorof* A) -> Any : A)]) (Array A)]{ Like @racket[list*->array], but accepts nested vectors of elements. @examples[#:eval typed-eval @@ -725,10 +739,12 @@ The axis number @racket[axis] must be nonnegative and less than the number of @r @defparam[array-custom-printer print-array - (All (A) ((Array A) Symbol Output-Port (U Boolean Zero One) -> Any))]{ -A parameter that controls array printing. + (All (A) ((Array A) + Symbol + Output-Port + (U Boolean 0 1) -> Any))]{ +A parameter whose value is used to print subtypes of @racket[Array]. } -(All (A) ((Array A) Symbol Output-Port (U Boolean Zero One) -> Any)) @defproc[(print-array [arr (Array A)] [name Symbol] @@ -737,16 +753,19 @@ A parameter that controls array printing. Any]{ Prints an array using @racket[array] syntax, using @racket[name] instead of @racket['array] as the head form. This function is set as the value of @racket[array-custom-printer] when -@racketmodname[math/array] is required. +@racketmodname[math/array] is first required. Well-behaved @racket[Array] subtypes do not call this function directly to print themselves. They call the current @racket[array-custom-printer]: @interaction[#:eval typed-eval - ((array-custom-printer) - (array #[0 1 2 3]) - 'my-cool-array - (current-output-port) - #t)] + (eval:alts + ((array-custom-printer) + (array #[0 1 2 3]) + 'my-cool-array + (current-output-port) + #t) + (eval:result "" "(my-cool-array #[0 1 2 3])"))] +See @racket[prop:custom-write] for the meaning of the @racket[port] and @racket[mode] arguments. } @@ -846,6 +865,9 @@ an error. @examples[#:eval typed-eval (array-map (λ: ([x : String]) (string-append x "!")) (array #[#["Hello" "I"] #["Am" "Shouting"]])) + (array-map string-append + (array #[#["Hello" "I"] #["Am" "Shouting"]]) + (array "!")) (array-map + (index-array #(3 3 3)) (array 2)) (array-map + (index-array #(2 2)) (index-array #(3 3)))] @@ -874,66 +896,6 @@ optimizer to replace @racket[f] with something unchecked and type-specific (for replace @racket[*] with @racket[unsafe-fl*]), at the expense of code size. } -@deftogether[(@defform[(array-abs arr)] - @defform[(array-round arr)] - @defform[(array-floor arr)] - @defform[(array-ceiling arr)] - @defform[(array-truncate arr)] - @defform[(array-sqrt arr)] - @defform[(array-sqr arr)] - @defform[(array-log arr)] - @defform[(array-exp arr)] - @defform[(array-sin arr)] - @defform[(array-cos arr)] - @defform[(array-tan arr)] - @defform[(array-asin arr)] - @defform[(array-acos arr)] - @defform[(array-atan arr)] - @defform[(array-inexact->exact arr)] - @defform[(array-exact->inexact arr)] - @defform[(array-fl arr)])]{ -Equivalent to @racket[(array-map f arr)], where @racket[f] is respectively -@racket[abs], -@racket[round], -@racket[floor], -@racket[ceiling], -@racket[truncate], -@racket[sqrt], -@racket[sqr], -@racket[log], -@racket[exp], -@racket[sin], -@racket[cos], -@racket[tan], -@racket[asin], -@racket[acos], -@racket[atan], -@racket[inexact->exact], -@racket[exact->inexact], -or @racket[fl]. - -Instead of @racket[array-fl], you might consider using @racket[array->flarray], which returns -a mutable array of flonums, stored unboxed in an flvector. -} - -@deftogether[(@defform[(array-magnitude arr)] - @defform[(array-angle arr)] - @defform[(array-conjugate arr)] - @defform[(array-real-part arr)] - @defform[(array-imag-part arr)] - @defform[(array-fc arr)])]{ -Equivalent to @racket[(array-map f arr)], where @racket[f] is respectively -@racket[magnitude], -@racket[angle], -@racket[conjugate], -@racket[real-part], -@racket[imag-part], -or @racket[number->float-complex]. - -Instead of @racket[array-fc], you might consider using @racket[array->fcarray], which returns -a mutable array of float-complex numbers, stored unboxed in a pair of flvectors. -} - @deftogether[(@defform[(array+ arrs ...)] @defform[(array* arrs ...)] @defform[(array- arr0 arrs ...)] @@ -953,10 +915,24 @@ returns sensible answers, so do @racket[(array+)] and @racket[(array*)]. Equivalent to @racket[(array* arr (array x))], but faster. } -@deftogether[(@defform[(array-expt arr0 arr1)] - @defform[(array-make-rectangular arr0 arr1)])]{ -Equivalent to @racket[(array-map expt arr0 arr1)] and -@racket[(array-map make-rectangular arr0 arr1)], respectively. +@deftogether[(@defform[(array-abs arr)] + @defform[(array-sqr arr)] + @defform[(array-sqrt arr)] + @defform[(array-conjugate arr)])]{ +Equivalent to @racket[(array-map f arr)], where @racket[f] is respectively +@racket[abs], +@racket[sqr], +@racket[sqrt], +or @racket[conjugate]. +} + +@deftogether[(@defform[(array-real-part arr)] + @defform[(array-imag-part arr)] + @defform[(array-make-rectangular arr0 arr1)] + @defform[(array-magnitude arr)] + @defform[(array-angle arr)] + @defform[(array-make-polar arr0 arr1)])]{ +Conversions to and from complex numbers, lifted to operate on arrays. } @deftogether[(@defform[(array< arr0 arr1 arrs ...)] @@ -972,7 +948,7 @@ Equivalent to @racket[(array-map f arr0 arr1 arrs ...)], where @racket[f] is res @defform[(array-and arr ...)] @defform[(array-or arr ...)] @defform[(array-if cond-arr true-arr false-err)])]{ -@racket[not], @racket[and], @racket[or] and @racket[if] lifted to operate on arrays. +Boolean operators lifted to operate on arrays. The short-cutting behavior of @racket[array-and], @racket[array-or] and @racket[array-if] can keep array arguments' elements from being referred to (and thus computed). However, @@ -1061,6 +1037,7 @@ If @racket[idxs] and @racket[vals] do not have the same shape, they are @tech{br (define idxs (array #['#(0 0) '#(1 1)])) (array-indexes-set! arr idxs (array -1)) arr] +When indexes are repeated in @racket[idxs], they are mutated in some unspecified order. } @defproc[(array-slice-ref [arr (Array A)] [specs (Listof Slice-Spec)]) (Array A)]{ @@ -1068,17 +1045,19 @@ Returns a transformation of @racket[arr] according to the list of slice specific @racket[specs]. See @secref{array:slice-specs} for documentation and examples. } -@defproc[(slice-indexes-array [ds In-Indexes] [specs (Listof Slice-Spec)]) (Array Indexes)]{ -Returns the indexes of the elements that would be retrieved if an array with shape @racket[ds] -were sliced according to @racket[specs]. -Equivalent to @racketblock[(array-slice-ref (indexes-array ds) specs)] -} - @defproc[(array-slice-set! [arr (Settable-Array A)] [specs (Listof Slice-Spec)] [vals (Array A)]) Void]{ Like @racket[array-indexes-set!], but for slice specifications. Equivalent to @racketblock[(let ([idxs (slice-indexes-array (array-shape arr) specs)]) (array-indexes-set! arr idxs vals))] +When a slice specification refers to an element in @racket[arr] more than once, the element is +mutated more than once in some unspecified order. +} + +@defproc[(slice-indexes-array [ds In-Indexes] [specs (Listof Slice-Spec)]) (Array Indexes)]{ +Returns the indexes of the elements that would be retrieved if an array with shape @racket[ds] +were sliced according to @racket[specs]. +Equivalent to @racketblock[(array-slice-ref (indexes-array ds) specs)] } @subsection[#:tag "array:slice-specs"]{Slice Specifications} @@ -1217,10 +1196,10 @@ Removing the second axis by collapsing it to the row with index @racket[1]: Removing the second-to-last axis (which for @racket[arr] is the same as the second): @interaction[#:eval typed-eval (array-slice-ref arr (list ::... 1 (::)))] -All of these examples can be done using @racket[array-axis-ref], but including it in slice -specifications can be handy: -@interaction[#:eval typed-eval - (array-slice-ref arr (list ::... 1 (:: #f #f 2)))] +All of these examples can be done using @racket[array-axis-ref]. However, removing an axis +relative to the dimension of the array (e.g. the second-to-last axis) is easier to do +using @racket[array-slice-ref], and it is sometimes convenient to combine axis removal with +other slice operations. @subsubsection{@racket[Slice-New-Axis]: add an axis} @@ -1570,10 +1549,31 @@ Every fold, including @racket[array-axis-fold], is ultimately defined using @section{Other Array Operations} -@;{array-fft -array-inverse-fft -array-axis-fft -array-axis-inverse-fft +@subsection{Fast Fourier Transform} + +@margin-note{Wikipedia: @hyperlink["http://wikipedia.org/wiki/Discrete_Fourier_transform"]{Discrete +Fourier Transform}} +@defproc[(array-axis-fft [arr (Array Number)] [k Integer]) (Array Float-Complex)]{ +Performs a discrete Fourier transform on axis @racket[k] of @racket[arr]. The length of @racket[k] +must be an integer power of two. (See @racket[power-of-two?].) The scaling convention is +determined by the parameter @racket[dft-convention], which defaults to the convention used in +signal processing. + +The transform is done in parallel using @racket[max-math-threads] threads. +} + +@defproc[(array-axis-inverse-fft [arr (Array Number)] [k Integer]) (Array Float-Complex)]{ +The inverse of @racket[array-axis-fft], performed by parameterizing the forward transform on +@racket[(dft-inverse-convention)]. +} + +@defproc[(array-fft [arr (Array Number)]) FCArray]{ +Performs a discrete Fourier transform on each axis of @racket[arr] using @racket[array-axis-fft]. +} + +@defproc[(array-inverse-fft [arr (Array Number)]) FCArray]{ +The inverse of @racket[array-fft], performed by parameterizing the forward transform on +@racket[(dft-inverse-convention)]. } @@ -1585,100 +1585,143 @@ array-axis-inverse-fft @subsection{Flonum Arrays} @defidform[FlArray]{ +The type of @deftech{flonum arrays}, a subtype of @racket[(Settable-Array Flonum)] that stores its +elements in an @racket[FlVector]. A flonum array is always strict. +} + +@defform[(flarray #[#[...] ...])]{ +Like @racket[array], but creates @tech{flonum arrays}. The listed elements must be real numbers, +and may be exact. + +@examples[#:eval typed-eval + (flarray 0.0) + (flarray #['x]) + (flarray #[#[1 2] #[3 4]])] } @defproc[(array->flarray [arr (Array Real)]) FlArray]{ -Returns an flonum array that has approximately the same elements as @racket[arr]. -The elements may lose precision during the conversion. +Returns a flonum array that has approximately the same elements as @racket[arr]. Exact +elements will likely lose precision during conversion. } -@;{ -FlArray -flarray-scale -flarray-round -flarray-floor -flarray-map -flarray-truncate -flarray-abs -flarray-sqr -flarray-sqrt -flarray-log -flarray-exp -flarray-sin -flarray-cos -flarray-tan -flarray-asin -flarray-acos -flarray-atan -flarray+ -flarray* -flarray-data -flarray -flarray-ceiling -array->flarray -inline-flarray-map -flarray- -flarray/ -flarray-expt -flarray-min -flarray-max -flarray= -flarray< -flarray<= -flarray> -flarray>= +@defproc[(flarray-data [arr FlArray]) FlVector]{ +Returns the elements of @racket[arr] in a flonum vector, in row-major order. +@examples[#:eval typed-eval + (flvector->list (flarray-data (flarray #[#[1 2] #[3 4]])))] +} + +@defproc[(flarray-map [f (Flonum ... -> Flonum)] [arrs FlArray] ...) FlArray]{ +Maps the function @racket[f] over the arrays @racket[arrs]. If the arrays do not have the same +shape, they are @tech{broadcast} first. If the arrays do have the same shape, this operation can +be quite fast. + +The function @racket[f] is meant to accept the same number of arguments as the number of its +following flonum array arguments. However, a current limitation in Typed Racket requires @racket[f] +to accept @italic{any} number of arguments. To map a single-arity function such as @racket[fl+], +for now, use @racket[inline-flarray-map] or @racket[array-map]. +} + +@defform[(inline-flarray-map f arrs ...) + #:contracts ([f (Flonum ... -> Flonum)] + [arrs FlArray])]{ +Like @racket[inline-array-map], but for flonum arrays. + +@bold{This is currently unavailable in untyped Racket.} +} + +@deftogether[(@defproc[(flarray+ [arr0 FlArray] [arr1 FlArray]) FlArray] + @defproc[(flarray* [arr0 FlArray] [arr1 FlArray]) FlArray] + @defproc*[([(flarray- [arr FlArray]) FlArray] + [(flarray- [arr0 FlArray] [arr1 FlArray]) FlArray])] + @defproc*[([(flarray/ [arr FlArray]) FlArray] + [(flarray/ [arr0 FlArray] [arr1 FlArray]) FlArray])] + @defproc[(flarray-min [arr0 FlArray] [arr1 FlArray]) FlArray] + @defproc[(flarray-max [arr0 FlArray] [arr1 FlArray]) FlArray] + @defproc[(flarray-scale [arr FlArray] [x Flonum]) FlArray] + @defproc[(flarray-abs [arr FlArray]) FlArray] + @defproc[(flarray-sqr [arr FlArray]) FlArray] + @defproc[(flarray-sqrt [arr FlArray]) FlArray])]{ +Arithmetic lifted to flonum arrays. } @subsection{Float-Complex Arrays} @defidform[FCArray]{ +The type of @deftech{float-complex arrays}, a subtype of @racket[(Settable-Array Float-Complex)] +that stores its elements in a pair of @racket[FlVector]s. A float-complex array is always strict. +} + +@defform[(fcarray #[#[...] ...])]{ +Like @racket[array], but creates @tech{float-complex arrays}. The listed elements must be numbers, +and may be exact. + +@examples[#:eval typed-eval + (fcarray 0.0) + (fcarray #['x]) + (fcarray #[#[1 2+1i] #[3 4+3i]])] } @defproc[(array->fcarray [arr (Array Number)]) FCArray]{ -Returns a float-complex array that has approximately the same elements as @racket[arr]. -The elements may lose precision during the conversion. +Returns a float-complex array that has approximately the same elements as @racket[arr]. Exact +elements will likely lose precision during conversion. } -@;{ -FCArray -fcarray-fft -inline-fcarray-map -fcarray-imag-data -fcarray-real-data -fcarray -array->fcarray -fcarray-scale -fcarray-sqr -fcarray-sqrt -fcarray-conjugate -fcarray-magnitude -fcarray-angle -fcarray-exp -fcarray-sin -fcarray-cos -fcarray-tan -fcarray-acos -fcarray+ -fcarray* -fcarray- -fcarray/ -fcarray-expt -fcarray-real-part -fcarray-imag-part -fcarray-make-rectangular -fcarray-map -fcarray-log -fcarray-asin -fcarray-atan +@deftogether[(@defproc[(fcarray-real-data [arr FCArray]) FlVector] + @defproc[(fcarray-imag-data [arr FCArray]) FlVector])]{ +Return the real and imaginary parts of @racket[arr]'s elements in flonum vectors, in row-major order. +@examples[#:eval typed-eval + (define arr (fcarray #[#[1 2+1i] #[3 4+3i]])) + (flvector->list (fcarray-real-data arr)) + (flvector->list (fcarray-imag-data arr))] } +@defproc[(fcarray-map [f (Float-Complex ... -> Float-Complex)] [arrs FCArray] ...) FCArray]{ +Maps the function @racket[f] over the arrays @racket[arrs]. If the arrays do not have the same +shape, they are @tech{broadcast} first. If the arrays do have the same shape, this operation can +be quite fast. + +The function @racket[f] is meant to accept the same number of arguments as the number of its +following float-complex array arguments. However, a current limitation in Typed Racket requires +@racket[f] to accept @italic{any} number of arguments. To map a single-arity function, for now, +use @racket[inline-fcarray-map] or @racket[array-map]. +} + +@defform[(inline-fcarray-map f arrs ...) + #:contracts ([f (Float-Complex ... -> Float-Complex)] + [arrs FCArray])]{ +Like @racket[inline-array-map], but for float-complex arrays. + +@bold{This is currently unavailable in untyped Racket.} +} + +@deftogether[(@defproc[(fcarray+ [arr0 FCArray] [arr1 FCArray]) FCArray] + @defproc[(fcarray* [arr0 FCArray] [arr1 FCArray]) FCArray] + @defproc*[([(fcarray- [arr FCArray]) FCArray] + [(fcarray- [arr0 FCArray] [arr1 FCArray]) FCArray])] + @defproc*[([(fcarray/ [arr FCArray]) FCArray] + [(fcarray/ [arr0 FCArray] [arr1 FCArray]) FCArray])] + @defproc[(fcarray-scale [arr FCArray] [z Float-Complex]) FCArray] + @defproc[(fcarray-sqr [arr FCArray]) FCArray] + @defproc[(fcarray-sqrt [arr FCArray]) FCArray] + @defproc[(fcarray-conjugate [arr FCArray]) FCArray])]{ +Arithmetic lifted to float-complex arrays. +} + +@deftogether[(@defproc[(fcarray-real-part [arr FCArray]) FlArray] + @defproc[(fcarray-imag-part [arr FCArray]) FlArray] + @defproc[(fcarray-make-rectangular [arr0 FlArray] [arr1 FlArray]) FCArray] + @defproc[(fcarray-magnitude [arr FCArray]) FlArray] + @defproc[(fcarray-angle [arr FCArray]) FlArray] + @defproc[(fcarray-make-polar [arr0 FlArray] [arr1 FlArray]) FCArray])]{ +Conversions to and from complex numbers, lifted to flonum and float-complex arrays. +} @;{==================================================================================================} +@;{ @section{Unsafe Operations} -@;{ unsafe-array-ref unsafe-array-set! unsafe-array-proc diff --git a/collects/math/scribblings/math-base.scrbl b/collects/math/scribblings/math-base.scrbl index 83bc27473f..de3b450069 100644 --- a/collects/math/scribblings/math-base.scrbl +++ b/collects/math/scribblings/math-base.scrbl @@ -22,7 +22,7 @@ as well as providing the values document below. In general, the functions provided by @racketmodname[math/base] are @deftech{elementary} functions, or those functions that can be defined in terms of a finite number of arithmetic operations, logarithms, exponentials, trigonometric functions, and constants. -For others, see @racketmodname[math/special-functions]. +For others, see @racketmodname[math/special-functions] and @racketmodname[math/distributions]. @section{Constants} diff --git a/collects/math/scribblings/math-bigfloat.scrbl b/collects/math/scribblings/math-bigfloat.scrbl index 9b6d74a078..79b27ee7c5 100644 --- a/collects/math/scribblings/math-bigfloat.scrbl +++ b/collects/math/scribblings/math-bigfloat.scrbl @@ -31,6 +31,8 @@ for Windows and Mac OS X, is installed on most Linux systems, and is @hyperlink["http://www.mpfr.org/ports.html"]{easy to install} on all major platforms. See @secref{loading} for details. +@local-table-of-contents[] + @section[#:tag "quick"]{Quick Start} @itemlist[#:style 'ordered diff --git a/collects/math/scribblings/math-distributions.scrbl b/collects/math/scribblings/math-distributions.scrbl index 00306ee24a..82aa14a8d0 100644 --- a/collects/math/scribblings/math-distributions.scrbl +++ b/collects/math/scribblings/math-distributions.scrbl @@ -216,10 +216,10 @@ ordered distribution returns as random samples. Additionally, its @tech{cdf} acc (ordered-dist? (discrete-dist '(a b c))) (ordered-dist? (normal-dist))] -The median is stored in an @racket[ordered-dist] to allow interval probabilities to be calculated +The median is stored in an @racket[ordered-dist] to allow interval probabilities to be computed accurately. For example, for @racket[d = (normal-dist)], whose median is @racket[0.0], -@racket[(real-dist-prob d -2.0 -1.0)] is calculated using lower-tail probabilities, and -@racket[(real-dist-prob d 1.0 2.0)] is calculated using upper-tail probabilities. +@racket[(real-dist-prob d -2.0 -1.0)] is computed using lower-tail probabilities, and +@racket[(real-dist-prob d 1.0 2.0)] is computed using upper-tail probabilities. } @defidform[Real-Dist]{ @@ -280,6 +280,9 @@ The first four are synonyms for @racket[ordered-dist-cdf], @racket[ordered-dist- (Discrete-Dist A)])] @defproc[(discrete-dist-values [d (Discrete-Dist A)]) (Listof A)] @defproc[(discrete-dist-probs [d (Discrete-Dist A)]) (Listof Positive-Flonum)])]{ +@bold{Warning:} The types and functions for discrete distributions will probably change soon, +as part of overhauling how @racket[math/statistics] handles weighted samples. + Represents families of unordered, discrete distributions over values of type @racket[A], with equality decided by @racket[equal?]. @@ -311,8 +314,7 @@ or over extended integers. The most common definitions use the extended reals, s @tech{distribution object} constructors return objects of type @racket[Real-Dist]. (Another reason is that the extended integers correspond with the type -@racket[(U Integer +inf.0 -inf.0)]. Values of this type are difficult to use as integers: in -practice, they are the same as values of type @racket[Real].) +@racket[(U Integer +inf.0 -inf.0)]. Values of this type have little support in Racket's library.) This leaves us with a quandary and two design decisions users should be aware of. The quandary is that, when an integer distribution is defined over the reals, it has a @tech{cdf}, but @italic{no @@ -422,7 +424,7 @@ independent events. @section{Real Distribution Families} -The @tech{distribution object} constructors documented in this section return unique, useful +The @tech{distribution object} constructors documented in this section return uniquely defined distributions for the largest possible parameter domain. This usually means that they return distributions for a larger domain than their mathematical counterparts are defined on. @@ -523,18 +525,13 @@ Represents the Cauchy distribution family parameterized by mode and scale. Represents the family of distributions whose densities are Dirac delta functions. @examples[#:eval untyped-eval - ((dist-pdf (delta-dist)) 0) - ((dist-pdf (delta-dist)) 1) + (pdf (delta-dist) 0) + (pdf (delta-dist) 1) (plot (for/list ([μ (in-list '(-1 0 1))] [i (in-naturals)]) (function (dist-cdf (delta-dist μ)) #:color i #:style i #:label (format "δ(~a)" μ))) #:x-min -2 #:x-max 2 #:y-label "probability")] -When a distribution with a scale parameter has scale zero, it behaves like a delta distribution: -@interaction[#:eval untyped-eval - ((dist-pdf (normal-dist 0 0)) 0) - ((dist-pdf (normal-dist 0 0)) 1) - (plot (function (dist-cdf (normal-dist 0 0)) -1e-300 1e-300))] } @subsection{Exponential Distributions} @@ -548,7 +545,7 @@ When a distribution with a scale parameter has scale zero, it behaves like a del Represents the exponential distribution family parameterized by mean, or scale. @bold{Warning:} The exponential distribution family is often parameterized by @italic{rate}, which -is the reciprocal of mean or scale. If you have rates, construct exponential distributions using +is the reciprocal of mean or scale. Construct exponential distributions from rates using @racketblock[(exponential-dist (/ 1.0 rate))] @examples[#:eval untyped-eval @@ -579,7 +576,7 @@ Represents the gamma distribution family parameterized by shape and scale. The @ parameter must be nonnegative. @bold{Warning:} The gamma distribution family is often parameterized by shape and @italic{rate}, -which is the reciprocal of scale. If you have rates, construct gamma distributions using +which is the reciprocal of scale. Construct gamma distributions from rates using @racketblock[(gamma-dist shape (/ 1.0 rate))] @examples[#:eval untyped-eval @@ -648,8 +645,7 @@ and scale. In this parameterization, the variance is @racket[(* 1/3 (sqr (* pi s Represents the normal distribution family parameterized by mean and standard deviation. @bold{Warning:} The normal distribution family is often parameterized by mean and @italic{variance}, -which is the square of standard deviation. If you have variances, construct normal distributions -using +which is the square of standard deviation. Construct normal distributions from variances using @racketblock[(normal-dist mean (sqrt var))] @examples[#:eval untyped-eval @@ -776,7 +772,7 @@ distribution centered at @racket[x]. @section[#:tag "dist:flonum"]{Low-Level Distribution Functions} The following functions are provided for users who need lower overhead than that of -@tech{distribution objects}, such as (currently) untyped Racket users, and library writers who +@tech{distribution objects}, such as untyped Racket users (currently), and library writers who are implementing their own distribution abstractions. Because applying these functions is meant to be fast, none of them have optional arguments. In @@ -854,7 +850,7 @@ Low-level flonum functions used to implement @racket[cauchy-dist]. @defproc[(fldelta-inv-cdf [mean Flonum] [p Flonum] [log? Any] [1-p? Any]) Flonum])]{ Low-level flonum functions used to implement @racket[delta-dist]. -If you need delta-distributed random samples, use @racket[(make-flvector n mean)]. +To get delta-distributed random samples, use @racket[(make-flvector n mean)]. } @deftogether[ diff --git a/collects/math/scribblings/math-flonum.scrbl b/collects/math/scribblings/math-flonum.scrbl index ac6f5a238c..90135ffb43 100644 --- a/collects/math/scribblings/math-flonum.scrbl +++ b/collects/math/scribblings/math-flonum.scrbl @@ -2,12 +2,15 @@ @(require scribble/eval racket/sandbox - (for-label racket/base + (for-label racket/base racket/vector racket/list math plot - (only-in typed/racket/base Flonum Real Boolean Any Listof Integer)) + (only-in typed/racket/base + -> + Flonum Integer Index Real Boolean Any Listof Vectorof FlVector)) "utils.rkt") @(define untyped-eval (make-untyped-math-eval)) +@interaction-eval[#:eval untyped-eval (require racket/list)] @title[#:tag "flonum"]{Flonums} @(author-neil) @@ -17,6 +20,8 @@ For convenience, @racketmodname[math/flonum] re-exports @racketmodname[racket/flonum] as well as providing the functions document below. +@local-table-of-contents[] + @section{Additional Flonum Functions} @defproc[(fl [x Real]) Flonum]{ @@ -61,6 +66,8 @@ The @racket[sum] function does the same for heterogenous lists of reals. Worst-case time complexity is O(@italic{n}@superscript{2}), though the pathological inputs needed to observe quadratic time are exponentially improbable and are hard to generate purposely. Expected time complexity is O(@italic{n} log(@italic{n})). + +See @racket[flvector-sums] for a variant that computes all the partial sums in @racket[xs]. } @deftogether[(@defproc[(flsinh [x Flonum]) Flonum] @@ -496,4 +503,89 @@ it is not zero and @racket[((abs x) . <= . +max-subnormal.0)]. @examples[#:eval untyped-eval +max-subnormal.0] } +@section{Additional Flonum Vector Functions} + +@defproc[(build-flvector [n Integer] [proc (Index -> Flonum)]) FlVector]{ +Creates a length-@racket[n] flonum vector by applying @racket[proc] to the indexes +from @racket[0] to @racket[(- n 1)]. Analogous to @racket[build-vector]. +@examples[#:eval untyped-eval + (flvector->list (build-flvector 10 fl))] +} + +@defform[(inline-build-flvector n proc) + #:contracts ([n Integer] + [proc (Index -> Flonum)])]{ +Like @racket[build-flvector], but always inlined. This increases speed at the expense of code size. +} + +@defproc[(flvector-map [proc (Flonum Flonum ... -> Flonum)] [xs FlVector] [xss FlVector] ...) + FlVector]{ +Applies @racket[proc] to the corresponding elements of @racket[xs] and @racket[xss]. Analogous to +@racket[vector-map]. + +The @racket[proc] is meant to accept the same number of arguments as the number of its following +flonum vector arguments. However, a current limitation in Typed Racket requires @racket[proc] +to accept @italic{any} number of arguments. To map a single-arity function such as @racket[fl+] +over the corresponding number of flonum vectors, for now, use @racket[inline-flvector-map]. +} + +@defform[(inline-flvector-map proc xs xss ...) + #:contracts ([proc (Flonum Flonum ... -> Flonum)] + [xs FlVector] + [xss FlVector])]{ +Like @racket[flvector-map], but always inlined. +} + +@defproc[(flvector-copy! [dest FlVector] + [dest-start Integer] + [src FlVector] + [src-start Integer 0] + [src-end Integer (flvector-length src)]) + Void]{ +Like @racket[vector-copy!], but for flonum vectors. +} + +@deftogether[(@defproc[(list->flvector [vs (Listof Real)]) FlVector] + @defproc[(flvector->list [xs FlVector]) (Listof Flonum)] + @defproc[(vector->flvector [vs (Vectorof Real)]) FlVector] + @defproc[(flvector->vector [xs FlVector]) (Vectorof Flonum)])]{ +Convert between lists and flonum vectors, and between vectors and flonum vectors. +} + +@deftogether[(@defproc[(flvector+ [xs FlVector] [ys FlVector]) FlVector] + @defproc[(flvector* [xs FlVector] [ys FlVector]) FlVector] + @defproc*[([(flvector- [xs FlVector]) FlVector] + [(flvector- [xs FlVector] [ys FlVector]) FlVector])] + @defproc*[([(flvector/ [xs FlVector]) FlVector] + [(flvector/ [xs FlVector] [ys FlVector]) FlVector])] + @defproc[(flvector-scale [xs FlVector] [y Flonum]) FlVector] + @defproc[(flvector-abs [xs FlVector]) FlVector] + @defproc[(flvector-sqr [xs FlVector]) FlVector] + @defproc[(flvector-sqrt [xs FlVector]) FlVector] + @defproc[(flvector-min [xs FlVector] [ys FlVector]) FlVector] + @defproc[(flvector-max [xs FlVector] [ys FlVector]) FlVector])]{ +Arithmetic lifted to operate on flonum vectors. +} + +@defproc[(flvector-sum [xs FlVector]) Flonum]{ +Like @racket[flsum], but operates on flonum vectors. In fact, @racket[flsum] is defined in terms +of @racket[flvector-sum]. +} + +@defproc[(flvector-sums [xs FlVector]) FlVector]{ +Computes the partial sums of the elements in @racket[xs] in a way that incurs rounding error only +once for each partial sum. +@examples[#:eval untyped-eval + (flvector->list + (flvector-sums + (flvector 1.0 1e-16 1e-16 1e-16 1e-16 1e100 -1e100)))] +Compare the same example computed by direct summation: +@interaction[#:eval untyped-eval + (rest + (reverse + (foldl (λ (x xs) (cons (+ x (first xs)) xs)) + (list 0.0) + '(1.0 1e-16 1e-16 1e-16 1e-16 1e100 -1e100))))] +} + @(close-eval untyped-eval) diff --git a/collects/math/scribblings/math-special-functions.scrbl b/collects/math/scribblings/math-special-functions.scrbl index 3b4de9b702..221eda6c81 100644 --- a/collects/math/scribblings/math-special-functions.scrbl +++ b/collects/math/scribblings/math-special-functions.scrbl @@ -2,7 +2,7 @@ @(require scribble/eval racket/sandbox - (for-label racket/base + (for-label racket/base racket/unsafe/ops math plot (only-in typed/racket/base Flonum Real Integer Natural Zero Positive-Integer Exact-Rational diff --git a/collects/math/scribblings/math-statistics.scrbl b/collects/math/scribblings/math-statistics.scrbl index 7390b634cd..4e905f7ffb 100644 --- a/collects/math/scribblings/math-statistics.scrbl +++ b/collects/math/scribblings/math-statistics.scrbl @@ -17,10 +17,11 @@ @defmodule[math/statistics] -xxx intro +@bold{Warning:} The documentation for this module is delayed while the +types for weighted samples get reworked. Please wait. -something about accepting weighted samples whenever it makes sense -(time it doesn't make sense: autocorrelation) +@;{something about accepting weighted samples whenever it makes sense +(time it doesn't make sense: autocorrelation)} @local-table-of-contents[] diff --git a/collects/math/scribblings/math-utils.scrbl b/collects/math/scribblings/math-utils.scrbl new file mode 100644 index 0000000000..50569399c4 --- /dev/null +++ b/collects/math/scribblings/math-utils.scrbl @@ -0,0 +1,45 @@ +#lang scribble/manual + +@(require scribble/eval + racket/sandbox + (for-label racket/base racket/future + math + (only-in typed/racket/base + Real Boolean Integer Natural Number Listof + Positive-Flonum Float-Complex Any List Positive-Integer)) + "utils.rkt") + +@(define untyped-eval (make-untyped-math-eval)) + +@title[#:tag "utils"]{Stuff That Doesn't Belong Anywhere Else} +@(author-neil) + +@defmodule[math/utils] + +@defparam[max-math-threads num Positive-Integer]{ +The maximum number of threads a parallelized @racketmodname[math] function +will use. The default value is @racket[(max 1 (processor-count))]. +} + +@defparam[dft-convention lst (List Real Real)]{ +A parameter controlling the convention used for scaling discrete Fourier transforms, such as those +performed by @racket[array-fft]. The default value is @racket['(1 -1)], which represents the convention +used in signal processing. + +In general, if @racket[lst] is @racket[(list a b)] and @racket[n] is the length of a transformed +array axis or vector, then +@itemlist[@item{Each sum is scaled by @racket[(expt n (/ (- a 1) 2))].} + @item{Each exponential in the sum has its argument scaled by @racket[b].}] +Conveniently, a Fourier transform with convention @racket[(list (- a) (- b))] is the inverse +of a Fourier transform with convention @racket[(list a b)]. + +See Mathematica's +@hyperlink["http://reference.wolfram.com/mathematica/tutorial/FourierTransforms.html"]{documentation +on @tt{Fourier}}, from which this excellent idea was stolen. +} + +@defproc[(dft-inverse-convention) (List Real Real)]{ +Returns the convention used for inverse Fourier transforms, given the current convention. +} + +@(close-eval untyped-eval) diff --git a/collects/math/scribblings/math.scrbl b/collects/math/scribblings/math.scrbl index c6c6f9b06d..4969547574 100644 --- a/collects/math/scribblings/math.scrbl +++ b/collects/math/scribblings/math.scrbl @@ -15,10 +15,10 @@ for working with numbers and collections of numbers. These include @item{Special functions} @item{@racket[Bigfloat]s, or arbitrary-precision floating-point numbers} @item{Probability distributions} - @item{Statistical functions} + @item{Statistical functions (currently undergoing refactoring)} @item{Number-theoretic functions} @item{@racket[Array]s for storing and transforming large rectangular data sets} - @item{Linear algebra functions} + @item{Linear algebra functions (currently under review)} ] With this library, we hope to support a wide variety of applied mathematics in @@ -42,3 +42,4 @@ be used in untyped Racket. Exceptions and performance warnings are in @bold{bold @include-section["math-array.scrbl"] @include-section["math-statistics.scrbl"] @include-section["math-distributions.scrbl"] +@include-section["math-utils.scrbl"] diff --git a/collects/math/utils.rkt b/collects/math/utils.rkt new file mode 100644 index 0000000000..d7e563c155 --- /dev/null +++ b/collects/math/utils.rkt @@ -0,0 +1,6 @@ +#lang racket + +(require "private/parameters.rkt") + +(provide (all-from-out + "private/parameters.rkt"))