;;; cfl.ms
;;; Copyright 1984-2017 Cisco Systems, Inc.
;;; 
;;; Licensed under the Apache License, Version 2.0 (the "License");
;;; you may not use this file except in compliance with the License.
;;; You may obtain a copy of the License at
;;; 
;;; http://www.apache.org/licenses/LICENSE-2.0
;;; 
;;; Unless required by applicable law or agreed to in writing, software
;;; distributed under the License is distributed on an "AS IS" BASIS,
;;; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
;;; See the License for the specific language governing permissions and
;;; limitations under the License.

(define *fuzz* 1e-14)

(define ~=
   (lambda (x y)
      (or (= x y)
          (and (fl~= (inexact (real-part x))
                     (inexact (real-part y)))
               (fl~= (inexact (imag-part x))
                     (inexact (imag-part y)))))))

(define fl~=
   (lambda (x y)
      (cond
         [(and (fl>= (flabs x) 2.0) (fl>= (flabs y) 2.0))
          (fl~= (fl/ x 2.0) (fl/ y 2.0))]
         [(and (fl< 0.0 (flabs x) 1.0) (fl< 0.0 (flabs y) 1.0))
          (fl~= (fl* x 2.0) (fl* y 2.0))]
         [else (let ([d (flabs (fl- x y))])
                  (or (fl<= d *fuzz*)
                      (begin (printf "fl~~=: ~s~%" d) #f)))])))

(define cfl~=
   (lambda (x y)
      (and (fl~= (cfl-real-part x) (cfl-real-part y))
           (fl~= (cfl-imag-part x) (cfl-imag-part y)))))

(define zero 0.0)
(define a 1.1)
(define b +1.1i)
(define c 1.1+1.1i)
(define aa 1.21)
(define ab +1.21i)
(define ac 1.21+1.21i)
(define bb -1.21)
(define bc -1.21+1.21i)
(define cc +2.42i)

(mat cflonum?
    (not (cflonum? 3))
    (not (cflonum? 18/2))
    (not (cflonum? 1+0i))
    (not (cflonum? 23084982309482034820348023423048230482304))
    (not (cflonum? 203480234802384/23049821))
    (not (cflonum? -3/4))
    (not (cflonum? -1))
    (not (cflonum? 0))
    (not (cflonum? -12))
    (cflonum? 3.5)
    (cflonum? 1.8e-10)
    (cflonum? -3e5)
    (cflonum? -1231.2344)
    (cflonum? 3+5.0i)
    (cflonum? 1.8e10@10)
    (cflonum? -3e5+1.0i)
    (cflonum? -1.0i)
    (cflonum? +1.0i)
    (not (cflonum? 'a))
    (not (cflonum? "hi"))
    (not (cflonum? (cons 3 4)))
    (cflonum? a)
    (cflonum? b)
    (cflonum? c)
 )

(mat fl-make-rectangular
   (error? (fl-make-rectangular 3 'a))
   (error? (fl-make-rectangular 'b 4))
   (error? (fl-make-rectangular 3 -4))
   (eqv? (fl-make-rectangular 3.0 -4.0) 3.0-4.0i)
   (eqv? (fl-make-rectangular a a) c)
 )

(mat cfl-real-part
    (error? (cfl-real-part 'a))
    (error? (cfl-real-part 3/2))
    (eqv? (cfl-real-part 3.2) 3.2)
    (eqv? (cfl-real-part -1.0+2.0i) -1.0)
    (eqv? (cfl-real-part a) a)
    (eqv? (cfl-real-part c) a)
    (eqv? (cfl-real-part b) zero)
 )

(mat cfl-imag-part
    (error? (cfl-imag-part 'a))
    (error? (cfl-imag-part -3))
    (eqv? (cfl-imag-part 3.2) zero)
    (eqv? (cfl-imag-part -1.0+2.0i) 2.0)
    (eqv? (cfl-imag-part a) zero)
    (eqv? (cfl-imag-part c) a)
    (eqv? (cfl-imag-part b) a)
 )

(mat cfl-conjugate
    (error? (cfl-conjugate 'a))
    (eqv? (cfl-conjugate 3.2) 3.2)
    (eqv? (cfl-conjugate 3.2+2.0i) 3.2-2.0i)
    (eqv? (cfl-conjugate a) a)
    (eqv? (cfl-conjugate c) (+ a (- b)))
    (eqv? (cfl-conjugate b) -1.1i)
 )

(mat conjugate
    (error? (conjugate 'a))
    (eqv? (conjugate 3.2) 3.2)
    (eqv? (conjugate 3.2+2.0i) 3.2-2.0i)
 )

(mat cfl-magnitude-squared
    (error? (cfl-magnitude-squared 'a))
    (eqv? (cfl-magnitude-squared 3.2) (fl* 3.2 3.2))
    (eqv? (cfl-magnitude-squared 3.5-2.0i) 16.25)
    (fl~= (cfl-magnitude-squared 3.5@2.0) 12.25)
 )

(mat magnitude-squared
    (error? (magnitude-squared 'a))
    (eqv? (magnitude-squared 3.5) 12.25)
    (eqv? (magnitude-squared 3.5-2.0i) 16.25)
    (fl~= (magnitude-squared 3.5@2.0) 12.25)
 )

(mat cfl+
    (error? (cfl+ 'a))
    (error? (cfl+ 'a 3))
    (error? (cfl+ 'a 3 4))
    (eqv? (cfl+) zero)
    (eqv? (cfl+ a) a)
    (eqv? (cfl+ b) b)
    (eqv? (cfl+ c) c)
    (eqv? (cfl+ a b) c)
    (cfl~= (cfl+ a b c) (cfl+ a (cfl+ b c)))
    (cfl~= (cfl+ a b c a b c) (cfl+ (cfl+ a b c) (cfl+ a b c)))
    (cfl~= (cfl+ 1+2.0i 3.0) 4.0+2.0i)
    (cfl~= (cfl+ 1.0+2.2i -3.7+5.3i) -2.7+7.5i)
    (cfl~= (cfl+ 1.0+2.2i -3.7) -2.7+2.2i)
    (cfl~= (cfl+ 1.0 -3.7+5.3i) -2.7+5.3i)
    (cfl~= (cfl+ 1.0+2.2i +5.3i) 1.0+7.5i)
    (cfl~= (cfl+ +2.2i -3.7+5.3i) -3.7+7.5i)
    (cfl~= (cfl+ 26.0 2.0) 28.0)
    (test-cp0-expansion eqv? '(cfl+) zero)
    (test-cp0-expansion eqv? `(cfl+ ,a) a)
    (test-cp0-expansion eqv? `(cfl+ ,b) b)
    (test-cp0-expansion eqv? `(cfl+ ,c) c)
    (test-cp0-expansion eqv? `(cfl+ ,a ,b) c)
    (test-cp0-expansion cfl~= `(cfl+ ,a ,b ,c) (cfl+ a (cfl+ b c)))
    (test-cp0-expansion cfl~= `(cfl+ ,a ,b ,c ,a ,b ,c) (cfl+ (cfl+ a b c) (cfl+ a b c)))
    (test-cp0-expansion cfl~= '(cfl+ 1+2.0i 3.0) 4.0+2.0i)
    (test-cp0-expansion cfl~= '(cfl+ 1.0+2.2i -3.7+5.3i) -2.7+7.5i)
    (test-cp0-expansion cfl~= '(cfl+ 1.0+2.2i -3.7) -2.7+2.2i)
    (test-cp0-expansion cfl~= '(cfl+ 1.0 -3.7+5.3i) -2.7+5.3i)
    (test-cp0-expansion cfl~= '(cfl+ 1.0+2.2i +5.3i) 1.0+7.5i)
    (test-cp0-expansion cfl~= '(cfl+ +2.2i -3.7+5.3i) -3.7+7.5i)
    (test-cp0-expansion cfl~= '(cfl+ 26.0 2.0) 28.0)
 )

(mat cfl-
    (error? (cfl- 'a))
    (error? (cfl- 'a 3))
    (error? (cfl- 'a 3 4))
    (error? (cfl-))
    (eqv? (cfl- a) -1.1)
    (eqv? (cfl- b) -0.0-1.1i)
    (eqv? (cfl- c) -1.1-1.1i)
    (eqv? (cfl- a a) zero)
    (cfl~= (cfl- b b) zero)
    (cfl~= (cfl- c c) zero)
    (eqv? (cfl- c a) b)
    (cfl~= (cfl- c b) a)
    (cfl~= (cfl- a b c) (cfl- (cfl- a b) c))
    (cfl~= (cfl- a b c a b c) (cfl- a (cfl+ b c a b c)))
    (cfl~= (cfl- 1+2.0i 3.0) -2.0+2.0i)
    (cfl~= (cfl- 1.0+2.2i -3.7+5.3i) 4.7-3.1i)
    (cfl~= (cfl- 1.0+2.2i -3.7) 4.7+2.2i)
    (cfl~= (cfl- 1.0 -3.7+5.3i) 4.7-5.3i)
    (cfl~= (cfl- 1.0+2.2i +5.3i) 1.0-3.1i)
    (cfl~= (cfl- +2.2i -3.7+5.3i) 3.7-3.1i)
    (cfl~= (cfl- 26.0 2.0) 24.0)
   (andmap
     (lambda (a)
       (andmap
         (lambda (b)
           (andmap
             (lambda (c) (eqv? (cfl- a b c) (cfl- (cfl- a b) c)))
             '(0.0 -0.0)))
         '(0.0 -0.0)))
     '(0.0 -0.0))
   (let ()
     (define-syntax ff
       (syntax-rules ()
         [(_ k1 k2) (lambda (x) (eqv? (cfl- k1 x k2) (cfl- (cfl- k1 x) k2)))]))
     (andmap
       (lambda (p) (and (p +0.0) (p -0.0)))
       (list (ff +0.0 +0.0) (ff +0.0 -0.0) (ff -0.0 +0.0) (ff -0.0 -0.0))))
   (error? (cfl- 3.0 5.4 'a))
   (error? (cfl- 'a 3.0 5.4))
   (error? (cfl- 3.0 'a 5.4))
   (eqv? (cfl- 5.0 4.0 3.0 2.0) -4.0)
   (eqv? (cfl- 5.0 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0) -2.0)
   (cfl~= (cfl- 1e30 1e30 7.0) -7.0)
   
   (test-cp0-expansion eqv? `(cfl- ,a) -1.1)
   (test-cp0-expansion eqv? `(cfl- ,b) -0.0-1.1i)
   (test-cp0-expansion eqv? `(cfl- ,c) -1.1-1.1i)
   (test-cp0-expansion eqv? `(cfl- ,a ,a) zero)
   (test-cp0-expansion cfl~= `(cfl- ,b ,b) zero)
   (test-cp0-expansion cfl~= `(cfl- ,c ,c) zero)
   (test-cp0-expansion eqv? `(cfl- ,c ,a) b)
   (test-cp0-expansion cfl~= `(cfl- ,c ,b) a)
   (test-cp0-expansion cfl~= `(cfl- ,a ,b ,c) (cfl- (cfl- a b) c))
   (test-cp0-expansion cfl~= `(cfl- ,a ,b ,c ,a ,b ,c) (cfl- a (cfl+ b c a b c)))
   (test-cp0-expansion cfl~= '(cfl- 1+2.0i 3.0) -2.0+2.0i)
   (test-cp0-expansion cfl~= '(cfl- 1.0+2.2i -3.7+5.3i) 4.7-3.1i)
   (test-cp0-expansion cfl~= '(cfl- 1.0+2.2i -3.7) 4.7+2.2i)
   (test-cp0-expansion cfl~= '(cfl- 1.0 -3.7+5.3i) 4.7-5.3i)
   (test-cp0-expansion cfl~= '(cfl- 1.0+2.2i +5.3i) 1.0-3.1i)
   (test-cp0-expansion cfl~= '(cfl- +2.2i -3.7+5.3i) 3.7-3.1i)
   (test-cp0-expansion cfl~= '(cfl- 26.0 2.0) 24.0)
   (test-cp0-expansion eqv? '(cfl- 5.0 4.0 3.0 2.0) -4.0)
   (test-cp0-expansion eqv? '(cfl- 5.0 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0) -2.0)
   (test-cp0-expansion cfl~= '(cfl- 1e30 1e30 7.0) -7.0)
 )

(mat cfl*
    (error? (cfl* 'a))
    (error? (cfl* 'a 3))
    (error? (cfl* 'a 3 4))
    (eqv? (cfl*) 1.0)
    (eqv? (cfl* a) a)
    (eqv? (cfl* b) b)
    (eqv? (cfl* c) c)
    (eqv? (cfl* zero a) zero)
    (cfl~= (cfl* zero b) zero)
    (cfl~= (cfl* zero c) zero)
    (cfl~= (cfl* a a) aa)
    (cfl~= (cfl* a b) ab)
    (cfl~= (cfl* a c) ac)
    (cfl~= (cfl* b b) bb)
    (cfl~= (cfl* b c) bc)
    (cfl~= (cfl* c c) cc)
    (cfl~= (cfl* a b c) (cfl* a (cfl* b c)))
    (cfl~= (cfl* a b c a b c) (cfl* (cfl* a b c) (cfl* a b c)))
    (cfl~= (cfl* 1+2.0i 3.0) 3.0+6.0i)
    (cfl~= (cfl* 1.0+2.0i 3.0+4.0i) -5.0+10.0i)
    (cfl~= (cfl* 1.0+2.0i 3.0) 3.0+6.0i)
    (cfl~= (cfl* -2.0 3.0+4.0i) -6.0-8.0i)
    (cfl~= (cfl* 1.0+2.0i +4.0i) -8.0+4.0i)
    (cfl~= (cfl* +2.0i 3.0+4.0i) -8.0+6.0i)
    (cfl~= (cfl* 26.0 2.0) 52.0)
    (test-cp0-expansion eqv? '(cfl*) 1.0)
    (test-cp0-expansion eqv? `(cfl* ,a) a)
    (test-cp0-expansion eqv? `(cfl* ,b) b)
    (test-cp0-expansion eqv? `(cfl* ,c) c)
    (test-cp0-expansion eqv? `(cfl* ,zero ,a) zero)
    (test-cp0-expansion cfl~= `(cfl* ,zero ,b) zero)
    (test-cp0-expansion cfl~= `(cfl* ,zero ,c) zero)
    (test-cp0-expansion cfl~= `(cfl* ,a ,a) aa)
    (test-cp0-expansion cfl~= `(cfl* ,a ,b) ab)
    (test-cp0-expansion cfl~= `(cfl* ,a ,c) ac)
    (test-cp0-expansion cfl~= `(cfl* ,b ,b) bb)
    (test-cp0-expansion cfl~= `(cfl* ,b ,c) bc)
    (test-cp0-expansion cfl~= `(cfl* ,c ,c) cc)
    (test-cp0-expansion cfl~= `(cfl* ,a ,b ,c) (cfl* a (cfl* b c)))
    (test-cp0-expansion cfl~= `(cfl* ,a ,b ,c ,a ,b ,c) (cfl* (cfl* a b c) (cfl* a b c)))
    (test-cp0-expansion cfl~= '(cfl* 1+2.0i 3.0) 3.0+6.0i)
    (test-cp0-expansion cfl~= '(cfl* 1.0+2.0i 3.0+4.0i) -5.0+10.0i)
    (test-cp0-expansion cfl~= '(cfl* 1.0+2.0i 3.0) 3.0+6.0i)
    (test-cp0-expansion cfl~= '(cfl* -2.0 3.0+4.0i) -6.0-8.0i)
    (test-cp0-expansion cfl~= '(cfl* 1.0+2.0i +4.0i) -8.0+4.0i)
    (test-cp0-expansion cfl~= '(cfl* +2.0i 3.0+4.0i) -8.0+6.0i)
    (test-cp0-expansion cfl~= '(cfl* 26.0 2.0) 52.0)
 )

(mat cfl/
    (error? (cfl/ 'a))
    (error? (cfl/ 'a 3))
    (error? (cfl/ 'a 3 4))
    (error? (cfl/))
    (fl~= (cfl/ a) (fl/ a))
    (eqv? (cfl/ zero a) zero)
    (cfl~= (cfl/ zero b) zero)
    (cfl~= (cfl/ zero c) zero)
    (cfl~= (cfl/ a a) 1.0)
    (cfl~= (cfl/ b b) 1.0)
    (cfl~= (cfl/ c c) 1.0)
    (cfl~= (cfl/ aa a) a)
    (cfl~= (cfl/ ab b) a)
    (cfl~= (cfl/ ab a) b)
    (cfl~= (cfl/ ac c) a)
    (cfl~= (cfl/ ac a) c)
    (cfl~= (cfl/ bc c) b)
    (cfl~= (cfl/ bc b) c)
    (cfl~= (cfl/ cc c) c)
    (cfl~= (cfl/ a b c) (cfl/ (cfl/ a b) c))
    (cfl~= (cfl/ a b c a b c) (cfl/ a (cfl* b c a b c)))
    (cfl~= (cfl/ 3+6.0i 3.0) 1.0+2.0i)
    (cfl~= (cfl/ -5.0+10.0i 1.0+2.0i) 3.0+4.0i)
    (cfl~= (cfl/ -6.0-8.0i -2.0) 3.0+4.0i)
    (cfl~= (cfl/ 26.0 3.0-2.0i) 6.0+4.0i)
    (cfl~= (cfl/ -8.0+6.0i +2.0i) 3.0+4.0i)
    (cfl~= (cfl/ +26.0i 3.0+2.0i) 4.0+6.0i)
    (cfl~= (cfl/ 26.0 2.0) 13.0)
   (andmap
     (lambda (a)
       (andmap
         (lambda (b)
           (andmap
             (lambda (c) (eqv? (cfl/ a b c) (cfl/ (cfl/ a b) c)))
             '(1e300 1e250)))
         '(1e300 1e250)))
     '(1e300 1e250))
   (error? (cfl/ 3.0 5.4 'a))
   (error? (cfl/ 'a 3.0 5.4))
   (error? (cfl/ 3.0 'a 5.4))
   (eqv? (cfl/ 16.0 2.0 -2.0 2.0) -2.0)
   (eqv? (cfl/ 16.0 2.0 -2.0 2.0 4.0 1.0 -1.0) 0.5)
   (test-cp0-expansion eqv? `(cfl/ ,zero ,a) zero)
   (test-cp0-expansion eqv? '(cfl/ 16.0 2.0 -2.0 2.0) -2.0)
   (test-cp0-expansion eqv? '(cfl/ 16.0 2.0 -2.0 2.0 4.0 1.0 -1.0) 0.5)
   (test-cp0-expansion cfl~= `(cfl/ ,zero ,b) zero)
   (test-cp0-expansion cfl~= `(cfl/ ,zero ,c) zero)
   (test-cp0-expansion cfl~= `(cfl/ ,a ,a) 1.0)
   (test-cp0-expansion cfl~= `(cfl/ ,b ,b) 1.0)
   (test-cp0-expansion cfl~= `(cfl/ ,c ,c) 1.0)
   (test-cp0-expansion cfl~= `(cfl/ ,aa ,a) a)
   (test-cp0-expansion cfl~= `(cfl/ ,ab ,b) a)
   (test-cp0-expansion cfl~= `(cfl/ ,ab ,a) b)
   (test-cp0-expansion cfl~= `(cfl/ ,ac ,c) a)
   (test-cp0-expansion cfl~= `(cfl/ ,ac ,a) c)
   (test-cp0-expansion cfl~= `(cfl/ ,bc ,c) b)
   (test-cp0-expansion cfl~= `(cfl/ ,bc ,b) c)
   (test-cp0-expansion cfl~= `(cfl/ ,cc ,c) c)
   (test-cp0-expansion cfl~= `(cfl/ ,a ,b ,c) (cfl/ (cfl/ a b) c))
   (test-cp0-expansion cfl~= `(cfl/ ,a ,b ,c ,a ,b ,c) (cfl/ a (cfl* b c a b c)))
   (test-cp0-expansion cfl~= '(cfl/ 3+6.0i 3.0) 1.0+2.0i)
   (test-cp0-expansion cfl~= '(cfl/ -5.0+10.0i 1.0+2.0i) 3.0+4.0i)
   (test-cp0-expansion cfl~= '(cfl/ -6.0-8.0i -2.0) 3.0+4.0i)
   (test-cp0-expansion cfl~= '(cfl/ 26.0 3.0-2.0i) 6.0+4.0i)
   (test-cp0-expansion cfl~= '(cfl/ -8.0+6.0i +2.0i) 3.0+4.0i)
   (test-cp0-expansion cfl~= '(cfl/ +26.0i 3.0+2.0i) 4.0+6.0i)
   (test-cp0-expansion cfl~= '(cfl/ 26.0 2.0) 13.0)
 )

(mat cfl=
    (error? (cfl= 'a))
    (error? (cfl= 'a 3))
    (error? (cfl= 'a 3 4))
    (error? (cfl=))
    (cfl= a a)
    (cfl= b b)
    (cfl= c c)
    (cfl= (- c c) zero)
    (cfl= (+ a b) c)
    (not (cfl= a b))
    (cfl= 1.1+1.1i c)
    (cfl= c 1.1+1.1i c)
    (not (cfl= c 1.1+1.1i c a))
    (not (cfl= 3+6.0i 3.0))
    (not (cfl= 3+6.0i +6.0i))
    (cfl= 1.0+2.0i 1.0+2.0i)
    (cfl= 5.4 5.4)
 )