;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;

(in-package :maxima)

;; hacked up version of a small part of rational functions to support bigfloats,
;; using Ray's bigfloat class, so as to masquerade bigfloats as
;; lisp atoms, hence they are allowed as polynomial coefficients
;; most changes are 1-3 lines in a larger function.
;; cplus, ctimes are rewritten
;; rl():= load ("c:/lisp/rat3abfloat.lisp");

#| examples:
z:rat(x+3.1b0) 				;
z^2					;
ratdiff(%,x) 				;
%+4*x^3					;
%+1.0					;
%+1.0b0					;
%/3					;
z*2.0b0					;
exp(z) 					;
rat(2.0b0)/rat(3.0b0) 			;
(z+1)/z					;
subst(10,x,z)				;
ev(z,x=1.0b0)				;
...
maybe not what you expect?
...
rat(3.0b0)/2				;

...
definitely broken.  put known bugs here..
...



|#
;;;;;;;;;changes for bigfloat

;; from simp
(defmfun mnump (x)
  (or (numberp x)
      (typep x 'bigfloat::bigfloat) ;; new
      (and (not (atom x)) (not (atom (car x)))
	   (member (caar x) '(rat bigfloat)) )))

;; now the bigfloat is an atom.
(defvar $keepfloat t)
(setf $keepfloat t)

;; from rat3c

(defun cgcd (a b)
  (cond (modulus 1)
	((and $keepfloat (or (floatp a) (floatp b) 
			     (typep a 'bigfloat::bigfloat) (typep b 'bigfloat::bigfloat) ;; new
			     )) 1)
	(t (gcd a b))))

;; from rat3e
(defun prep1 (x &aux temp)
  (cond ((floatp x)
	 (cond ($keepfloat (cons x 1.0))
	       ((prepfloat x))))
	((bigfloatp x)(cons (bigfloat::bigfloat x) 1)) ;; new encapsulate as atom that is number!
	((integerp x) (cons (cmod x) 1))
	((rationalp x)
	      (if (null modulus)
		  (cons  (numerator x) (denominator x))
		  (cquotient (numerator x) (denominator x))))
	((atom x) (cond ((assolike x genpairs))
			(t (newsym x))))
	((and $ratfac (assolike x genpairs)))
	((eq (caar x) 'mplus)
	 (cond ($ratfac
		(setq x (mapcar #'prep1 (cdr x)))
		(cond ((every #'frpoly? x)
		       (cons (mfacpplus (mapl #'(lambda (x) (rplaca x (caar x))) x)) 1))
		      (t (do ((a (car x) (facrplus a (car l)))
			      (l (cdr x) (cdr l)))
			     ((null l) a)))))
	       (t (do ((a (prep1 (cadr x)) (ratplus a (prep1 (car l))))
		       (l (cddr x) (cdr l)))
		      ((null l) a)))))
	((eq (caar x) 'mtimes)
	 (do ((a (savefactors (prep1 (cadr x)))
		 (rattimes a (savefactors (prep1 (car l))) sw))
	      (l (cddr x) (cdr l))
	      (sw (not (and $norepeat (member 'ratsimp (cdar x) :test #'eq)))))
	     ((null l) a)))
	((eq (caar x) 'mexpt)
	 (newvarmexpt x (caddr x) t))
	((eq (caar x) 'mquotient)
	 (ratquotient (savefactors (prep1 (cadr x)))
		      (savefactors (prep1 (caddr x)))))
	((eq (caar x) 'mminus)
	 (ratminus (prep1 (cadr x))))
	((eq (caar x) 'rat)
	 (cond (modulus (cons (cquotient (cmod (cadr x)) (cmod (caddr x))) 1))
	       (t (cons (cadr x) (caddr x)))))
	((eq (caar x) 'bigfloat)(bigfloat2rat x))
	((eq (caar x) 'mrat)
	 (cond ((and *withinratf* (member 'trunc (car x) :test #'eq))
		(throw 'ratf nil))
	       ((catch 'compatvl
		  (progn
		    (setq temp (compatvarl (caddar x) varlist (cadddr (car x)) genvar))
		    t))
		(cond ((member 'trunc (car x) :test #'eq)
		       (cdr ($taytorat x)))
		      ((and (not $keepfloat)
			    (or (pfloatp (cadr x)) (pfloatp (cddr x))))
		       (cdr (ratrep* ($ratdisrep x))))
		      ((sublis temp (cdr x)))))
	       (t (cdr (ratrep* ($ratdisrep x))))))
	((assolike x genpairs))
	(t (setq x (littlefr1 x))
	   (cond ((assolike x genpairs))
		 (t (newsym x))))))

(defun newvar1 (x) 
  (declare (special varlist vlist))
  (cond ((numberp x) nil)
	((bigfloatp x) nil) ;; new
	((memalike x varlist) nil)
	((memalike x vlist) nil)
	((atom x) (putonvlist x))
	((member (caar x)
	       '(mplus mtimes rat mdifference
		 mquotient mminus bigfloat) :test #'eq)
	 (mapc #'newvar1 (cdr x)))
	((eq (caar x) 'mexpt)
	 (newvarmexpt x (caddr x) nil))
	((eq (caar x) 'mrat)
	 (and *withinratf* (member 'trunc (cdddar x) :test #'eq) (throw 'ratf '%%))
	 (cond ($ratfac (mapc 'newvar3 (caddar x)))
	       (t (mapc #'newvar1 (reverse (caddar x))))))
	(t (cond (*fnewvarsw (setq x (littlefr1 x))
			     (mapc #'newvar1 (cdr x))
			     (or (memalike x vlist)
				 (memalike x varlist)
				 (putonvlist x)))
		 (t (putonvlist x))))))


(defun pdisrep (p)
  (cond ((atom p)
	 (if (typep p 'bigfloat::bigfloat)  ;new
	     (bigfloat::real-value p) p)) ;new
    (t(pdisrep+ (pdisrep2 (cdr p) (get (car p) 'disrep))))))

;; from ratmac

(defun pzerop (x)
  (cond
    ((fixnump x) (zerop x))
    ((consp x) nil)
    ((floatp x) (zerop x))
    ((typep x 'bigfloat::bigfloat)  nil) ; new but, uh,  something like  (meqp (bigfloat::real-value x) 0 )
  ))

;; CQUOTIENT

(defun cquotient (a b)
  (cond ((equal a 0) 0)
	((null modulus)
         (if (or (floatp a) (floatp b)
		 (typep a 'bigfloat::bigfloat) ;new
		 (typep b 'bigfloat::bigfloat) ;new
                 (zerop (rem a b)))
             (div a b)
             (rat-error "quotient is not exact")))
	(t (ctimes a (crecip b)))))

;; CEXPT
;;
;; Raise an coefficient to a positive integral power. BASE should be a
;; number. POW should be a non-negative integer.
(defun cexpt (base pow)
  (unless (typep pow '(integer 0))
	   (error "CEXPT only defined for non-negative integral exponents."))
  (cond ((typep base 'bigfloat::bigfloat) ;new
	 (bigfloat::bigfloat  (power (bigfloat::real-value base) pow))) ;new
	(t  
	 (if (not modulus)
	     (expt base pow)
	   (do ((pow (ash pow -1) (ash pow -1))
		(s (if (oddp pow) base 1)))
	       ((zerop pow) s)
	     (setq base (rem (* base base) modulus))
	     (when (oddp pow) (setq s (rem (* s base) modulus))))))))

(defun cplus       (a b) ;; new, whole function
  (if (typep a 'bigfloat::bigfloat) (setf a (bigfloat::real-value a)))
  (if (typep b 'bigfloat::bigfloat) (setf b (bigfloat::real-value b)))
  (setf a (add2 a b))
  (if (bigfloatp a)(bigfloat::bigfloat a) a))

(defun ctimes       (a b) ;; new, whole function
  (if (typep a 'bigfloat::bigfloat) (setf a (bigfloat::real-value a)))
  (if (typep b 'bigfloat::bigfloat) (setf b (bigfloat::real-value b)))
  (setf a (mul2 a b))
  (if (bigfloatp a)(bigfloat::bigfloat a) a))

(defun cdifference (a b) (add a (mul -1 b))) ;new

;; from rat3e
#+ignore
(defun cdisrep (x &aux n d sign)
  (setf n (car x) d (cdr x))
  (format t "~%n=~s d=~s, type-of n =~s" n d (type-of n))
  (cond 
   ((or (if (typep n 'bigfloat::bigfloat) (setf n (bigfloat::real-value n)) nil) ;new
	(if (typep d 'bigfloat::bigfloat) (setf d (bigfloat::real-value d)) nil))

    (if (= d 1) n 
      `((mquotient) ,n ,d))) ;; new.. something better could be here..
   ((or (equal 1 d) (floatp d)) (pdisrep n))
   ((pzerop n) 0)
   (t (setq sign (cond ($ratexpand (setq n (pdisrep n)) 1)
		       ((pminusp n)
			(setq n (pdisrep (pminus n))) -1)
		       (t (setq n (pdisrep n)) 1)))
      (setq d (pdisrep d))
      (cond ((and (numberp n) (numberp d))
	     (list '(rat) (* sign n) d))

	    ((and $ratdenomdivide $ratexpand
		  (not (atom n))
		  (eq (caar n) 'mplus))
	     (fancydis n d))
	    ((numberp d)
	     (list '(mtimes ratsimp)
		   (list '(rat) sign d) n))
	    ((equal sign -1)
	     (cons '(mtimes ratsimp)
		   (cond ((numberp n)
			  (list (* n -1)
				(list '(mexpt ratsimp) d -1)))
			 (t (list sign n (list '(mexpt ratsimp) d -1))))))
	    ((equal n 1)
	     (list '(mexpt ratsimp) d -1))
	    (t (list '(mtimes ratsimp) n
		     (list '(mexpt ratsimp) d -1)))))))

(defun ratinvert (y) 
  (let ((n (car y))(d (cdr y)))
    (cond ((typep n 'bigfloat::bigfloat)
	   (cons (pctimes (bigfloat::bigfloat(div 1 (bigfloat::real-value n))) d) 
		 1))
	 
	 (t (ratalgdenom
   (cond ((pzerop n) (rat-error "quotient by zero"))
	 ((and modulus (pcoefp n))
	  (cons (pctimes (crecip n) d) 1))
	 ((and $keepfloat (floatp n))
	  (cons (pctimes (/ n) d) 1))
	 ((pminusp n) (cons (pminus d) (pminus n)))
	 (t (cons d n ))))))))

(defvar $dont_reduce_fracs t)

(defmfun ratreduce (x y &aux b) 
  (cond ((pzerop y) (rat-error "quotient by zero"))
	((pzerop x) (rzero))
	((equal y 1) (cons x 1))
	((and $keepfloat (pcoefp y) (or $float (floatp y) (pfloatp x)))
	 (cons (pctimes (quotient 1.0 y) x) 1))
	($dont_reduce_fracs 
	 (cons x y))			; new
	;; maybe poly/integer should be changed to
	;; bigfloat(1/integer)* poly ?
	(t (setq b (pgcdcofacts x y))
	   (setq b (ratalgdenom (rplacd (cdr b) (caddr b))))
	   (cond ((and modulus (pcoefp (cdr b)))
		  (cons (pctimes (crecip (cdr b)) (car b)) 1))
		 ((pminusp (cdr b))
		  (cons (pminus (car b)) (pminus (cdr b))))
		 (t b)))))

(defun cdisrep (x &aux sign)
  (cond ((typep x 'bigfloat::bigfloat) 
	 (bigfloat::real-value x))
	((and (consp (car x))(eq (caar x) 'bigfloat)) x)
	(t
  (let ((n (car x))
	(d (cdr x)))
    
    (cond 
     ((or (equal 1 d) (floatp d)) (pdisrep n))
     ((or (if (typep n 'bigfloat::bigfloat) 
	       (setf n (bigfloat::real-value n)))
	  (if (typep d 'bigfloat::bigfloat)
	      (setf d (bigfloat::real-value d))))
      (div (pdisrep n)(pdisrep d)))
     ((pzerop n) 0)
     (t (setq sign (cond ($ratexpand (setq n (pdisrep n)) 1)
			 ((pminusp n)
			  (setq n (pdisrep (pminus n))) -1)
			 (t (setq n (pdisrep n)) 1)))
	(setq d (pdisrep d))
	(cond ((and (numberp n) (numberp d))
	       (list '(rat) (* sign n) d))
	      ((and $ratdenomdivide $ratexpand
		    (not (atom n))
		    (eq (caar n) 'mplus))
	       (fancydis n d))
	      ((numberp d)
	       (list '(mtimes ratsimp)
		     (list '(rat) sign d) n))
	      ((equal sign -1)
	       (cons '(mtimes ratsimp)
		     (cond ((numberp n)
			    (list (* n -1)
				  (list '(mexpt ratsimp) d -1)))
			   (t (list sign n (list '(mexpt ratsimp) d -1))))))
	      ((equal n 1)
	       (list '(mexpt ratsimp) d -1))
	      (t (list '(mtimes ratsimp) n
		       (list '(mexpt ratsimp) d -1))))))))))
;; from numeric.lisp
(in-package :bigfloat)

#+ignore
(defmethod print-object ((x bigfloat) stream) ;+ sign redundant, except for complex??
  (let ((r (cdr (real-value x))))
    (multiple-value-bind (sign output-list)
	(if (cl:minusp (first r))
	    (values "-" (maxima::fpformat (maxima::bcons (list (cl:- (first r)) (second r)))))
	    (values "" (maxima::fpformat (maxima::bcons r))))  ;was "+"
      (format stream "~A~{~D~}" sign output-list))))

(defmethod print-object ((x bigfloat) stream)
  (let ((r  (real-value x)))
    (format stream "#~s" r))) ; signal that it is the structure

(defun bigfloat (re &optional im)
  "Convert RE to a BIGFLOAT.  If IM is given, return a COMPLEX-BIGFLOAT"
  (cond (im
	 (make-instance 'complex-bigfloat
			:real (intofp re)
			:imag (intofp im)))
	((cl:realp re)
	 (make-instance 'bigfloat :real (intofp re)))
	((cl:complexp re)
	 (make-instance 'complex-bigfloat
			:real (intofp (cl:realpart re))
			:imag (intofp (cl:imagpart re))))
	((maxima::$bfloatp re)
	 ;;(make-instance 'bigfloat :real (intofp re))
	 (make-instance 'bigfloat :real re) ;new
	 )
	((maxima::complex-number-p re 'maxima::bigfloat-or-number-p)
	 (make-instance 'complex-bigfloat
			:real (intofp (maxima::$realpart re))
			:imag (intofp (maxima::$imagpart re))))
	((typep re 'bigfloat)
	 ;; Done this way so the new bigfloat updates the precision of
	 ;; the given bigfloat, if necessary.
	 (make-instance 'bigfloat :real (real-value re)))
	((typep re 'complex-bigfloat)
	 ;; Done this way so the new bigfloat updates the precision of
	 ;; the given bigfloat, if necessary.
	 (make-instance 'complex-bigfloat
			:real (real-value re)
			:imag (imag-value re)))
	(t
	 (make-instance 'bigfloat :real (intofp re)))))