;; -*-theme-d-*-

;; Copyright (C) 2014  Tommi Höynälänmaa
;; Distributed under GNU Lesser General Public License version 3,
;; see file doc/LGPL-3.


(define-body (standard-library complex)


  (import (standard-library bitwise-arithmetic)
	  (standard-library basic-math)
	  (standard-library string-utilities))


  (define-simple-virtual-method complex
    (((re1 <real-number>) (im1 <real-number>)) <complex> (pure))
    (raise-simple 'complex:dispatch-error))

  (define-simple-virtual-method complex
    (((re1 <real>) (im1 <real>)) <complex> (pure))
    (create <complex> re1 im1))

  (define-simple-virtual-method complex
    (((re1 <real>) (im1 <integer>)) <complex> (pure))
    (create <complex> re1 (integer->real im1)))

  (define-simple-virtual-method complex
    (((re1 <real>) (im1 <rational>)) <complex> (pure))
    (create <complex> re1 (rational->real im1)))

  (define-simple-virtual-method complex
    (((re1 <integer>) (im1 <real>)) <complex> (pure))
    (create <complex> (integer->real re1) im1))

  (define-simple-virtual-method complex
    (((re1 <integer>) (im1 <integer>)) <complex> (pure))
    (create <complex> (integer->real re1) (integer->real im1)))

  (define-simple-virtual-method complex
    (((re1 <integer>) (im1 <rational>)) <complex> (pure))
    (create <complex> (integer->real re1) (rational->real im1)))

  (define-simple-virtual-method complex
    (((re1 <rational>) (im1 <real>)) <complex> (pure))
    (create <complex> (rational->real re1) im1))

  (define-simple-virtual-method complex
    (((re1 <rational>) (im1 <integer>)) <complex> (pure))
    (create <complex> (rational->real re1) (integer->real im1)))

  (define-simple-virtual-method complex
    (((re1 <rational>) (im1 <rational>)) <complex> (pure))
    (create <complex> (rational->real re1) (rational->real im1)))


  (define-simple-method real->complex (((r <real>)) <complex> (pure))
    (create <complex> r 0.0))

  (define-simple-method integer->complex (((n <integer>)) <complex> (pure))
    (create <complex> (integer->real n) 0.0))

  (define-simple-method rational->complex (((rat <rational>)) <complex> (pure))
    (create <complex> (rational->real rat) 0.0))


  (define-simple-method simplify-complex
      (((cx <complex>)) (:union <complex> <real>) pure)
    (if (= (field-ref cx 'im) 0.0)
	(field-ref cx 're)
	cx))

  
  (define-simple-method complex=? (((cx1 <complex>) (cx2 <complex>)) <boolean>
				 pure)
    (and (real=? (field-ref cx1 're) (field-ref cx2 're))
	 (real=? (field-ref cx1 'im) (field-ref cx2 'im))))

  (define-simple-method complex= (((cx1 <complex>) (cx2 <complex>)) <boolean>
				pure)
    (and (real= (field-ref cx1 're) (field-ref cx2 're))
	 (real= (field-ref cx1 'im) (field-ref cx2 'im))))

  (define-simple-method complex-integer=
      (((cx <complex>) (i <integer>)) <boolean> pure)
    (and (real-integer= (field-ref cx 're) i)
	 (real= (field-ref cx 'im) 0.0)))

  (define-simple-method integer-complex=
      (((i <integer>) (cx <complex>)) <boolean> pure)
    (and (real-integer= (field-ref cx 're) i)
	 (real= (field-ref cx 'im) 0.0)))

  (define-simple-method complex-real=
      (((cx <complex>) (r <real>)) <boolean> pure)
    (and (real= (field-ref cx 're) r)
	 (real= (field-ref cx 'im) 0.0)))

  (define-simple-method real-complex=
      (((r <real>) (cx <complex>)) <boolean> pure)
    (and (real= (field-ref cx 're) r)
	 (real= (field-ref cx 'im) 0.0)))

  (define-simple-method complex-rational=
      (((cx <complex>) (rat <rational>)) <boolean> pure)
    (and (real-rational= (field-ref cx 're) rat)
	 (real= (field-ref cx 'im) 0.0)))

  (define-simple-method rational-complex=
      (((rat <rational>) (cx <complex>)) <boolean> pure)
    (and (real-rational= (field-ref cx 're) rat)
	 (real= (field-ref cx 'im) 0.0)))

  (include-virtual-methods equal? complex=?)
  (include-virtual-methods = complex=)
  (include-virtual-methods = complex-integer=)
  (include-virtual-methods = integer-complex=)
  (include-virtual-methods = complex-real=)
  (include-virtual-methods = real-complex=)
  (include-virtual-methods = complex-rational=)
  (include-virtual-methods = rational-complex=)


  (define-simple-method real-part (((cx <complex>)) <real> (pure))
    (field-ref cx 're))


  (define-simple-method imag-part (((cx <complex>)) <real> (pure))
    (field-ref cx 'im))


  (define-simple-method make-polar (((mgn <real>) (phase <real>))
				  <complex>
				  (pure))
    (create <complex> (* mgn (r-cos phase)) (* mgn (r-sin phase))))


  (define-simple-method c-abs (((cx <complex>)) <real> (pure))
    (let ((re1 (field-ref cx 're))
	  (im1 (field-ref cx 'im)))
      (r-sqrt (real+ (real* re1 re1) (real* im1 im1)))))


  (define-simple-method c-angle (((cx <complex>)) <real> (pure))
    (r-atan2 (field-ref cx 'im) (field-ref cx 're)))


  (define-simple-method c-neg (((cx <complex>)) <complex> (pure))
    (create <complex> (r-neg (field-ref cx 're)) (r-neg (field-ref cx 'im))))


  (define-simple-method c-square (((cx <complex>)) <complex> (pure))
    (let ((re1 (field-ref cx 're))
	  (im1 (field-ref cx 'im)))
      (create <complex> (real- (r-square re1) (r-square im1))
	    (integer-real* 2 (real* re1 im1)))))


  (define-simple-method c-conj (((cx <complex>)) <complex> pure)
    (create <complex> (field-ref cx 're) (- (field-ref cx 'im))))

  
  (include-virtual-methods abs c-abs)
  (include-virtual-methods - c-neg)
  (include-virtual-methods square c-square)
  (include-virtual-methods conj c-conj)

  (define-simple-method complex+
      (((cx1 <complex>) (cx2 <complex>)) <complex> (pure))
    (create <complex>
	  (real+ (field-ref cx1 're) (field-ref cx2 're))
	  (real+ (field-ref cx1 'im) (field-ref cx2 'im))))

  (include-virtual-methods + complex+)


  (define-simple-method complex-
      (((cx1 <complex>) (cx2 <complex>)) <complex> (pure))
    (create <complex>
	  (real- (field-ref cx1 're) (field-ref cx2 're))
	  (real- (field-ref cx1 'im) (field-ref cx2 'im))))

  (include-virtual-methods - complex-)


  (define-simple-method complex*
      (((cx1 <complex>) (cx2 <complex>)) <complex> (pure))
    (let ((re1 (field-ref cx1 're))
	  (im1 (field-ref cx1 'im))
	  (re2 (field-ref cx2 're))
	  (im2 (field-ref cx2 'im)))
      (create <complex>
	    (real- (real* re1 re2) (real* im1 im2))
	    (real+ (real* re1 im2) (real* im1 re2)))))

  (include-virtual-methods * complex*)



  (define-simple-method complex/
      (((cx1 <complex>) (cx2 <complex>)) <complex> (pure))
    (let* ((re1 (field-ref cx1 're))
	   (im1 (field-ref cx1 'im))
	   (re2 (field-ref cx2 're))
	   (im2 (field-ref cx2 'im))
	   (denom (real+ (r-square re2) (r-square im2)))
	   (x (real+ (real* re1 re2) (real* im1 im2)))
	   (y (real- (real* im1 re2) (real* re1 im2))))
      (create <complex> (real/ x denom) (real/ y denom))))


  (include-virtual-methods / complex/)


  (define-simple-method real-complex+
      (((r <real>) (cx <complex>)) <complex> (pure))
    (create <complex> (real+ r (field-ref cx 're)) (field-ref cx 'im)))

  (define-simple-method complex-real+
      (((cx <complex>) (r <real>)) <complex> (pure))
    (create <complex> (real+ r (field-ref cx 're)) (field-ref cx 'im)))

  (define-simple-method integer-complex+
      (((n <integer>) (cx <complex>)) <complex> (pure))
    (create <complex> (integer-real+ n (field-ref cx 're)) (field-ref cx 'im)))

  (define-simple-method complex-integer+
      (((cx <complex>) (n <integer>)) <complex> (pure))
    (create <complex> (integer-real+ n (field-ref cx 're)) (field-ref cx 'im)))

  (define-simple-method rational-complex+
      (((rat <rational>) (cx <complex>)) <complex> (pure))
    (create <complex>
	  (rational-real+ rat (field-ref cx 're))
	  (field-ref cx 'im)))

  (define-simple-method complex-rational+
      (((cx <complex>) (rat <rational>)) <complex> (pure))
    (create <complex>
	  (rational-real+ rat (field-ref cx 're))
	  (field-ref cx 'im)))


  (define-simple-method real-complex-
      (((r <real>) (cx <complex>)) <complex> (pure))
    (create <complex> (real- r (field-ref cx 're)) (r-neg (field-ref cx 'im))))

  (define-simple-method complex-real-
      (((cx <complex>) (r <real>)) <complex> (pure))
    (create <complex> (real- (field-ref cx 're) r) (field-ref cx 'im)))

  (define-simple-method integer-complex-
      (((n <integer>) (cx <complex>)) <complex> (pure))
    (create <complex> (integer-real- n (field-ref cx 're))
	  (r-neg (field-ref cx 'im))))

  (define-simple-method complex-integer-
      (((cx <complex>) (n <integer>)) <complex> (pure))
    (create <complex> (real-integer- (field-ref cx 're) n) (field-ref cx 'im)))

  (define-simple-method rational-complex-
      (((rat <rational>) (cx <complex>)) <complex> (pure))
    (create <complex>
	  (rational-real- rat (field-ref cx 're))
	  (r-neg (field-ref cx 'im))))

  (define-simple-method complex-rational-
      (((cx <complex>) (rat <rational>)) <complex> (pure))
    (create <complex> (real-rational- (field-ref cx 're) rat) (field-ref cx 'im)))


  (define-simple-method real-complex*
      (((r <real>) (cx <complex>)) <complex> (pure))
    (create <complex> (real* r (field-ref cx 're)) (real* r (field-ref cx 'im))))

  (define-simple-method complex-real*
      (((cx <complex>) (r <real>)) <complex> (pure))
    (create <complex> (real* r (field-ref cx 're)) (real* r (field-ref cx 'im))))

  (define-simple-method integer-complex*
      (((n <integer>) (cx <complex>)) <complex> (pure))
    (create <complex>
	  (integer-real* n (field-ref cx 're))
	  (integer-real* n (field-ref cx 'im))))

  (define-simple-method complex-integer*
      (((cx <complex>) (n <integer>)) <complex> (pure))
    (create <complex>
	  (integer-real* n (field-ref cx 're))
	  (integer-real* n (field-ref cx 'im))))

  (define-simple-method rational-complex*
      (((rat <rational>) (cx <complex>)) <complex> (pure))
    (create <complex>
	  (rational-real* rat (field-ref cx 're))
	  (rational-real* rat (field-ref cx 'im))))

  (define-simple-method complex-rational*
      (((cx <complex>) (rat <rational>)) <complex> (pure))
    (create <complex>
	  (rational-real* rat (field-ref cx 're))
	  (rational-real* rat (field-ref cx 'im))))


  (define-simple-method real-complex/
      (((r <real>) (cx <complex>)) <complex> (pure))
    (let* ((re2 (field-ref cx 're))
	   (im2 (field-ref cx 'im))
	   (denom (real+ (r-square re2) (r-square im2)))
	   (x (real* r re2))
	   (y (r-neg (real* r im2))))
      (create <complex> (real/ x denom) (real/ y denom))))

  (define-simple-method complex-real/
      (((cx <complex>) (r <real>)) <complex> (pure))
    (create <complex> (real/ (field-ref cx 're) r) (real/ (field-ref cx 'im) r)))

  (define-simple-method integer-complex/
      (((n <integer>) (cx <complex>)) <complex> (pure))
    (let* ((re2 (field-ref cx 're))
	   (im2 (field-ref cx 'im))
	   (denom (real+ (r-square re2) (r-square im2)))
	   (x (integer-real* n re2))
	   (y (r-neg (integer-real* n im2))))
      (create <complex> (real/ x denom) (real/ y denom))))

  (define-simple-method complex-integer/
      (((cx <complex>) (n <integer>)) <complex> (pure))
    (create <complex>
	  (real-integer/ (field-ref cx 're) n)
	  (real-integer/ (field-ref cx 'im) n)))

  (define-simple-method rational-complex/
      (((rat <rational>) (cx <complex>)) <complex> (pure))
    (let* ((re2 (field-ref cx 're))
	   (im2 (field-ref cx 'im))
	   (denom (real+ (r-square re2) (r-square im2)))
	   (x (rational-real* rat re2))
	   (y (r-neg (rational-real* rat im2))))
      (create <complex> (real/ x denom) (real/ y denom))))

  (define-simple-method complex-rational/
      (((cx <complex>) (rat <rational>)) <complex> (pure))
    (create <complex>
	  (real-rational/ (field-ref cx 're) rat)
	  (real-rational/ (field-ref cx 'im) rat)))

  (include-virtual-methods + real-complex+)
  (include-virtual-methods + complex-real+)
  (include-virtual-methods + integer-complex+)
  (include-virtual-methods + complex-integer+)
  (include-virtual-methods + rational-complex+)
  (include-virtual-methods + complex-rational+)

  (include-virtual-methods - real-complex-)
  (include-virtual-methods - complex-real-)
  (include-virtual-methods - integer-complex-)
  (include-virtual-methods - complex-integer-)
  (include-virtual-methods - rational-complex-)
  (include-virtual-methods - complex-rational-)

  (include-virtual-methods * real-complex*)
  (include-virtual-methods * complex-real*)
  (include-virtual-methods * integer-complex*)
  (include-virtual-methods * complex-integer*)
  (include-virtual-methods * rational-complex*)
  (include-virtual-methods * complex-rational*)

  (include-virtual-methods / real-complex/)
  (include-virtual-methods / complex-real/)
  (include-virtual-methods / integer-complex/)
  (include-virtual-methods / complex-integer/)
  (include-virtual-methods / rational-complex/)
  (include-virtual-methods / complex-rational/)


  (define-simple-method c-sqrt (((cx <complex>)) <complex> (pure))
    (let ((abs-value (c-abs cx))
	  (phase (c-angle cx)))
      (make-polar (r-sqrt abs-value) (/ phase 2))))


  (define-simple-method c-exp (((cx <complex>)) <complex> (pure))
    (let ((coeff (r-exp (field-ref cx 're)))
	  (im1 (field-ref cx 'im)))
      (create <complex> (* coeff (r-cos im1)) (* coeff (r-sin im1)))))


  (define-simple-method c-log (((cx <complex>)) <complex> (pure))
    (let ((abs-value (c-abs cx))
	  (phase (c-angle cx)))
      (create <complex> (r-log abs-value) phase)))


  (define-simple-method c-log10 (((cx <complex>)) <complex> (pure))
    (/ (c-log cx) (r-log 10.0)))


  (define-simple-method c-expt (((base <complex>) (exponent <complex>))
			      <complex>
			      (pure))
    (c-exp (* exponent (c-log base))))


  (define-simple-method c-sin (((cx <complex>)) <complex> (pure))
    (let* ((re1 (field-ref cx 're))
	   (im1 (field-ref cx 'im))
	   (exp-term1 (r-exp im1))
	   (exp-term2 (r-exp (- im1)))
	   (cos-term (r-cos re1))
	   (sin-term (r-sin re1))
	   (a (* (- exp-term2 exp-term1) cos-term))
	   (b (* (+ exp-term2 exp-term1) sin-term)))
      (create <complex> (/ b 2) (- (/ a 2)))))


  (define-simple-method c-cos (((cx <complex>)) <complex> (pure))
    (let* ((re1 (field-ref cx 're))
	   (im1 (field-ref cx 'im))
	   (exp-term1 (r-exp im1))
	   (exp-term2 (r-exp (- im1)))
	   (cos-term (r-cos re1))
	   (sin-term (r-sin re1))
	   (a (* (+ exp-term2 exp-term1) cos-term))
	   (b (* (- exp-term2 exp-term1) sin-term)))
      (create <complex> (/ a 2) (/ b 2))))


  (define-simple-method c-tan (((cx <complex>)) <complex> (pure))
    (let* ((re1 (field-ref cx 're))
	   (im1 (field-ref cx 'im))
	   (exp-term1 (r-exp im1))
	   (exp-term2 (r-exp (- im1)))
	   (cos-term (r-cos re1))
	   (sin-term (r-sin re1))
	   (a (* (- exp-term2 exp-term1) cos-term))
	   (b (* (+ exp-term2 exp-term1) sin-term))
	   (c (* (+ exp-term2 exp-term1) cos-term))
	   (d (* (- exp-term2 exp-term1) sin-term)))
      (/ (create <complex> b (- a)) (create <complex> c d))))


  (define-simple-method c-asin (((cx <complex>)) <complex> (pure))
    (let* ((unit (create <complex> 1.0 0.0))
	   (imag-unit (create <complex> 0.0 1.0))
	   (term1 (* imag-unit cx))
	   (sq (c-square cx))
	   (term2 (c-sqrt (- unit sq)))
	   (log-term (c-log (+ term1 term2))))
      (* (create <complex> 0.0 -1.0) log-term)))


  (define-simple-method c-acos (((cx <complex>)) <complex> (pure))
    (let* ((unit (create <complex> 1.0 0.0))
	   (imag-unit (create <complex> 0.0 1.0))
	   (sq (c-square cx))
	   (term2 (* imag-unit (c-sqrt (- unit sq))))
	   (log-term (c-log (- cx term2))))
      (* imag-unit log-term)))


  (define-simple-method c-atan (((cx <complex>)) <complex> (pure))
    (let* ((unit (create <complex> 1.0 0.0))
	   (imag-unit (create <complex> 0.0 1.0))
	   (iz (create <complex> (- (field-ref cx 'im)) (field-ref cx 're)))
	   (term1 (c-log (- unit iz)))
	   (term2 (c-log (+ unit iz))))
      (* (create <complex> 0.0 0.5) (- term1 term2))))


  (define-simple-method c-sinh (((cx <complex>)) <complex> (pure))
    (let ((term1 (c-exp cx))
	  (term2 (c-exp (- cx))))
      (* 0.5 (- term1 term2))))


  (define-simple-method c-cosh (((cx <complex>)) <complex> (pure))
    (let ((term1 (c-exp cx))
	  (term2 (c-exp (- cx))))
      (* 0.5 (+ term1 term2))))


  (define-simple-method c-tanh (((cx <complex>)) <complex> (pure))
    (let* ((exp-term1 (c-exp cx))
	   (exp-term2 (c-exp (- cx))))
      (/ (- exp-term1 exp-term2) (+ exp-term1 exp-term2))))


  (define-simple-method c-asinh (((cx <complex>)) <complex> (pure))
    (let* ((unit (create <complex> 1.0 0.0))
	   (sq (c-square cx))
	   (term2 (c-sqrt (+ sq unit))))
      (c-log (+ cx term2))))


  (define-simple-method c-acosh (((cx <complex>)) <complex> (pure))
    (let* ((unit (create <complex> 1.0 0.0))
	   (sq (c-square cx))
	   (term2 (c-sqrt (- sq unit))))
      (c-log (+ cx term2))))


  (define-simple-method c-atanh (((cx <complex>)) <complex> (pure))
    (let* ((unit (create <complex> 1.0 0.0)))
      (* 0.5 (c-log (/ (+ unit cx) (- unit cx))))))
	   

  (define-simple-method r-complex-log (((r <real>)) <complex> pure)
    (if (>= r 0.0)
	(complex (r-log r) 0.0)
	(complex (r-log (r-abs r)) gl-r-pi)))

 
  (define-simple-method complex-real-expt (((base <complex>) (exponent <real>))
					 <complex> pure)
    (c-exp (* exponent (c-log base))))


  (define-simple-method real-complex-expt (((base <real>) (exponent <complex>))
					 <complex> pure)
    (c-exp (* exponent (r-complex-log base))))


  (define-simple-method r-log-neg (((r <real>)) <complex> pure)
    (complex (r-log (r-abs r)) gl-r-pi))


  (define-simple-method r-log10-neg (((r <real>)) <complex> pure)
    (complex (r-log10 (r-abs r)) gl-r-pi/ln10))


  (define-simple-method r-complex-expt (((base <real>) (exponent <real>))
				      <complex> pure)
    (c-exp (* exponent (r-complex-log base))))


  (define-simple-method c-exp2 (((exponent <complex>)) <complex> pure)
    (c-exp (* exponent gl-r-ln2)))


  ;; Return 1.0+0.0i also for 0^0.
  (define-simple-method c-nonneg-int-expt
      (((cx-base <complex>) (i-exponent <integer>)) <complex> pure)
    (assert (>= i-exponent 0))
    (let-mutable ((i-cur-exponent <integer> i-exponent)
		  (cx-cur-base <complex> cx-base)
		  (cx-sum <complex> (complex 1.0 0.0)))
      (until ((or (equal? i-cur-exponent 0) (equal? i-cur-exponent 1)))
	(if (not (equal? (bitwise-and i-cur-exponent 1) 0))
	    (set! cx-sum (* cx-sum cx-cur-base)))
	(set! cx-cur-base (square cx-cur-base))
	(set! i-cur-exponent (ash i-cur-exponent -1)))
      (cond
       ((= i-cur-exponent 0) cx-sum)
       ((= i-cur-exponent 1) (* cx-sum cx-cur-base))
       (else (raise-simple 'c-nonneg-int-expt:internal-error)))))


  (define-simple-method c-int-expt (((cx-base <complex>) (i-expt <integer>))
				  <complex> pure)
    (if (>= i-expt 0)
	(c-nonneg-int-expt cx-base i-expt)
	(/ 1.0 (c-nonneg-int-expt cx-base (i-neg i-expt)))))


  (define-simple-method complex->string
      (((cx <complex>)) <string> (pure))
    (let* ((re1 (field-ref cx 're))
	   (im1 (field-ref cx 'im))
	   (ch (if (and (>= im1 0) (finite? im1)) "+" "")))
      (string-append (real->string re1)
		     ch
		     (real->string im1)
		     "i")))

  
  (include-virtual-methods object->string complex->string))