;; -*- Mode:Common-Lisp;Package:mma; Base:10 -*-
(in-package :mma)
;;; Extended numbers ...

;;; generalized Complex numbers stored as structures containing two 
;;; generalized numbers
;;; which are interpreted as real and imaginary parts of a quantity.
;;; Often, each part is a bigfloat, but could be some CL numeric
;;; quantity (not a complex, please.)

;; declare the structure of complex 
;; we use a rectangular form, but could presumably use (r, theta)
;; if we wished to change a few of the routines.

(shadow '(complex) :mma)

(deftype supercomplex () '(or lisp::complex mma::complex))  

;; if we shadow the old "complex", we can get away with this...

(defstruct (complex  ;; this is in the mma package only
	    (:constructor complex (real imag))
	    (:print-function complexprintfunction))
real imag)

(defun complexprintfunction (x s pl)
  (declare (ignore pl))
  (format s "(~s + ~s I)" (Realpart x)(Imagpart x)))

(defun Realpart(x)(typecase x
			    (lisp::number (realpart x))
			    (complex (complex-real x))
			    (t `(Re ,x))))
(defun Imagpart(x)(typecase x
			    (lisp::number (imagpart x))
			    (complex (complex-imag x))
			    (t `(Im ,x))))

;; test for equality of two complexes

(defun complex-eql(x y)
  (and (eql (complex-real x)(complex-real y))
       (eql (complex-imag x)(complex-imag y))))

;; takes whatever x is, and returns a complex bigfloat.  
(defun complex-bigfloat-convert(x)
;;; somehow check conversion possibility?
	 (make-complex
		   (bigfloat-convert (Realpart x))
		   (bigfloat-convert (Imagpart x))))
	 
;;; re-examine this... what we might want to do is determine
;;; the types of a and b, and depending on that, use various addition
;;; routines. contract the answer, sometimes..
;;;; Rearing its head --   Inheritance again..  see ~/newarith/gen.l

(defun complex-+ (a b)
;;; we want check for Imagpart = 0, and then return 0
(let ((r  (bigfloat-+ (Realpart a)(Realpart b)))
      (i  (bigfloat-+ (Imagpart a)(Imagpart b))))
(if (bigfloat-zerop i) r
  (make-complex r i))


;;;; rewrite beyond here...


;; cbigfloat-* multiplies two cbigfloats of arbitrary (perhaps different)
;; precisions, and return a cbigfloat of global (bigfloat-bin-prec)
;; precision.
;; (note:  (a+bi)*(c+di)=  (a*c-b*d) +(a*d+b*c)i )

(defun complex-* (r s) 
  (let ((a (cbf-real r))
	(b (cbf-imag r))
	(c (cbf-real s))
	(d (cbf-imag s)))
    (make-cbig
     (bigfloat-- (bigfloat-* a c)(bigfloat-* b d))
     (bigfloat-+ (bigfloat-* a d)(bigfloat-* b c)))))

;; cbigfloat-/ divides two cbigfloats of arbitrary (perhaps different)
;; precisions, and return a bigfloat of global (bigfloat-bin-prec)
;; precision.
;; (a+bi)/(c+di) = 1/(c^2+d^2) * ( (ac+bd) + (bc-ad) i).

(defun cbigfloat-/ (r s)
    (let ((a (cbf-real r))
	(b (cbf-imag r))
	(c (cbf-real s))
	(d (cbf-imag s))
	denom)
      (setq denom (bigfloat-+ (bigfloat-* c c)(bigfloat-* d d)))
      (make-cbig
     (bigfloat-/(bigfloat-+ (bigfloat-* a c)(bigfloat-* b d))denom)
     (bigfloat-/(bigfloat-- (bigfloat-* b c)(bigfloat-* a d))denom))))
    

;; coerce-cbigfloat is like the CL function coerce, but allows
;; the first argument to be of type bigfloat.  It converts the
;; first argument to  ratio, double-float, single-float (= float),
;; integer (= bignum or fixnum)

(defun coerce-cbigfloat(x typ) 
  (cond ;; convert bigfloat to bigfloat of global precision
        ((eq typ 'cbigfloat) 
	 (cbigfloat-convert x))
	(complex (coerce-bigfloat (cbf-real x) typ)
		 (coerce-bigfloat (cbf-imag x) typ)))


;; compute x-y or -x  (if y is missing)

(defun cbigfloat-- (x &optional (y nil))
  (cond (y (cbigfloat-+ x (cbigfloat-- y)));; compute x-y
	((cbigfloat-zerop x) x)
	(t (makecbig
	    (bigfloat-- (cbf-real x))
	    (bigfloat-- (cbf-imag x))))))



;; compute p^n for bigfloat p.  nn is positive or negative integer.
;; Note that we do NOT allow bigfloat exponent here. (use log/exp for that)

(defun cbigfloat-expt (p nn) 
  (format t "cbigfloat-expt not implemented yet")
)

(defun cbigfloat-abs (x) 
  (let ((c (cbf-real x))(d (cbf-imag x)))
    (bigfloat-+ (bigfloat-* c c)(bigfloat-* d d))))


;;how to print these guys?