[racket-dev] [plt] Push #26129: master branch updated
Of course, on part #1, I meant "TR's optimizer replaces them with
`unsafe-fl+' and `unsafe-fl*'...".
Neil ⊥
On 01/21/2013 10:04 PM, ntoronto at racket-lang.org wrote:
> ntoronto has updated `master' from a0f910c3dc to f42cc6f14a.
> http://git.racket-lang.org/plt/a0f910c3dc..f42cc6f14a
>
> =====[ One Commit ]=====================================================
> Directory summary:
> 76.5% collects/math/private/matrix/
> 11.5% collects/math/private/vector/
> 9.7% collects/math/scribblings/
>
> ~~~~~~~~~~
>
> f42cc6f Neil Toronto <ntoronto at racket-lang.org> 2013-01-21 21:55
> :
> | Fixed major performance issue with matrix arithmetic; please merge to 5.3.2
> |
> | The fix consists of three parts:
> |
> | 1. Rewriting `inline-matrix*'. The material change here is that the
> | expansion now contains only direct applications of `+' and `*'.
> | TR's optimizer replaces them with `unsafe-fx+' and `unsafe-fx*',
> | which keeps intermediate flonum values from being boxed.
> |
> | 2. Making the types of all functions that operate on (Matrix Number)
> | values more precise. Now TR can prove that matrix operations preserve
> | inexactness. For example, matrix-conjugate : (Matrix Flonum) ->
> | (Matrix Flonum) and three other cases for Real, Float-Complex, and
> | Number.
> |
> | 3. Changing the return types of some functions that used to return
> | things like (Matrix (U A 0)). Now that we worry about preserving
> | inexactness, we can't have `matrix-upper-triangle' always return a
> | matrix that contains exact zeros. It now accepts an optional `zero'
> | argument of type A.
> :
> M collects/math/matrix.rkt | 22 +++
> M collects/math/private/matrix/matrix-basic.rkt | 173 +++++++++++++------
> M collects/math/private/matrix/matrix-expt.rkt | 22 ++-
> M collects/math/private/matrix/matrix-lu.rkt | 18 +-
> M collects/math/private/matrix/matrix-qr.rkt | 23 ++-
> M collects/math/private/matrix/matrix-solve.rkt | 38 ++--
> M collects/math/private/matrix/matrix-subspace.rkt | 16 +-
> M collects/math/private/matrix/utils.rkt | 66 ++++++-
> M collects/math/private/vector/vector-mutate.rkt | 98 +++++++----
> M collects/math/scribblings/math-matrix.scrbl | 76 ++++----
> M collects/math/tests/matrix-strictness-tests.rkt | 2 +-
> M .../math/private/matrix/matrix-constructors.rkt | 49 ++++--
> M .../math/private/matrix/matrix-gauss-elim.rkt | 24 ++-
> M .../math/private/matrix/matrix-gram-schmidt.rkt | 40 ++++-
> M .../math/private/matrix/matrix-operator-norm.rkt | 23 ++-
> M .../private/matrix/typed-matrix-arithmetic.rkt | 43 +++--
> M .../private/matrix/untyped-matrix-arithmetic.rkt | 92 +++++++---
>
> =====[ Overall Diff ]===================================================
>
> collects/math/matrix.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/matrix.rkt
> +++ NEW/collects/math/matrix.rkt
> @@ -17,6 +17,10 @@
> (except-in "private/matrix/matrix-constructors.rkt"
> vandermonde-matrix)
> (except-in "private/matrix/matrix-basic.rkt"
> + matrix-1norm
> + matrix-2norm
> + matrix-inf-norm
> + matrix-norm
> matrix-dot
> matrix-cos-angle
> matrix-angle
> @@ -29,6 +33,9 @@
> (except-in "private/matrix/matrix-subspace.rkt"
> matrix-col-space)
> (except-in "private/matrix/matrix-operator-norm.rkt"
> + matrix-op-1norm
> + matrix-op-2norm
> + matrix-op-inf-norm
> matrix-basis-cos-angle
> matrix-basis-angle)
> ;;"private/matrix/matrix-qr.rkt" ; all use require/untyped-contract
> @@ -77,6 +84,11 @@
> (require/untyped-contract
> (begin (require "private/matrix/matrix-types.rkt"))
> "private/matrix/matrix-basic.rkt"
> + [matrix-1norm ((Matrix Number) -> Nonnegative-Real)]
> + [matrix-2norm ((Matrix Number) -> Nonnegative-Real)]
> + [matrix-inf-norm ((Matrix Number) -> Nonnegative-Real)]
> + [matrix-norm (case-> ((Matrix Number) -> Nonnegative-Real)
> + ((Matrix Number) Real -> Nonnegative-Real))]
> [matrix-dot
> (case-> ((Matrix Number) -> Nonnegative-Real)
> ((Matrix Number) (Matrix Number) -> Number))]
> @@ -113,6 +125,9 @@
> (require/untyped-contract
> (begin (require "private/matrix/matrix-types.rkt"))
> "private/matrix/matrix-operator-norm.rkt"
> + [matrix-op-1norm ((Matrix Number) -> Nonnegative-Real)]
> + [matrix-op-2norm ((Matrix Number) -> Nonnegative-Real)]
> + [matrix-op-inf-norm ((Matrix Number) -> Nonnegative-Real)]
> [matrix-basis-cos-angle
> ((Matrix Number) (Matrix Number) -> Number)]
> [matrix-basis-angle
> @@ -167,6 +182,10 @@
> ;; matrix-constructors.rkt
> vandermonde-matrix
> ;; matrix-basic.rkt
> + matrix-1norm
> + matrix-2norm
> + matrix-inf-norm
> + matrix-norm
> matrix-dot
> matrix-cos-angle
> matrix-angle
> @@ -179,6 +198,9 @@
> ;; matrix-subspace.rkt
> matrix-col-space
> ;; matrix-operator-norm.rkt
> + matrix-op-1norm
> + matrix-op-2norm
> + matrix-op-inf-norm
> matrix-basis-cos-angle
> matrix-basis-angle
> ;; matrix-qr.rkt
>
> collects/math/private/matrix/matrix-basic.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-basic.rkt
> +++ NEW/collects/math/private/matrix/matrix-basic.rkt
> @@ -130,29 +130,37 @@
> (define i (unsafe-vector-ref js 0))
> (proc ((inst vector Index) i i))))))
>
> -(: matrix-upper-triangle (All (A) ((Matrix A) -> (Matrix (U A 0)))))
> -(define (matrix-upper-triangle M)
> - (define-values (m n) (matrix-shape M))
> - (define proc (unsafe-array-proc M))
> - (array-default-strict
> - (unsafe-build-array
> - ((inst vector Index) m n)
> - (λ: ([ij : Indexes])
> - (define i (unsafe-vector-ref ij 0))
> - (define j (unsafe-vector-ref ij 1))
> - (if (i . fx<= . j) (proc ij) 0)))))
> -
> -(: matrix-lower-triangle (All (A) ((Matrix A) -> (Matrix (U A 0)))))
> -(define (matrix-lower-triangle M)
> - (define-values (m n) (matrix-shape M))
> - (define proc (unsafe-array-proc M))
> - (array-default-strict
> - (unsafe-build-array
> - ((inst vector Index) m n)
> - (λ: ([ij : Indexes])
> - (define i (unsafe-vector-ref ij 0))
> - (define j (unsafe-vector-ref ij 1))
> - (if (i . fx>= . j) (proc ij) 0)))))
> +(: matrix-upper-triangle (All (A) (case-> ((Matrix A) -> (Matrix (U A 0)))
> + ((Matrix A) A -> (Matrix A)))))
> +(define matrix-upper-triangle
> + (case-lambda
> + [(M) (matrix-upper-triangle M 0)]
> + [(M zero)
> + (define-values (m n) (matrix-shape M))
> + (define proc (unsafe-array-proc M))
> + (array-default-strict
> + (unsafe-build-array
> + ((inst vector Index) m n)
> + (λ: ([ij : Indexes])
> + (define i (unsafe-vector-ref ij 0))
> + (define j (unsafe-vector-ref ij 1))
> + (if (i . fx<= . j) (proc ij) zero))))]))
> +
> +(: matrix-lower-triangle (All (A) (case-> ((Matrix A) -> (Matrix (U A 0)))
> + ((Matrix A) A -> (Matrix A)))))
> +(define matrix-lower-triangle
> + (case-lambda
> + [(M) (matrix-lower-triangle M 0)]
> + [(M zero)
> + (define-values (m n) (matrix-shape M))
> + (define proc (unsafe-array-proc M))
> + (array-default-strict
> + (unsafe-build-array
> + ((inst vector Index) m n)
> + (λ: ([ij : Indexes])
> + (define i (unsafe-vector-ref ij 0))
> + (define j (unsafe-vector-ref ij 1))
> + (if (i . fx>= . j) (proc ij) zero))))]))
>
> ;; ===================================================================================================
> ;; Embiggenment (this is a perfectly cromulent word)
> @@ -184,42 +192,63 @@
> ;; ===================================================================================================
> ;; Inner product space (entrywise norm)
>
> -(: matrix-1norm ((Matrix Number) -> Nonnegative-Real))
> -(define (matrix-1norm a)
> +(: nonstupid-magnitude (case-> (Flonum -> Nonnegative-Flonum)
> + (Real -> Nonnegative-Real)
> + (Float-Complex -> Nonnegative-Flonum)
> + (Number -> Nonnegative-Real)))
> +(define (nonstupid-magnitude x)
> + (if (real? x) (abs x) (magnitude x)))
> +
> +(: matrix-1norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Number) -> Nonnegative-Real)))
> +(define (matrix-1norm M)
> (parameterize ([array-strictness #f])
> - (array-all-sum (array-magnitude a))))
> + (array-all-sum (inline-array-map nonstupid-magnitude M))))
>
> -(: matrix-2norm ((Matrix Number) -> Nonnegative-Real))
> -(define (matrix-2norm a)
> +(: matrix-2norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Number) -> Nonnegative-Real)))
> +(define (matrix-2norm M)
> (parameterize ([array-strictness #f])
> - (let ([a (array-strict (array-magnitude a))])
> + (let ([M (array-strict (inline-array-map nonstupid-magnitude M))])
> ;; Compute this divided by the maximum to avoid underflow and overflow
> - (define mx (array-all-max a))
> + (define mx (array-all-max M))
> (cond [(and (rational? mx) (positive? mx))
> - (* mx (sqrt (array-all-sum
> - (inline-array-map (λ: ([x : Nonnegative-Real]) (sqr (/ x mx))) a))))]
> + (* mx (sqrt (array-all-sum (inline-array-map (λ (x) (sqr (/ x mx))) M))))]
> [else mx]))))
>
> -(: matrix-inf-norm ((Matrix Number) -> Nonnegative-Real))
> -(define (matrix-inf-norm a)
> +(: matrix-inf-norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Number) -> Nonnegative-Real)))
> +(define (matrix-inf-norm M)
> (parameterize ([array-strictness #f])
> - (array-all-max (array-magnitude a))))
> + (array-all-max (inline-array-map nonstupid-magnitude M))))
>
> -(: matrix-p-norm ((Matrix Number) Positive-Real -> Nonnegative-Real))
> -(define (matrix-p-norm a p)
> +(: matrix-p-norm (case-> ((Matrix Flonum) Positive-Real -> Nonnegative-Flonum)
> + ((Matrix Real) Positive-Real -> Nonnegative-Real)
> + ((Matrix Float-Complex) Positive-Real -> Nonnegative-Flonum)
> + ((Matrix Number) Positive-Real -> Nonnegative-Real)))
> +(define (matrix-p-norm M p)
> (parameterize ([array-strictness #f])
> - (let ([a (array-strict (array-magnitude a))])
> + (let ([M (array-strict (inline-array-map nonstupid-magnitude M))])
> ;; Compute this divided by the maximum to avoid underflow and overflow
> - (define mx (array-all-max a))
> + (define mx (array-all-max M))
> (cond [(and (rational? mx) (positive? mx))
> - (assert
> - (* mx (expt (array-all-sum
> - (inline-array-map (λ: ([x : Nonnegative-Real]) (expt (/ x mx) p)) a))
> - (/ p)))
> - (make-predicate Nonnegative-Real))]
> + (* mx (expt (array-all-sum (inline-array-map (λ (x) (expt (abs (/ x mx)) p)) M))
> + (/ p)))]
> [else mx]))))
>
> -(: matrix-norm (case-> ((Matrix Number) -> Nonnegative-Real)
> +(: matrix-norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Flonum) Real -> Nonnegative-Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> + ((Matrix Real) Real -> Nonnegative-Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Float-Complex) Real -> Nonnegative-Flonum)
> + ((Matrix Number) -> Nonnegative-Real)
> ((Matrix Number) Real -> Nonnegative-Real)))
> ;; Computes the p norm of a matrix
> (define (matrix-norm a [p 2])
> @@ -230,8 +259,12 @@
> [(p . > . 1) (matrix-p-norm a p)]
> [else (raise-argument-error 'matrix-norm "Real >= 1" 1 a p)]))
>
> -(: matrix-dot (case-> ((Matrix Real) -> Nonnegative-Real)
> +(: matrix-dot (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Flonum) (Matrix Flonum) -> Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> ((Matrix Real) (Matrix Real) -> Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) -> Nonnegative-Real)
> ((Matrix Number) (Matrix Number) -> Number)))
> ;; Computes the Frobenius inner product of a matrix with itself or of two matrices
> @@ -256,20 +289,30 @@
> (λ: ([js : Indexes])
> (* (aproc js) (conjugate (bproc js)))))))]))
>
> -(: matrix-cos-angle (case-> ((Matrix Real) (Matrix Real) -> Real)
> +(: matrix-cos-angle (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum)
> + ((Matrix Real) (Matrix Real) -> Real)
> + ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) (Matrix Number) -> Number)))
> (define (matrix-cos-angle M N)
> (/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N))))
>
> -(: matrix-angle (case-> ((Matrix Real) (Matrix Real) -> Real)
> +(: matrix-angle (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum)
> + ((Matrix Real) (Matrix Real) -> Real)
> + ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) (Matrix Number) -> Number)))
> (define (matrix-angle M N)
> (acos (matrix-cos-angle M N)))
>
> (: matrix-normalize
> - (All (A) (case-> ((Matrix Real) -> (Matrix Real))
> + (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Flonum) Real -> (Matrix Flonum))
> + ((Matrix Flonum) Real (-> A) -> (U A (Matrix Flonum)))
> + ((Matrix Real) -> (Matrix Real))
> ((Matrix Real) Real -> (Matrix Real))
> ((Matrix Real) Real (-> A) -> (U A (Matrix Real)))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Real -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Real (-> A) -> (U A (Matrix Float-Complex)))
> ((Matrix Number) -> (Matrix Number))
> ((Matrix Number) Real -> (Matrix Number))
> ((Matrix Number) Real (-> A) -> (U A (Matrix Number))))))
> @@ -291,19 +334,25 @@
> (define (matrix-transpose a)
> (array-axis-swap (ensure-matrix 'matrix-transpose a) 0 1))
>
> -(: matrix-conjugate (case-> ((Matrix Real) -> (Matrix Real))
> +(: matrix-conjugate (case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Real) -> (Matrix Real))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> ((Matrix Number) -> (Matrix Number))))
> (define (matrix-conjugate a)
> (array-conjugate (ensure-matrix 'matrix-conjugate a)))
>
> -(: matrix-hermitian (case-> ((Matrix Real) -> (Matrix Real))
> +(: matrix-hermitian (case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Real) -> (Matrix Real))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> ((Matrix Number) -> (Matrix Number))))
> (define (matrix-hermitian a)
> (array-default-strict
> (parameterize ([array-strictness #f])
> (array-axis-swap (array-conjugate (ensure-matrix 'matrix-hermitian a)) 0 1))))
>
> -(: matrix-trace (case-> ((Matrix Real) -> Real)
> +(: matrix-trace (case-> ((Matrix Flonum) -> Flonum)
> + ((Matrix Real) -> Real)
> + ((Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) -> Number)))
> (define (matrix-trace a)
> (cond [(square-matrix? a)
> @@ -349,15 +398,23 @@
> [else (fail)])]))]
> [else (fail)])]))
>
> -(: make-matrix-normalize (Real -> (case-> ((Matrix Real) -> (U #f (Matrix Real)))
> +(: make-matrix-normalize (Real -> (case-> ((Matrix Flonum) -> (U #f (Matrix Flonum)))
> + ((Matrix Real) -> (U #f (Matrix Real)))
> + ((Matrix Float-Complex) -> (U #f (Matrix Float-Complex)))
> ((Matrix Number) -> (U #f (Matrix Number))))))
> (define ((make-matrix-normalize p) M)
> (matrix-normalize M p (λ () #f)))
>
> (: matrix-normalize-rows
> - (All (A) (case-> ((Matrix Real) -> (Matrix Real))
> + (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Flonum) Real -> (Matrix Flonum))
> + ((Matrix Flonum) Real (-> A) -> (U A (Matrix Flonum)))
> + ((Matrix Real) -> (Matrix Real))
> ((Matrix Real) Real -> (Matrix Real))
> ((Matrix Real) Real (-> A) -> (U A (Matrix Real)))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Real -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Real (-> A) -> (U A (Matrix Float-Complex)))
> ((Matrix Number) -> (Matrix Number))
> ((Matrix Number) Real -> (Matrix Number))
> ((Matrix Number) Real (-> A) -> (U A (Matrix Number))))))
> @@ -371,9 +428,15 @@
> (matrix-map-rows (make-matrix-normalize p) M fail)]))
>
> (: matrix-normalize-cols
> - (All (A) (case-> ((Matrix Real) -> (Matrix Real))
> + (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Flonum) Real -> (Matrix Flonum))
> + ((Matrix Flonum) Real (-> A) -> (U A (Matrix Flonum)))
> + ((Matrix Real) -> (Matrix Real))
> ((Matrix Real) Real -> (Matrix Real))
> ((Matrix Real) Real (-> A) -> (U A (Matrix Real)))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Real -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Real (-> A) -> (U A (Matrix Float-Complex)))
> ((Matrix Number) -> (Matrix Number))
> ((Matrix Number) Real -> (Matrix Number))
> ((Matrix Number) Real (-> A) -> (U A (Matrix Number))))))
>
> collects/math/private/matrix/matrix-constructors.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-constructors.rkt
> +++ NEW/collects/math/private/matrix/matrix-constructors.rkt
> @@ -1,8 +1,10 @@
> #lang typed/racket/base
>
> (require racket/fixnum
> + racket/flonum
> racket/list
> racket/vector
> + math/base
> "matrix-types.rkt"
> "../unsafe.rkt"
> "../array/array-struct.rkt"
> @@ -22,8 +24,14 @@
> ;; ===================================================================================================
> ;; Basic constructors
>
> -(: identity-matrix (Integer -> (Matrix (U 0 1))))
> -(define (identity-matrix m) (diagonal-array 2 m 1 0))
> +(: identity-matrix (All (A) (case-> (Integer -> (Matrix (U 1 0)))
> + (Integer A -> (Matrix (U A 0)))
> + (Integer A A -> (Matrix A)))))
> +(define identity-matrix
> + (case-lambda
> + [(m) (diagonal-array 2 m 1 0)]
> + [(m one) (diagonal-array 2 m one 0)]
> + [(m one zero) (diagonal-array 2 m one zero)]))
>
> (: make-matrix (All (A) (Integer Integer A -> (Matrix A))))
> (define (make-matrix m n x)
> @@ -60,9 +68,12 @@
> (cond [(= i (unsafe-vector-ref js 1)) (unsafe-vector-ref vs i)]
> [else zero])))]))
>
> -(: diagonal-matrix (All (A) ((Listof A) -> (Matrix (U A 0)))))
> -(define (diagonal-matrix xs)
> - (diagonal-matrix/zero xs 0))
> +(: diagonal-matrix (All (A) (case-> ((Listof A) -> (Matrix (U A 0)))
> + ((Listof A) A -> (Matrix A)))))
> +(define diagonal-matrix
> + (case-lambda
> + [(xs) (diagonal-matrix/zero xs 0)]
> + [(xs zero) (diagonal-matrix/zero xs zero)]))
>
> ;; ===================================================================================================
> ;; Block diagonal matrices
> @@ -129,21 +140,29 @@
> [else
> (block-diagonal-matrix/zero* as zero)])))
>
> -(: block-diagonal-matrix (All (A) ((Listof (Matrix A)) -> (Matrix (U A 0)))))
> -(define (block-diagonal-matrix as)
> - (block-diagonal-matrix/zero as 0))
> +(: block-diagonal-matrix (All (A) (case-> ((Listof (Matrix A)) -> (Matrix (U A 0)))
> + ((Listof (Matrix A)) A -> (Matrix A)))))
> +(define block-diagonal-matrix
> + (case-lambda
> + [(as) (block-diagonal-matrix/zero as 0)]
> + [(as zero) (block-diagonal-matrix/zero as zero)]))
>
> ;; ===================================================================================================
> ;; Special matrices
>
> -(: expt-hack (case-> (Real Integer -> Real)
> - (Number Integer -> Number)))
> -;; Stop using this when TR correctly derives expt : Real Integer -> Real
> -(define (expt-hack x n)
> - (cond [(real? x) (assert (expt x n) real?)]
> +(: sane-expt (case-> (Flonum Index -> Flonum)
> + (Real Index -> Real)
> + (Float-Complex Index -> Float-Complex)
> + (Number Index -> Number)))
> +(define (sane-expt x n)
> + (cond [(flonum? x) (flexpt x (real->double-flonum n))]
> + [(real? x) (real-part (expt x n))] ; remove `real-part' when expt : Real Index -> Real
> + [(float-complex? x) (number->float-complex (expt x n))]
> [else (expt x n)]))
>
> -(: vandermonde-matrix (case-> ((Listof Real) Integer -> (Matrix Real))
> +(: vandermonde-matrix (case-> ((Listof Flonum) Integer -> (Matrix Flonum))
> + ((Listof Real) Integer -> (Matrix Real))
> + ((Listof Float-Complex) Integer -> (Matrix Float-Complex))
> ((Listof Number) Integer -> (Matrix Number))))
> (define (vandermonde-matrix xs n)
> (cond [(empty? xs)
> @@ -151,4 +170,4 @@
> [(or (not (index? n)) (zero? n))
> (raise-argument-error 'vandermonde-matrix "Positive-Index" 1 xs n)]
> [else
> - (array-axis-expand (list->array xs) 1 n expt-hack)]))
> + (array-axis-expand (list->array xs) 1 n sane-expt)]))
>
> collects/math/private/matrix/matrix-expt.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-expt.rkt
> +++ NEW/collects/math/private/matrix/matrix-expt.rkt
> @@ -3,11 +3,16 @@
> (require "matrix-types.rkt"
> "matrix-constructors.rkt"
> "matrix-arithmetic.rkt"
> - "utils.rkt")
> + "utils.rkt"
> + "../array/array-struct.rkt"
> + "../array/utils.rkt"
> + "../unsafe.rkt")
>
> (provide matrix-expt)
>
> -(: matrix-expt/ns (case-> ((Matrix Real) Positive-Integer -> (Matrix Real))
> +(: matrix-expt/ns (case-> ((Matrix Flonum) Positive-Integer -> (Matrix Flonum))
> + ((Matrix Real) Positive-Integer -> (Matrix Real))
> + ((Matrix Float-Complex) Positive-Integer -> (Matrix Float-Complex))
> ((Matrix Number) Positive-Integer -> (Matrix Number))))
> (define (matrix-expt/ns a n)
> (define n/2 (quotient n 2))
> @@ -22,10 +27,19 @@
> ;; m = n - 1
> (matrix* a (matrix-expt/ns a m))))))
>
> -(: matrix-expt (case-> ((Matrix Real) Integer -> (Matrix Real))
> +(: matrix-expt (case-> ((Matrix Flonum) Integer -> (Matrix Flonum))
> + ((Matrix Real) Integer -> (Matrix Real))
> + ((Matrix Float-Complex) Integer -> (Matrix Float-Complex))
> ((Matrix Number) Integer -> (Matrix Number))))
> (define (matrix-expt a n)
> (cond [(not (square-matrix? a)) (raise-argument-error 'matrix-expt "square-matrix?" 0 a n)]
> [(negative? n) (raise-argument-error 'matrix-expt "Natural" 1 a n)]
> - [(zero? n) (identity-matrix (square-matrix-size a))]
> + [(zero? n)
> + (define proc (unsafe-array-proc a))
> + (array-default-strict
> + (unsafe-build-array
> + (array-shape a)
> + (λ: ([ij : Indexes])
> + (define x (proc ij))
> + (if (= (unsafe-vector-ref ij 0) (unsafe-vector-ref ij 1)) (one* x) (zero* x)))))]
> [else (call/ns (λ () (matrix-expt/ns a n)))]))
>
> collects/math/private/matrix/matrix-gauss-elim.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-gauss-elim.rkt
> +++ NEW/collects/math/private/matrix/matrix-gauss-elim.rkt
> @@ -16,10 +16,18 @@
> (define-type Pivoting (U 'first 'partial))
>
> (: matrix-gauss-elim
> - (case-> ((Matrix Real) -> (Values (Matrix Real) (Listof Index)))
> + (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Listof Index)))
> + ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Listof Index)))
> + ((Matrix Flonum) Any Any -> (Values (Matrix Flonum) (Listof Index)))
> + ((Matrix Flonum) Any Any Pivoting -> (Values (Matrix Flonum) (Listof Index)))
> + ((Matrix Real) -> (Values (Matrix Real) (Listof Index)))
> ((Matrix Real) Any -> (Values (Matrix Real) (Listof Index)))
> ((Matrix Real) Any Any -> (Values (Matrix Real) (Listof Index)))
> ((Matrix Real) Any Any Pivoting -> (Values (Matrix Real) (Listof Index)))
> + ((Matrix Float-Complex) -> (Values (Matrix Float-Complex) (Listof Index)))
> + ((Matrix Float-Complex) Any -> (Values (Matrix Float-Complex) (Listof Index)))
> + ((Matrix Float-Complex) Any Any -> (Values (Matrix Float-Complex) (Listof Index)))
> + ((Matrix Float-Complex) Any Any Pivoting -> (Values (Matrix Float-Complex) (Listof Index)))
> ((Matrix Number) -> (Values (Matrix Number) (Listof Index)))
> ((Matrix Number) Any -> (Values (Matrix Number) (Listof Index)))
> ((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index)))
> @@ -52,17 +60,25 @@
> (vector-swap! rows i p)
> ;; Possibly unitize the new current row
> (let ([pivot (if unitize-pivot?
> - (begin (vector-scale! (unsafe-vector-ref rows i) (/ pivot))
> - 1)
> + (begin (vector-scale! (unsafe-vector-ref rows i) (/ 1 pivot))
> + (/ pivot pivot))
> pivot)])
> (elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
> (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
>
> (: matrix-row-echelon
> - (case-> ((Matrix Real) -> (Matrix Real))
> + (case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Flonum) Any -> (Matrix Flonum))
> + ((Matrix Flonum) Any Any -> (Matrix Flonum))
> + ((Matrix Flonum) Any Any Pivoting -> (Matrix Flonum))
> + ((Matrix Real) -> (Matrix Real))
> ((Matrix Real) Any -> (Matrix Real))
> ((Matrix Real) Any Any -> (Matrix Real))
> ((Matrix Real) Any Any Pivoting -> (Matrix Real))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Any -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Any Any -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) Any Any Pivoting -> (Matrix Float-Complex))
> ((Matrix Number) -> (Matrix Number))
> ((Matrix Number) Any -> (Matrix Number))
> ((Matrix Number) Any Any -> (Matrix Number))
>
> collects/math/private/matrix/matrix-gram-schmidt.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-gram-schmidt.rkt
> +++ NEW/collects/math/private/matrix/matrix-gram-schmidt.rkt
> @@ -16,7 +16,9 @@
> (provide matrix-gram-schmidt
> matrix-basis-extension)
>
> -(: find-nonzero-vector (case-> ((Vectorof (Vectorof Real)) -> (U #f Index))
> +(: find-nonzero-vector (case-> ((Vectorof (Vectorof Flonum)) -> (U #f Index))
> + ((Vectorof (Vectorof Real)) -> (U #f Index))
> + ((Vectorof (Vectorof Float-Complex)) -> (U #f Index))
> ((Vectorof (Vectorof Number)) -> (U #f Index))))
> (define (find-nonzero-vector vss)
> (define n (vector-length vss))
> @@ -28,7 +30,10 @@
> [else #f]))]))
>
> (: subtract-projections!
> - (case-> ((Vectorof (Vectorof Real)) Nonnegative-Fixnum Index (Vectorof Real) -> Void)
> + (case-> ((Vectorof (Vectorof Flonum)) Nonnegative-Fixnum Index (Vectorof Flonum) -> Void)
> + ((Vectorof (Vectorof Real)) Nonnegative-Fixnum Index (Vectorof Real) -> Void)
> + ((Vectorof (Vectorof Float-Complex)) Nonnegative-Fixnum Index (Vectorof Float-Complex)
> + -> Void)
> ((Vectorof (Vectorof Number)) Nonnegative-Fixnum Index (Vectorof Number) -> Void)))
> (define (subtract-projections! rows i m row)
> (let loop ([#{i : Nonnegative-Fixnum} i])
> @@ -36,7 +41,9 @@
> (vector-sub-proj! (unsafe-vector-ref rows i) row #f)
> (loop (fx+ i 1)))))
>
> -(: matrix-gram-schmidt/ns (case-> ((Matrix Real) Any Integer -> (Array Real))
> +(: matrix-gram-schmidt/ns (case-> ((Matrix Flonum) Any Integer -> (Array Flonum))
> + ((Matrix Real) Any Integer -> (Array Real))
> + ((Matrix Float-Complex) Any Integer -> (Array Float-Complex))
> ((Matrix Number) Any Integer -> (Array Number))))
> ;; Performs Gram-Schmidt orthogonalization on M, assuming the rows before `start' are already
> ;; orthogonal
> @@ -60,31 +67,46 @@
> [else
> (matrix-transpose (vector*->matrix (list->vector (reverse bs))))]))]
> [else
> - (make-array (vector (matrix-num-rows M) 0) 0)]))
> + (make-array (vector (matrix-num-rows M) 0)
> + ;; Value won't be in the matrix, but this satisfies TR:
> + (zero* (unsafe-vector2d-ref rows 0 0)))]))
>
> -(: matrix-gram-schmidt (case-> ((Matrix Real) -> (Array Real))
> +(: matrix-gram-schmidt (case-> ((Matrix Flonum) -> (Array Flonum))
> + ((Matrix Flonum) Any -> (Array Flonum))
> + ((Matrix Flonum) Any Integer -> (Array Flonum))
> + ((Matrix Real) -> (Array Real))
> ((Matrix Real) Any -> (Array Real))
> ((Matrix Real) Any Integer -> (Array Real))
> + ((Matrix Float-Complex) -> (Array Float-Complex))
> + ((Matrix Float-Complex) Any -> (Array Float-Complex))
> + ((Matrix Float-Complex) Any Integer -> (Array Float-Complex))
> ((Matrix Number) -> (Array Number))
> ((Matrix Number) Any -> (Array Number))
> ((Matrix Number) Any Integer -> (Array Number))))
> (define (matrix-gram-schmidt M [normalize? #f] [start 0])
> (call/ns (λ () (matrix-gram-schmidt/ns M normalize? start))))
>
> -(: matrix-basis-extension/ns (case-> ((Matrix Real) -> (Array Real))
> +(: matrix-basis-extension/ns (case-> ((Matrix Flonum) -> (Array Flonum))
> + ((Matrix Real) -> (Array Real))
> + ((Matrix Float-Complex) -> (Array Float-Complex))
> ((Matrix Number) -> (Array Number))))
> (define (matrix-basis-extension/ns B)
> (define-values (m n) (matrix-shape B))
> + (define x00 (matrix-ref B 0 0))
> + (define zero (zero* x00))
> + (define one (one* x00))
> (cond [(n . < . m)
> - (define S (matrix-gram-schmidt (matrix-augment (list B (identity-matrix m))) #f n))
> + (define S (matrix-gram-schmidt (matrix-augment (list B (identity-matrix m one zero))) #f n))
> (define R (submatrix S (::) (:: n #f)))
> (matrix-augment (take (sort/key (matrix-cols R) > matrix-norm) (- m n)))]
> [(n . = . m)
> - (make-array (vector m 0) 0)]
> + (make-array (vector m 0) zero)]
> [else
> (raise-argument-error 'matrix-extend-row-basis "matrix? with width < height" B)]))
>
> -(: matrix-basis-extension (case-> ((Matrix Real) -> (Array Real))
> +(: matrix-basis-extension (case-> ((Matrix Flonum) -> (Array Flonum))
> + ((Matrix Real) -> (Array Real))
> + ((Matrix Float-Complex) -> (Array Float-Complex))
> ((Matrix Number) -> (Array Number))))
> (define (matrix-basis-extension B)
> (call/ns (λ () (matrix-basis-extension/ns B))))
>
> collects/math/private/matrix/matrix-lu.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-lu.rkt
> +++ NEW/collects/math/private/matrix/matrix-lu.rkt
> @@ -8,15 +8,22 @@
> "../unsafe.rkt"
> "../vector/vector-mutate.rkt"
> "../array/mutable-array.rkt"
> - "../array/array-struct.rkt")
> + "../array/array-struct.rkt"
> + "../array/array-pointwise.rkt")
>
> (provide matrix-lu)
>
> ;; An LU factorization exists iff Gaussian elimination can be done without row swaps.
>
> (: matrix-lu
> - (All (A) (case-> ((Matrix Real) -> (Values (Matrix Real) (Matrix Real)))
> + (All (A) (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Matrix Flonum)))
> + ((Matrix Flonum) (-> A) -> (Values (U A (Matrix Flonum)) (Matrix Flonum)))
> + ((Matrix Real) -> (Values (Matrix Real) (Matrix Real)))
> ((Matrix Real) (-> A) -> (Values (U A (Matrix Real)) (Matrix Real)))
> + ((Matrix Float-Complex) -> (Values (Matrix Float-Complex)
> + (Matrix Float-Complex)))
> + ((Matrix Float-Complex) (-> A) -> (Values (U A (Matrix Float-Complex))
> + (Matrix Float-Complex)))
> ((Matrix Number) -> (Values (Matrix Number) (Matrix Number)))
> ((Matrix Number) (-> A) -> (Values (U A (Matrix Number)) (Matrix Number))))))
> (define matrix-lu
> @@ -28,7 +35,7 @@
> (define L
> (parameterize ([array-strictness #f])
> ;; Construct L in a weird way to prove to TR that it has the right type
> - (array->mutable-array (matrix-scale M (ann 0 Real)))))
> + (array->mutable-array (inline-array-map zero* M))))
> ;; Going to fill in the lower triangle by banging values into `ys'
> (define ys (mutable-array-data L))
> (let loop ([#{i : Nonnegative-Fixnum} 0])
> @@ -51,12 +58,13 @@
> ;; Add row i, scaled
> (vector-scaled-add! (unsafe-vector-ref rows l)
> (unsafe-vector-ref rows i)
> - (- y_li)))
> + (* -1 y_li)))
> (l-loop (fx+ l 1))]
> [else
> (loop (fx+ i 1))]))])]
> [else
> ;; L's lower triangle has been filled; now fill the diagonal with 1s
> (for: ([i : Integer (in-range 0 m)])
> - (vector-set! ys (+ (* i m) i) 1))
> + (define j (+ (* i m) i))
> + (vector-set! ys j (one* (vector-ref ys j))))
> (values L (vector*->matrix rows))]))]))
>
> collects/math/private/matrix/matrix-operator-norm.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-operator-norm.rkt
> +++ NEW/collects/math/private/matrix/matrix-operator-norm.rkt
> @@ -42,33 +42,46 @@ See "How to Measure Errors" in the LAPACK manual for more details:
> matrix-orthonormal?
> )
>
> -(: matrix-op-1norm ((Matrix Number) -> Nonnegative-Real))
> +(: matrix-op-1norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Number) -> Nonnegative-Real)))
> ;; When M is a column matrix, this is equivalent to matrix-1norm
> (define (matrix-op-1norm M)
> (parameterize ([array-strictness #f])
> (assert (apply max (map matrix-1norm (matrix-cols M))) nonnegative?)))
>
> -(: matrix-op-2norm ((Matrix Number) -> Nonnegative-Real))
> +(: matrix-op-2norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Number) -> Nonnegative-Real)))
> ;; When M is a column matrix, this is equivalent to matrix-2norm
> (define (matrix-op-2norm M)
> ;(matrix-max-singular-value M)
> ;(sqrt (matrix-max-eigenvalue M))
> (error 'unimplemented))
>
> -(: matrix-op-inf-norm ((Matrix Number) -> Nonnegative-Real))
> +(: matrix-op-inf-norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
> + ((Matrix Real) -> Nonnegative-Real)
> + ((Matrix Float-Complex) -> Nonnegative-Flonum)
> + ((Matrix Number) -> Nonnegative-Real)))
> ;; When M is a column matrix, this is equivalent to matrix-inf-norm
> (define (matrix-op-inf-norm M)
> (parameterize ([array-strictness #f])
> (assert (apply max (map matrix-1norm (matrix-rows M))) nonnegative?)))
>
> -(: matrix-basis-cos-angle (case-> ((Matrix Real) (Matrix Real) -> Real)
> +(: matrix-basis-cos-angle (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum)
> + ((Matrix Real) (Matrix Real) -> Real)
> + ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) (Matrix Number) -> Number)))
> ;; Returns the angle between the two subspaces spanned by the two given sets of column vectors
> (define (matrix-basis-cos-angle M R)
> ;(matrix-min-singular-value (matrix* (matrix-hermitian M) R))
> (error 'unimplemented))
>
> -(: matrix-basis-angle (case-> ((Matrix Real) (Matrix Real) -> Real)
> +(: matrix-basis-angle (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum)
> + ((Matrix Real) (Matrix Real) -> Real)
> + ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) (Matrix Number) -> Number)))
> ;; Returns the angle between the two subspaces spanned by the two given sets of column vectors
> (define (matrix-basis-angle M R)
>
> collects/math/private/matrix/matrix-qr.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-qr.rkt
> +++ NEW/collects/math/private/matrix/matrix-qr.rkt
> @@ -5,6 +5,7 @@
> "matrix-arithmetic.rkt"
> "matrix-constructors.rkt"
> "matrix-gram-schmidt.rkt"
> + "utils.rkt"
> "../array/array-transform.rkt"
> "../array/array-struct.rkt")
>
> @@ -24,21 +25,33 @@ produces matrices for which `matrix-orthogonal?' returns #t with eps <= 10*epsil
> independently of the matrix size.
> |#
>
> -(: matrix-qr/ns (case-> ((Matrix Real) Any -> (Values (Matrix Real) (Matrix Real)))
> +(: matrix-qr/ns (case-> ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Matrix Flonum)))
> + ((Matrix Real) Any -> (Values (Matrix Real) (Matrix Real)))
> + ((Matrix Float-Complex) Any -> (Values (Matrix Float-Complex)
> + (Matrix Float-Complex)))
> ((Matrix Number) Any -> (Values (Matrix Number) (Matrix Number)))))
> (define (matrix-qr/ns M full?)
> + (define x00 (matrix-ref M 0 0))
> + (define zero (zero* x00))
> + (define one (one* x00))
> (define B (matrix-gram-schmidt M #f))
> (define Q
> (matrix-gram-schmidt
> (cond [(or (square-matrix? B) (and (matrix? B) (not full?))) B]
> [(matrix? B) (array-append* (list B (matrix-basis-extension B)) 1)]
> - [full? (identity-matrix (matrix-num-rows M))]
> - [else (matrix-col (identity-matrix (matrix-num-rows M)) 0)])
> + [full? (identity-matrix (matrix-num-rows M) one zero)]
> + [else (matrix-col (identity-matrix (matrix-num-rows M) one zero) 0)])
> #t))
> - (values Q (matrix-upper-triangle (matrix* (matrix-hermitian Q) M))))
> + (values Q (matrix-upper-triangle (matrix* (matrix-hermitian Q) M) zero)))
>
> -(: matrix-qr (case-> ((Matrix Real) -> (Values (Matrix Real) (Matrix Real)))
> +(: matrix-qr (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Matrix Flonum)))
> + ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Matrix Flonum)))
> + ((Matrix Real) -> (Values (Matrix Real) (Matrix Real)))
> ((Matrix Real) Any -> (Values (Matrix Real) (Matrix Real)))
> + ((Matrix Float-Complex) -> (Values (Matrix Float-Complex)
> + (Matrix Float-Complex)))
> + ((Matrix Float-Complex) Any -> (Values (Matrix Float-Complex)
> + (Matrix Float-Complex)))
> ((Matrix Number) -> (Values (Matrix Number) (Matrix Number)))
> ((Matrix Number) Any -> (Values (Matrix Number) (Matrix Number)))))
> (define (matrix-qr M [full? #t])
>
> collects/math/private/matrix/matrix-solve.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-solve.rkt
> +++ NEW/collects/math/private/matrix/matrix-solve.rkt
> @@ -24,7 +24,9 @@
> ;; ===================================================================================================
> ;; Determinant
>
> -(: matrix-determinant (case-> ((Matrix Real) -> Real)
> +(: matrix-determinant (case-> ((Matrix Flonum) -> Flonum)
> + ((Matrix Real) -> Real)
> + ((Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) -> Number)))
> (define (matrix-determinant M)
> (define m (square-matrix-size M))
> @@ -41,20 +43,22 @@
> [else
> (matrix-determinant/row-reduction M)]))
>
> -(: matrix-determinant/row-reduction (case-> ((Matrix Real) -> Real)
> +(: matrix-determinant/row-reduction (case-> ((Matrix Flonum) -> Flonum)
> + ((Matrix Real) -> Real)
> + ((Matrix Float-Complex) -> Float-Complex)
> ((Matrix Number) -> Number)))
> (define (matrix-determinant/row-reduction M)
> (define m (square-matrix-size M))
> (define rows (matrix->vector* M))
> - (let loop ([#{i : Nonnegative-Fixnum} 0] [#{sign : Real} 1])
> + (let loop ([#{i : Nonnegative-Fixnum} 0] [#{sign : (U Positive-Fixnum Negative-Fixnum)} 1])
> (cond
> [(i . fx< . m)
> (define-values (p pivot) (find-partial-pivot rows m i i))
> (cond
> - [(zero? pivot) 0] ; no pivot means non-invertible matrix
> + [(zero? pivot) pivot] ; no pivot means non-invertible matrix
> [else
> (let ([sign (if (= i p) sign (begin (vector-swap! rows i p) ; swapping negates sign
> - (* -1 sign)))])
> + (if (= sign 1) -1 1)))])
> (elim-rows! rows m i i pivot (fx+ i 1)) ; adding scaled rows doesn't change it
> (loop (fx+ i 1) sign))])]
> [else
> @@ -72,8 +76,12 @@
> (and (square-matrix? M)
> (not (zero? (matrix-determinant M)))))
>
> -(: matrix-inverse (All (A) (case-> ((Matrix Real) -> (Matrix Real))
> +(: matrix-inverse (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Flonum) (-> A) -> (U A (Matrix Flonum)))
> + ((Matrix Real) -> (Matrix Real))
> ((Matrix Real) (-> A) -> (U A (Matrix Real)))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) (-> A) -> (U A (Matrix Float-Complex)))
> ((Matrix Number) -> (Matrix Number))
> ((Matrix Number) (-> A) -> (U A (Matrix Number))))))
> (define matrix-inverse
> @@ -81,7 +89,8 @@
> [(M) (matrix-inverse M (λ () (raise-argument-error 'matrix-inverse "matrix-invertible?" M)))]
> [(M fail)
> (define m (square-matrix-size M))
> - (define I (identity-matrix m))
> + (define x00 (matrix-ref M 0 0))
> + (define I (identity-matrix m (one* x00) (zero* x00)))
> (define-values (IM^-1 wps) (parameterize ([array-strictness #f])
> (matrix-gauss-elim (matrix-augment (list M I)) #t #t)))
> (cond [(and (not (empty? wps)) (= (first wps) m))
> @@ -91,11 +100,16 @@
> ;; ===================================================================================================
> ;; Solving linear systems
>
> -(: matrix-solve (All (A) (case->
> - ((Matrix Real) (Matrix Real) -> (Matrix Real))
> - ((Matrix Real) (Matrix Real) (-> A) -> (U A (Matrix Real)))
> - ((Matrix Number) (Matrix Number) -> (Matrix Number))
> - ((Matrix Number) (Matrix Number) (-> A) -> (U A (Matrix Number))))))
> +(: matrix-solve
> + (All (A) (case->
> + ((Matrix Flonum) (Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Flonum) (Matrix Flonum) (-> A) -> (U A (Matrix Flonum)))
> + ((Matrix Real) (Matrix Real) -> (Matrix Real))
> + ((Matrix Real) (Matrix Real) (-> A) -> (U A (Matrix Real)))
> + ((Matrix Float-Complex) (Matrix Float-Complex) -> (Matrix Float-Complex))
> + ((Matrix Float-Complex) (Matrix Float-Complex) (-> A) -> (U A (Matrix Float-Complex)))
> + ((Matrix Number) (Matrix Number) -> (Matrix Number))
> + ((Matrix Number) (Matrix Number) (-> A) -> (U A (Matrix Number))))))
> (define matrix-solve
> (case-lambda
> [(M B) (matrix-solve M B (λ () (raise-argument-error 'matrix-solve "matrix-invertible?" 0 M B)))]
>
> collects/math/private/matrix/matrix-subspace.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/matrix-subspace.rkt
> +++ NEW/collects/math/private/matrix/matrix-subspace.rkt
> @@ -35,7 +35,9 @@
> (cond [(= j0 j1) Bs]
> [else (cons (submatrix M (::) (:: j0 j1)) Bs)]))
>
> -(: matrix-col-space/ns (All (A) (case-> ((Matrix Real) -> (U #f (Matrix Real)))
> +(: matrix-col-space/ns (All (A) (case-> ((Matrix Flonum) -> (U #f (Matrix Flonum)))
> + ((Matrix Real) -> (U #f (Matrix Real)))
> + ((Matrix Float-Complex) -> (U #f (Matrix Float-Complex)))
> ((Matrix Number) -> (U #f (Matrix Number))))))
> (define (matrix-col-space/ns M)
> (define n (matrix-num-cols M))
> @@ -52,13 +54,19 @@
> (define next-j (first wps))
> (loop next-j (rest wps) (maybe-cons-submatrix M (fx+ j 1) next-j Bs))]))]))
>
> -(: matrix-col-space (All (A) (case-> ((Matrix Real) -> (Matrix Real))
> +(: matrix-col-space (All (A) (case-> ((Matrix Flonum) -> (Array Flonum))
> + ((Matrix Flonum) (-> A) -> (U A (Matrix Flonum)))
> + ((Matrix Real) -> (Array Real))
> ((Matrix Real) (-> A) -> (U A (Matrix Real)))
> - ((Matrix Number) -> (Matrix Number))
> + ((Matrix Float-Complex) -> (Array Float-Complex))
> + ((Matrix Float-Complex) (-> A) -> (U A (Matrix Float-Complex)))
> + ((Matrix Number) -> (Array Number))
> ((Matrix Number) (-> A) -> (U A (Matrix Number))))))
> (define matrix-col-space
> (case-lambda
> - [(M) (matrix-col-space M (λ () (make-array (vector 0 (matrix-num-cols M)) 0)))]
> + [(M) (matrix-col-space M (λ ()
> + (define x00 (matrix-ref M 0 0))
> + (make-array (vector 0 (matrix-num-cols M)) (zero* x00))))]
> [(M fail)
> (define S (parameterize ([array-strictness #f])
> (matrix-col-space/ns M)))
>
> collects/math/private/matrix/typed-matrix-arithmetic.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/typed-matrix-arithmetic.rkt
> +++ NEW/collects/math/private/matrix/typed-matrix-arithmetic.rkt
> @@ -65,46 +65,65 @@
> (loop (first arrs) (rest arrs)))])))]))
>
>
> -(: matrix*/ns (case-> ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
> - ((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
> +(: matrix*/ns
> + (case-> ((Matrix Flonum) (Listof (Matrix Flonum)) -> (Matrix Flonum))
> + ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
> + ((Matrix Float-Complex) (Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
> + ((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
> (define (matrix*/ns a as)
> (cond [(empty? as) a]
> [else (matrix*/ns (inline-matrix-multiply a (first as)) (rest as))]))
>
> -(: matrix* (case-> ((Matrix Real) (Matrix Real) * -> (Matrix Real))
> +(: matrix* (case-> ((Matrix Flonum) (Matrix Flonum) * -> (Matrix Flonum))
> + ((Matrix Real) (Matrix Real) * -> (Matrix Real))
> + ((Matrix Float-Complex) (Matrix Float-Complex) * -> (Matrix Float-Complex))
> ((Matrix Number) (Matrix Number) * -> (Matrix Number))))
> (define (matrix* a . as) (call/ns (λ () (matrix*/ns a as))))
>
>
> -(: matrix+/ns (case-> ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
> - ((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
> +(: matrix+/ns
> + (case-> ((Matrix Flonum) (Listof (Matrix Flonum)) -> (Matrix Flonum))
> + ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
> + ((Matrix Float-Complex) (Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
> + ((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
> (define (matrix+/ns a as)
> (cond [(empty? as) a]
> [else (matrix+/ns (inline-matrix+ a (first as)) (rest as))]))
>
> -(: matrix+ (case-> ((Matrix Real) (Matrix Real) * -> (Matrix Real))
> +(: matrix+ (case-> ((Matrix Flonum) (Matrix Flonum) * -> (Matrix Flonum))
> + ((Matrix Real) (Matrix Real) * -> (Matrix Real))
> + ((Matrix Float-Complex) (Matrix Float-Complex) * -> (Matrix Float-Complex))
> ((Matrix Number) (Matrix Number) * -> (Matrix Number))))
> (define (matrix+ a . as) (call/ns (λ () (matrix+/ns a as))))
>
>
> -(: matrix-/ns (case-> ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
> - ((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
> +(: matrix-/ns
> + (case-> ((Matrix Flonum) (Listof (Matrix Flonum)) -> (Matrix Flonum))
> + ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
> + ((Matrix Float-Complex) (Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
> + ((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
> (define (matrix-/ns a as)
> (cond [(empty? as) a]
> [else (matrix-/ns (inline-matrix- a (first as)) (rest as))]))
>
> -(: matrix- (case-> ((Matrix Real) (Matrix Real) * -> (Matrix Real))
> +(: matrix- (case-> ((Matrix Flonum) (Matrix Flonum) * -> (Matrix Flonum))
> + ((Matrix Real) (Matrix Real) * -> (Matrix Real))
> + ((Matrix Float-Complex) (Matrix Float-Complex) * -> (Matrix Float-Complex))
> ((Matrix Number) (Matrix Number) * -> (Matrix Number))))
> (define (matrix- a . as)
> - (call/ns (λ () (cond [(empty? as) (inline-matrix- a)]
> + (call/ns (λ () (cond [(empty? as) (inline-matrix-scale a -1)]
> [else (matrix-/ns a as)]))))
>
>
> -(: matrix-scale (case-> ((Matrix Real) Real -> (Matrix Real))
> +(: matrix-scale (case-> ((Matrix Flonum) Flonum -> (Matrix Flonum))
> + ((Matrix Real) Real -> (Matrix Real))
> + ((Matrix Float-Complex) Float-Complex -> (Matrix Float-Complex))
> ((Matrix Number) Number -> (Matrix Number))))
> (define (matrix-scale a x) (inline-matrix-scale a x))
>
> -(: matrix-sum (case-> ((Listof (Matrix Real)) -> (Matrix Real))
> +(: matrix-sum (case-> ((Listof (Matrix Flonum)) -> (Matrix Flonum))
> + ((Listof (Matrix Real)) -> (Matrix Real))
> + ((Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
> ((Listof (Matrix Number)) -> (Matrix Number))))
> (define (matrix-sum lst)
> (cond [(empty? lst) (raise-argument-error 'matrix-sum "nonempty List" lst)]
>
> collects/math/private/matrix/untyped-matrix-arithmetic.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/untyped-matrix-arithmetic.rkt
> +++ NEW/collects/math/private/matrix/untyped-matrix-arithmetic.rkt
> @@ -8,39 +8,67 @@
> inline-matrix-map
> matrix-map)
>
> -(module syntax-defs racket/base
> - (require (for-syntax racket/base)
> - (only-in typed/racket/base λ: : inst Index)
> +(module typed-multiply-defs typed/racket/base
> + (require racket/fixnum
> "matrix-types.rkt"
> "utils.rkt"
> + "../unsafe.rkt"
> "../array/array-struct.rkt"
> - "../array/array-fold.rkt"
> "../array/array-transform.rkt"
> + "../array/mutable-array.rkt"
> "../array/utils.rkt")
>
> (provide (all-defined-out))
>
> - ;(: matrix-multiply ((Matrix Number) (Matrix Number) -> (Matrix Number)))
> + (: matrix-multiply-data (All (A) ((Matrix A) (Matrix A) -> (Values Index Index Index
> + (Vectorof A) (Vectorof A)
> + (-> (Boxof A))))))
> + (define (matrix-multiply-data arr brr)
> + (let-values ([(m p n) (matrix-multiply-shape arr brr)])
> + (define arr-data (mutable-array-data (array->mutable-array arr)))
> + (define brr-data (mutable-array-data (array->mutable-array (parameterize ([array-strictness #f])
> + (array-axis-swap brr 0 1)))))
> + (define bx (make-thread-local-box (unsafe-vector-ref arr-data 0)))
> + (values m p n arr-data brr-data bx)))
> +
> + (: make-matrix-multiply (All (A) (Index Index Index (Index Index -> A) -> (Matrix A))))
> + (define (make-matrix-multiply m p n sum-loop)
> + (array-default-strict
> + (unsafe-build-array
> + ((inst vector Index) m n)
> + (λ: ([ij : Indexes])
> + (sum-loop (assert (fx* (unsafe-vector-ref ij 0) p) index?)
> + (assert (fx* (unsafe-vector-ref ij 1) p) index?))))))
> + ) ; module
> +
> +(module untyped-multiply-defs racket/base
> + (require (for-syntax racket/base)
> + racket/fixnum
> + racket/unsafe/ops
> + (only-in typed/racket/base λ: : Index let: Nonnegative-Fixnum)
> + (submod ".." typed-multiply-defs)
> + "matrix-types.rkt"
> + "utils.rkt"
> + "../array/array-struct.rkt")
> +
> + (provide (all-defined-out))
> +
> ;; This is a macro so the result can have as precise a type as possible
> - (define-syntax-rule (inline-matrix-multiply arr-expr brr-expr)
> - (let ([arr arr-expr]
> - [brr brr-expr])
> - (let-values ([(m p n) (matrix-multiply-shape arr brr)]
> - ;; Make arr strict because its elements are reffed multiple times
> - [(_) (array-strict! arr)])
> - (let (;; Extend arr in the center dimension
> - [arr-proc (unsafe-array-proc (array-axis-insert arr 1 n))]
> - ;; Transpose brr and extend in the leftmost dimension
> - [brr-proc (unsafe-array-proc
> - (array-axis-insert (array-strict (array-axis-swap brr 0 1)) 0 m))])
> - ;; The *transpose* of brr is traversed in row-major order when this result is traversed
> - ;; in row-major order (which is why the transpose is strict, not brr)
> - (array-axis-sum
> - (unsafe-build-array
> - ((inst vector Index) m n p)
> - (λ: ([js : Indexes])
> - (* (arr-proc js) (brr-proc js))))
> - 2)))))
> + (define-syntax-rule (inline-matrix-multiply arr brr)
> + (let-values ([(m p n arr-data brr-data bx) (matrix-multiply-data arr brr)])
> + (make-matrix-multiply
> + m p n
> + (λ: ([i : Index] [j : Index])
> + (let ([bx (bx)]
> + [v (* (unsafe-vector-ref arr-data i)
> + (unsafe-vector-ref brr-data j))])
> + (let: loop ([k : Nonnegative-Fixnum 1] [v v])
> + (cond [(k . fx< . p)
> + (loop (fx+ k 1)
> + (+ v (* (unsafe-vector-ref arr-data (fx+ i k))
> + (unsafe-vector-ref brr-data (fx+ j k)))))]
> + [else (set-box! bx v)]))
> + (unbox bx))))))
>
> (define-syntax (do-inline-matrix* stx)
> (syntax-case stx ()
> @@ -54,6 +82,19 @@
> (parameterize ([array-strictness #f])
> (do-inline-matrix* arr brrs ...))))
>
> + ) ; module
> +
> +(module syntax-defs racket/base
> + (require (for-syntax racket/base)
> + (only-in typed/racket/base λ: : inst Index)
> + (submod ".." typed-multiply-defs)
> + "matrix-types.rkt"
> + "utils.rkt"
> + "../array/array-struct.rkt"
> + "../array/utils.rkt")
> +
> + (provide (all-defined-out))
> +
> (define-syntax (inline-matrix-map stx)
> (syntax-case stx ()
> [(_ f arr-expr)
> @@ -117,5 +158,6 @@
>
> ) ; module
>
> -(require 'syntax-defs
> +(require 'untyped-multiply-defs
> + 'syntax-defs
> 'untyped-defs)
>
> collects/math/private/matrix/utils.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/matrix/utils.rkt
> +++ NEW/collects/math/private/matrix/utils.rkt
> @@ -3,6 +3,7 @@
> (require racket/performance-hint
> racket/string
> racket/fixnum
> + math/base
> "matrix-types.rkt"
> "../unsafe.rkt"
> "../array/array-struct.rkt"
> @@ -69,7 +70,9 @@
> (rational? (imag-part z)))])))
>
> (: find-partial-pivot
> - (case-> ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
> + (case-> ((Vectorof (Vectorof Flonum)) Index Index Index -> (Values Index Flonum))
> + ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
> + ((Vectorof (Vectorof Float-Complex)) Index Index Index -> (Values Index Float-Complex))
> ((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number))))
> ;; Find the element with maximum magnitude in a column
> (define (find-partial-pivot rows m i j)
> @@ -85,7 +88,9 @@
> [else (values p pivot)])))
>
> (: find-first-pivot
> - (case-> ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
> + (case-> ((Vectorof (Vectorof Flonum)) Index Index Index -> (Values Index Flonum))
> + ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
> + ((Vectorof (Vectorof Float-Complex)) Index Index Index -> (Values Index Float-Complex))
> ((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number))))
> ;; Find the first nonzero element in a column
> (define (find-first-pivot rows m i j)
> @@ -100,7 +105,10 @@
> (values i pivot)]))))
>
> (: elim-rows!
> - (case-> ((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void)
> + (case-> ((Vectorof (Vectorof Flonum)) Index Index Index Flonum Nonnegative-Fixnum -> Void)
> + ((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void)
> + ((Vectorof (Vectorof Float-Complex)) Index Index Index Float-Complex Nonnegative-Fixnum
> + -> Void)
> ((Vectorof (Vectorof Number)) Index Index Index Number Nonnegative-Fixnum -> Void)))
> (define (elim-rows! rows m i j pivot start)
> (define row_i (unsafe-vector-ref rows i))
> @@ -109,8 +117,8 @@
> (unless (l . fx= . i)
> (define row_l (unsafe-vector-ref rows l))
> (define x_lj (unsafe-vector-ref row_l j))
> - (unless (zero? x_lj)
> - (vector-scaled-add! row_l row_i (- (/ x_lj pivot)) j)
> + (unless (= x_lj 0)
> + (vector-scaled-add! row_l row_i (* -1 (/ x_lj pivot)) j)
> ;; Make sure the element below the pivot is zero
> (unsafe-vector-set! row_l j (- x_lj x_lj))))
> (loop (fx+ l 1)))))
> @@ -124,3 +132,51 @@
> (thnk))))
>
> ) ; begin-encourage-inline
> +
> +(: make-thread-local-box (All (A) (A -> (-> (Boxof A)))))
> +(define (make-thread-local-box contents)
> + (let: ([val : (Thread-Cellof (U #f (Boxof A))) (make-thread-cell #f)])
> + (λ () (or (thread-cell-ref val)
> + (let: ([v : (Boxof A) (box contents)])
> + (thread-cell-set! val v)
> + v)))))
> +
> +(: one (case-> (Flonum -> Nonnegative-Flonum)
> + (Real -> (U 1 Nonnegative-Flonum))
> + (Float-Complex -> Nonnegative-Flonum)
> + (Number -> (U 1 Nonnegative-Flonum))))
> +(define (one x)
> + (cond [(flonum? x) 1.0]
> + [(real? x) 1]
> + [(float-complex? x) 1.0]
> + [else 1]))
> +
> +(: zero (case-> (Flonum -> Flonum-Positive-Zero)
> + (Real -> (U 0 Flonum-Positive-Zero))
> + (Float-Complex -> Flonum-Positive-Zero)
> + (Number -> (U 0 Flonum-Positive-Zero))))
> +(define (zero x)
> + (cond [(flonum? x) 0.0]
> + [(real? x) 0]
> + [(float-complex? x) 0.0]
> + [else 0]))
> +
> +(: one* (case-> (Flonum -> Nonnegative-Flonum)
> + (Real -> (U 1 Nonnegative-Flonum))
> + (Float-Complex -> Float-Complex)
> + (Number -> (U 1 Nonnegative-Flonum Float-Complex))))
> +(define (one* x)
> + (cond [(flonum? x) 1.0]
> + [(real? x) 1]
> + [(float-complex? x) 1.0+0.0i]
> + [else 1]))
> +
> +(: zero* (case-> (Flonum -> Flonum-Positive-Zero)
> + (Real -> (U 0 Flonum-Positive-Zero))
> + (Float-Complex -> Float-Complex)
> + (Number -> (U 0 Flonum-Positive-Zero Float-Complex))))
> +(define (zero* x)
> + (cond [(flonum? x) 0.0]
> + [(real? x) 0]
> + [(float-complex? x) 0.0+0.0i]
> + [else 0]))
>
> collects/math/private/vector/vector-mutate.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/private/vector/vector-mutate.rkt
> +++ NEW/collects/math/private/vector/vector-mutate.rkt
> @@ -1,7 +1,7 @@
> #lang typed/racket/base
>
> (require racket/fixnum
> - racket/math
> + math/base
> math/private/unsafe)
>
> (provide vector-swap!
> @@ -14,9 +14,12 @@
> vector-zero!
> vector-zero?)
>
> -(: mag^2 (Number -> Nonnegative-Real))
> +(: mag^2 (case-> (Flonum -> Nonnegative-Flonum)
> + (Float-Complex -> Nonnegative-Flonum)
> + (Number -> Nonnegative-Real)))
> (define (mag^2 x)
> - (max 0 (real-part (* x (conjugate x)))))
> + (cond [(real? x) (sqr x)]
> + [else (abs (real-part (* x (conjugate x))))]))
>
> (: vector-swap! (All (A) ((Vectorof A) Integer Integer -> Void)))
> (define (vector-swap! vs i0 i1)
> @@ -35,7 +38,9 @@
> (loop (fx+ i 1)))
> (void)))))
>
> -(: vector-scale! (case-> ((Vectorof Real) Real -> Void)
> +(: vector-scale! (case-> ((Vectorof Flonum) Flonum -> Void)
> + ((Vectorof Real) Real -> Void)
> + ((Vectorof Float-Complex) Float-Complex -> Void)
> ((Vectorof Number) Number -> Void)))
> (define (vector-scale! vs v)
> (vector-generic-scale! vs v *))
> @@ -52,25 +57,38 @@
> (loop (fx+ i 1)))
> (void)))))
>
> -(: vector-scaled-add! (case-> ((Vectorof Real) (Vectorof Real) Real -> Void)
> - ((Vectorof Real) (Vectorof Real) Real Index -> Void)
> - ((Vectorof Number) (Vectorof Number) Number -> Void)
> - ((Vectorof Number) (Vectorof Number) Number Index -> Void)))
> +(: vector-scaled-add!
> + (case-> ((Vectorof Flonum) (Vectorof Flonum) Flonum -> Void)
> + ((Vectorof Flonum) (Vectorof Flonum) Flonum Index -> Void)
> + ((Vectorof Real) (Vectorof Real) Real -> Void)
> + ((Vectorof Real) (Vectorof Real) Real Index -> Void)
> + ((Vectorof Float-Complex) (Vectorof Float-Complex) Float-Complex -> Void)
> + ((Vectorof Float-Complex) (Vectorof Float-Complex) Float-Complex Index -> Void)
> + ((Vectorof Number) (Vectorof Number) Number -> Void)
> + ((Vectorof Number) (Vectorof Number) Number Index -> Void)))
> (define (vector-scaled-add! vs0 vs1 s [start 0])
> (vector-generic-scaled-add! vs0 vs1 s start + *))
>
> -(: vector-mag^2 (case-> ((Vectorof Real) -> Nonnegative-Real)
> +(: vector-mag^2 (case-> ((Vectorof Flonum) -> Nonnegative-Flonum)
> + ((Vectorof Real) -> Nonnegative-Real)
> + ((Vectorof Float-Complex) -> Nonnegative-Flonum)
> ((Vectorof Number) -> Nonnegative-Real)))
> (define (vector-mag^2 vs)
> (define n (vector-length vs))
> - (let loop ([#{i : Nonnegative-Fixnum} 0] [#{s : Nonnegative-Real} 0])
> - (if (i . fx>= . n) s (loop (fx+ i 1) (+ s (mag^2 (unsafe-vector-ref vs i)))))))
> + (cond [(fx= n 0) (raise-argument-error 'vector-mag^2 "nonempty Vector" vs)]
> + [else
> + (define s (mag^2 (unsafe-vector-ref vs 0)))
> + (let: loop ([i : Nonnegative-Fixnum 1] [s s])
> + (cond [(i . fx< . n) (loop (fx+ i 1) (+ s (mag^2 (unsafe-vector-ref vs i))))]
> + [else (abs s)]))]))
>
> -(: vector-dot (case-> ((Vectorof Real) (Vectorof Real) -> Real)
> +(: vector-dot (case-> ((Vectorof Flonum) (Vectorof Flonum) -> Flonum)
> + ((Vectorof Real) (Vectorof Real) -> Real)
> + ((Vectorof Float-Complex) (Vectorof Float-Complex) -> Float-Complex)
> ((Vectorof Number) (Vectorof Number) -> Number)))
> (define (vector-dot vs0 vs1)
> (define n (min (vector-length vs0) (vector-length vs1)))
> - (cond [(= n 0) 0]
> + (cond [(fx= n 0) (raise-argument-error 'vector-dot "nonempty Vector" 0 vs0 vs1)]
> [else
> (define v0 (unsafe-vector-ref vs0 0))
> (define v1 (unsafe-vector-ref vs1 0))
> @@ -81,7 +99,9 @@
> (loop (fx+ i 1) (+ s (* v0 (conjugate v1))))]
> [else s]))]))
>
> -(: vector-normalize! (case-> ((Vectorof Real) -> Nonnegative-Real)
> +(: vector-normalize! (case-> ((Vectorof Flonum) -> Nonnegative-Flonum)
> + ((Vectorof Real) -> Nonnegative-Real)
> + ((Vectorof Float-Complex) -> Nonnegative-Flonum)
> ((Vectorof Number) -> Nonnegative-Real)))
> (define (vector-normalize! vs)
> (define n (vector-length vs))
> @@ -93,31 +113,51 @@
> (loop (fx+ i 1)))))
> s)
>
> -(: vector-sub-proj! (case-> ((Vectorof Real) (Vectorof Real) Any -> Nonnegative-Real)
> - ((Vectorof Number) (Vectorof Number) Any -> Nonnegative-Real)))
> +(: one (case-> (Flonum -> Nonnegative-Flonum)
> + (Real -> Nonnegative-Real)
> + (Float-Complex -> Nonnegative-Flonum)
> + (Number -> Nonnegative-Real)))
> +(define (one x)
> + (cond [(flonum? x) 1.0]
> + [(real? x) 1]
> + [(float-complex? x) 1.0]
> + [else 1]))
> +
> +(: vector-sub-proj!
> + (case-> ((Vectorof Flonum) (Vectorof Flonum) Any -> Nonnegative-Flonum)
> + ((Vectorof Real) (Vectorof Real) Any -> Nonnegative-Real)
> + ((Vectorof Float-Complex) (Vectorof Float-Complex) Any -> Nonnegative-Flonum)
> + ((Vectorof Number) (Vectorof Number) Any -> Nonnegative-Real)))
> (define (vector-sub-proj! vs0 vs1 unit?)
> (define n (min (vector-length vs0) (vector-length vs1)))
> - (define t (if unit? 1 (vector-mag^2 vs1)))
> - (unless (and (zero? t) (exact? t))
> - (define s (/ (vector-dot vs0 vs1) t))
> - (let loop ([#{i : Nonnegative-Fixnum} 0])
> - (when (i . fx< . n)
> - (define v0 (unsafe-vector-ref vs0 i))
> - (define v1 (unsafe-vector-ref vs1 i))
> - (unsafe-vector-set! vs0 i (- v0 (* v1 s)))
> - (loop (fx+ i 1)))))
> - t)
> + (cond [(fx= n 0) (raise-argument-error 'vector-sub-proj! "nonempty Vector" 0 vs0 vs1)]
> + [else
> + (define t (if unit? (one (unsafe-vector-ref vs0 0)) (vector-mag^2 vs1)))
> + (unless (and (zero? t) (exact? t))
> + (define s (/ (vector-dot vs0 vs1) t))
> + (let loop ([#{i : Nonnegative-Fixnum} 0])
> + (when (i . fx< . n)
> + (define v0 (unsafe-vector-ref vs0 i))
> + (define v1 (unsafe-vector-ref vs1 i))
> + (unsafe-vector-set! vs0 i (- v0 (* v1 s)))
> + (loop (fx+ i 1)))))
> + t]))
>
> -(: vector-zero! (case-> ((Vectorof Real) -> Void)
> +(: vector-zero! (case-> ((Vectorof Flonum) -> Void)
> + ((Vectorof Real) -> Void)
> + ((Vectorof Float-Complex) -> Void)
> ((Vectorof Number) -> Void)))
> (define (vector-zero! vs)
> (define n (vector-length vs))
> (let loop ([#{i : Nonnegative-Fixnum} 0])
> (when (i . fx< . n)
> - (unsafe-vector-set! vs i 0)
> + (define x (unsafe-vector-ref vs i))
> + (unsafe-vector-set! vs i (- x x))
> (loop (fx+ i 1)))))
>
> -(: vector-zero? (case-> ((Vectorof Real) -> Boolean)
> +(: vector-zero? (case-> ((Vectorof Flonum) -> Boolean)
> + ((Vectorof Real) -> Boolean)
> + ((Vectorof Float-Complex) -> Boolean)
> ((Vectorof Number) -> Boolean)))
> (define (vector-zero? vs)
> (define n (vector-length vs))
>
> collects/math/scribblings/math-matrix.scrbl
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/scribblings/math-matrix.scrbl
> +++ NEW/collects/math/scribblings/math-matrix.scrbl
> @@ -45,12 +45,13 @@ Like all of @racketmodname[math], @racketmodname[math/matrix] is a work in progr
> Most of the basic algorithms are implemented, but some are still in planning.
> Possibly the most useful unimplemented algorithms are
> @itemlist[@item{LUP decomposition (currently, LU decomposition is implemented, in @racket[matrix-lu])}
> - @item{@racket[matrix-solve] for upper-triangular matrices}
> + @item{@racket[matrix-solve] for triangular matrices}
> @item{Singular value decomposition (SVD)}
> @item{Eigendecomposition}
> @item{Decomposition-based solvers}
> - @item{Pseudoinverse, least-squares fitting}]
> -All of these are planned for the next Racket release.
> + @item{Pseudoinverse and least-squares solving}]
> +All of these are planned for the next Racket release, as well as fast flonum-specific matrix
> +operations and LAPACK integration.
>
> @local-table-of-contents[]
>
> @@ -66,8 +67,7 @@ From the point of view of the functions in @racketmodname[math/matrix], a @defte
>
> Technically, a matrix's entries may be any type, and some fully polymorphic matrix functions such as
> @racket[matrix-row] and @racket[matrix-map] operate on any kind of matrix.
> -Other functions, such as @racket[matrix+], require their input matrices to contain either
> - at racket[Real] or @racket[Number] values.
> +Other functions, such as @racket[matrix+], require their matrix arguments to contain numeric values.
>
> @subsection[#:tag "matrix:function-types"]{Function Types}
>
> @@ -77,17 +77,22 @@ that a return value is a matrix.
>
> Most functions that implement matrix algorithms are documented as accepting @racket[(Matrix Number)]
> values. This includes @racket[(Matrix Real)], which is a subtype. Most of these functions have a more
> -precise type than is documented. For example, @racket[matrix-conjugate] actually has the type
> - at racketblock[(case-> ((Matrix Real) -> (Matrix Real))
> - ((Matrix Number) -> (Matrix Number)))]
> -even though it is documented as having the less precise type
> - at racket[((Matrix Number) -> (Matrix Number))].
> +precise type than is documented. For example, @racket[matrix-conjugate] has the type
> + at racketblock[(case-> ((Matrix Flonum) -> (Matrix Flonum))
> + ((Matrix Real) -> (Matrix Real))
> + ((Matrix Float-Complex) -> (Matrix Float-Complex))
> + ((Matrix Number) -> (Matrix Number)))]
> +but is documented as having the type @racket[((Matrix Number) -> (Matrix Number))].
>
> Precise function types allow Typed Racket to prove more facts about @racketmodname[math/matrix]
> -client programs. In particular, it is usually easy for it to prove that matrix expressions on real
> +client programs. In particular, it is usually easy for it to prove that operations on real
> matrices return real matrices:
> @interaction[#:eval typed-eval
> (matrix-conjugate (matrix [[1 2 3] [4 5 6]]))]
> +and that operations on inexact matrices return inexact matrices:
> + at interaction[#:eval typed-eval
> + (matrix-conjugate (matrix [[1.0+2.0i 2.0+3.0i 3.0+4.0i]
> + [4.0+5.0i 5.0+6.0i 6.0+7.0i]]))]
>
> @subsection[#:tag "matrix:failure"]{Failure Arguments}
>
> @@ -102,7 +107,7 @@ For example, the (simplified) type of @racket[matrix-inverse] is
> Thus, if a failure thunk is given, the call site is required to check for return values of type
> @racket[F] explicitly.
>
> -The default failure thunk, which raises an error, has type @racket[(-> Nothing)].
> +Default failure thunks usually raise an error, and have the type @racket[(-> Nothing)].
> For such failure thunks, @racket[(U F (Matrix Number))] is equivalent to @racket[(Matrix Number)],
> because @racket[Nothing] is part of every type. (In Racket, any expression may raise an error.)
> Thus, in this case, no explicit test for values of type @racket[F] is necessary (though of course they
> @@ -110,8 +115,8 @@ may be caught using @racket[with-handlers] or similar).
>
> @subsection[#:tag "matrix:broadcasting"]{Broadcasting}
>
> -Pointwise matrix operations do not @tech{broadcast} their arguments when given matrices with
> -different sizes:
> +Unlike array operations, pointwise matrix operations @bold{do not} @tech{broadcast} their arguments
> +when given matrices with different axis lengths:
> @interaction[#:eval typed-eval
> (matrix+ (identity-matrix 2) (matrix [[10]]))]
> If you need broadcasting, use array operations:
> @@ -217,8 +222,12 @@ Like @racket[matrix], but returns a @tech{column matrix}.
> (col-matrix [])]
> }
>
> - at defproc[(identity-matrix [n Integer]) (Matrix (U 0 1))]{
> -Returns an @racket[n]×@racket[n] identity matrix; @racket[n] must be positive.
> + at defproc[(identity-matrix [n Integer] [one A 1] [zero A 0]) (Matrix A)]{
> +Returns an @racket[n]×@racket[n] identity matrix, which has the value @racket[one] on the diagonal
> +and @racket[zero] everywhere else. The height/width @racket[n] must be positive.
> + at examples[#:eval typed-eval
> + (identity-matrix 3)
> + (identity-matrix 4 1.0+0.0i 0.0+0.0i)]
> }
>
> @defproc[(make-matrix [m Integer] [n Integer] [x A]) (Matrix A)]{
> @@ -233,22 +242,28 @@ both @racket[m] and @racket[n] must be positive.
> Analogous to @racket[build-array] (and defined in terms of it).
> }
>
> - at defproc[(diagonal-matrix [xs (Listof A)]) (Matrix (U A 0))]{
> -Returns a matrix with @racket[xs] along the diagonal and @racket[0] everywhere else.
> + at defproc[(diagonal-matrix [xs (Listof A)] [zero A 0]) (Matrix A)]{
> +Returns a matrix with @racket[xs] along the diagonal and @racket[zero] everywhere else.
> The length of @racket[xs] must be positive.
> + at examples[#:eval typed-eval
> + (diagonal-matrix '(1 2 3 4 5 6))
> + (diagonal-matrix '(1.0 2.0 3.0 4.0 5.0) 0.0)]
> }
>
> @define[block-diagonal-url]{http://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices}
>
> @margin-note*{@hyperlink[block-diagonal-url]{Wikipedia: Block-diagonal matrices}}
> - at defproc[(block-diagonal-matrix [Xs (Listof (Matrix A))]) (Matrix (U A 0))]{
> -Returns a matrix with matrices @racket[Xs] along the diagonal and @racket[0] everywhere else.
> + at defproc[(block-diagonal-matrix [Xs (Listof (Matrix A))] [zero A 0]) (Matrix A)]{
> +Returns a matrix with matrices @racket[Xs] along the diagonal and @racket[zero] everywhere else.
> The length of @racket[Xs] must be positive.
> @examples[#:eval typed-eval
> (block-diagonal-matrix (list (matrix [[6 7] [8 9]])
> (diagonal-matrix '(7 5 7))
> (col-matrix [1 2 3])
> - (row-matrix [4 5 6])))]
> + (row-matrix [4 5 6])))
> + (block-diagonal-matrix (list (make-matrix 2 2 2.0+3.0i)
> + (make-matrix 2 2 5.0+7.0i))
> + 0.0+0.0i)]
> }
>
> @define[vandermonde-url]{http://en.wikipedia.org/wiki/Vandermonde_matrix}
> @@ -257,7 +272,8 @@ The length of @racket[Xs] must be positive.
> @defproc[(vandermonde-matrix [xs (Listof Number)] [n Integer]) (Matrix Number)]{
> Returns an @racket[m]×@racket[n] Vandermonde matrix, where @racket[m = (length xs)].
> @examples[#:eval typed-eval
> - (vandermonde-matrix '(1 2 3 4) 5)]
> + (vandermonde-matrix '(1 2 3 4) 5)
> + (vandermonde-matrix '(5.2 3.4 2.0) 3)]
> Using a Vandermonde matrix to find a Lagrange polynomial (the polynomial of least degree that passes
> through a given set of points):
> @interaction[#:eval untyped-eval
> @@ -384,7 +400,8 @@ Computes @racket[(matrix* M ...)] with @racket[n] arguments, but more efficientl
> @examples[#:eval untyped-eval
> ; The 100th (and 101th) Fibonacci number:
> (matrix* (matrix-expt (matrix [[1 1] [1 0]]) 100)
> - (col-matrix [0 1]))]
> + (col-matrix [0 1]))
> + (->col-matrix (list (fibonacci 100) (fibonacci 99)))]
> }
>
> @defproc[(matrix-scale [M (Matrix Number)] [z Number]) (Matrix Number)]{
> @@ -459,18 +476,19 @@ Returns array of the entries on the diagonal of @racket[M].
> (matrix ([1 2 3] [4 5 6] [7 8 9])))]
> }
>
> - at deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)]) (Matrix (U A 0))]
> - @defproc[(matrix-lower-triangle [M (Matrix A)]) (Matrix (U A 0))])]{
> + at deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)] [zero A 0]) (Matrix A)]
> + @defproc[(matrix-lower-triangle [M (Matrix A)] [zero A 0]) (Matrix A)])]{
> The function @racket[matrix-upper-triangle] returns an upper
> -triangular matrix (entries below the diagonal are zero) with
> +triangular matrix (entries below the diagonal have the value @racket[zero]) with
> entries from the given matrix. Likewise the function
> @racket[matrix-lower-triangle] returns a lower triangular
> matrix.
> - at examples[#:eval untyped-eval
> + at examples[#:eval typed-eval
> (define M (array+ (array 1) (axis-index-array #(5 7) 1)))
> M
> (matrix-upper-triangle M)
> - (matrix-lower-triangle M)]
> + (matrix-lower-triangle M)
> + (matrix-lower-triangle (array->flarray M) 0.0)]
> }
>
> @deftogether[(@defproc[(matrix-rows [M (Matrix A)]) (Listof (Matrix A))]
> @@ -1010,7 +1028,7 @@ The norm used by @racket[matrix-relative-error] and @racket[matrix-absolute-erro
> The default value is @racket[matrix-op-inf-norm].
>
> Besides being a true norm, @racket[norm] should also be @deftech{submultiplicative}:
> - at racketblock[(norm (matrix* M0 M1)) <= (* (norm A) (norm B))]
> + at racketblock[(norm (matrix* M0 M1)) <= (* (norm M0) (norm M1))]
> This additional triangle-like inequality makes it possible to prove error bounds for formulas that
> involve matrix multiplication.
>
>
> collects/math/tests/matrix-strictness-tests.rkt
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> --- OLD/collects/math/tests/matrix-strictness-tests.rkt
> +++ NEW/collects/math/tests/matrix-strictness-tests.rkt
> @@ -131,7 +131,7 @@
>
> (for*: ([M (list nonstrict-2x2-arr strict-2x2-arr)])
> (check-false (array-strict? (matrix-scale M -1)))
> - (check-true (array-strict? (matrix-expt M 0)))
> + (check-false (array-strict? (matrix-expt M 0)))
> (check-false (array-strict? (matrix-expt (array-lazy M) 1)))
> (check-false (array-strict? (matrix-expt M 2)))
> (check-false (array-strict? (matrix-expt M 3)))
>
![[logo]](/logo.png){kind=logo}