diff --git a/math-lib/math/private/matrix/matrix-basic.rkt b/math-lib/math/private/matrix/matrix-basic.rkt index 1723e21..516881c 100644 --- a/math-lib/math/private/matrix/matrix-basic.rkt +++ b/math-lib/math/private/matrix/matrix-basic.rkt @@ -309,15 +309,73 @@ ((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)))) +(define matrix-cos-angle + (let () + (define (mag² [x : Number]) (if (real? x) (sqr x) (+ (sqr (real-part x)) (sqr (imag-part x))))) + (define (inf=>1 [x : Real]) : Flonum (if (eq? x +inf.0) 1. (if (eq? x -inf.0) -1. 0.))) + (: cinf=>1 (case-> (-> Real Flonum) + (-> Float-Complex Float-Complex) + (-> Number Number))) + (define (cinf=>1 x) + (if (real? x) + (inf=>1 x) + (make-flrectangular (inf=>1 (real-part x)) + (inf=>1 (imag-part x))))) + (: unit-bound (case-> (-> Flonum Flonum) + (-> Real Real))) + (define (unit-bound x) + (if (flonum? x) + (flmin (flmax -1. x) 1.) + (min (max -1 x) 1))) + (: inner (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 (inner A* B*) + (define nA (sqrt (array-all-sum (inline-array-map mag² A*)))) + (define nB (sqrt (array-all-sum (inline-array-map mag² B*)))) + + (define result (/ (matrix-dot A* B*) (* nA nB))) + + (if (real? result) + (unit-bound result) + (make-rectangular (unit-bound (real-part result)) + (unit-bound (imag-part result))))) + (λ (A B) + (define mA (array-strict (inline-array-map nonstupid-magnitude A))) + (define mB (array-strict (inline-array-map nonstupid-magnitude B))) + (define mxA (array-all-max mA)) + (define mxB (array-all-max mB)) + + (cond + [(and (rational? mxA) (positive? mxA) + (rational? mxB) (positive? mxB)) + (define A* (inline-array-map (λ (x) (/ x mxA)) A)) + (define B* (inline-array-map (λ (x) (/ x mxB)) B)) + + (inner A* B*)] + [(or (nan? mxA) (nan? mxB) (= 0 mxA) (= 0 mxB)) + (/ (matrix-dot A B) (* mxA mxB))] + [else + (define A* (if (rational? mxA) + (inline-array-map (λ (x) (/ x mxA)) A) + (array-strict (inline-array-map cinf=>1 A)))) + (define B* (if (rational? mxB) + (inline-array-map (λ (x) (/ x mxB)) B) + (array-strict (inline-array-map cinf=>1 B)))) + (inner A* B*)])))) (: 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))) +(define (matrix-angle A B) + (define ca (matrix-cos-angle A B)) + (cond + [(flonum? ca) (flacos ca)] + [(real? ca) + (if (eq? ca 1) 0 (flacos (fl ca)))] + [else (acos ca)])) (: matrix-normalize (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum)) diff --git a/math-lib/math/private/matrix/matrix-operator-norm.rkt b/math-lib/math/private/matrix/matrix-operator-norm.rkt index d6e9d7c..0c0f11e 100644 --- a/math-lib/math/private/matrix/matrix-operator-norm.rkt +++ b/math-lib/math/private/matrix/matrix-operator-norm.rkt @@ -79,9 +79,7 @@ See "How to Measure Errors" in the LAPACK manual for more details: ;(matrix-min-singular-value (matrix* (matrix-hermitian M) R)) (error 'unimplemented)) -(: matrix-basis-angle (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum) - ((Matrix Real) (Matrix Real) -> Real) - ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex) +(: matrix-basis-angle (case-> ((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) diff --git a/math-test/math/tests/matrix-tests.rkt b/math-test/math/tests/matrix-tests.rkt index 2821b7f..0b804be 100644 --- a/math-test/math/tests/matrix-tests.rkt +++ b/math-test/math/tests/matrix-tests.rkt @@ -661,7 +661,39 @@ (check-exn exn:fail:contract? (λ () (matrix-dot a (matrix [[1]])))) (check-exn exn:fail:contract? (λ () (matrix-dot (matrix [[1]]) a)))) -;; TODO: matrix-angle +;; matrix-angle +;; guarantee that matrix-cos-angle on real is in [-1 .. 1] except for zero? +(for: ([i (in-range 10)]) + (define n (+ 1 (random 3))) + (define m (+ 1 (random 4))) + (define A (random-matrix n m)) + (define B (random-matrix n m)) + (define ca (matrix-cos-angle A B)) + (check-true (or (and (rational? ca) (<= ca 1.)) + (matrix-zero? A) (matrix-zero? B)))) +(let ([A (matrix [[-1.65e+145 -9.55e+152]])] + [B (matrix [[-2.05e-130 -1.88e-122]])]) + (check-true (<= (abs (matrix-cos-angle A B)) 1.))) +;; check matrix-cos-angle doesn't overflow easily +(let ([A (matrix [[-1.20e-21+8.16e+173i -3.18e+280-4.15e+275i]])] + [B (matrix [[-2.51e-270-4.38e-68i -1.65e+232+3.06e+227i]])]) + (define ca (matrix-cos-angle A B)) + ;; 0.9999999995008538+3.159576900273938e-05i + (check-true (<= (magnitude ca) (flnext 1.)))) +(let ([A (matrix-scale (matrix [[1 1]]) 1e270+3.i)] + [B (matrix-scale (matrix [[9 0]]) 1e270+3.i)]) + (define ca (matrix-cos-angle A B)) + ;; 0.9999999995008538+3.159576900273938e-05i + (check-true (<= (magnitude (- ca (make-rectangular (sqrt 1/2) 0.))) epsilon.0))) +;; check matrix-cos-angle treats axes with inf as principal axes +(check-equal? (matrix-cos-angle (matrix [[1 5] [2 9]]) + (matrix [[2 8] [3 -inf.0]])) + (matrix-cos-angle (matrix [[1 5] [2 9]]) + (matrix [[0. 0.] [0. -1.0]]))) +(check-equal? (matrix-cos-angle (matrix [[2 9]]) + (matrix [[3 -8+inf.0i]])) + (matrix-cos-angle (matrix [[2 9]]) + (matrix [[0. -0+1.0i]]))) ;; TODO: matrix-normalize