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

(in-package :maxima)


;;; REVISED 1/2016 to use decimal fpprec, input, output  and inside.
;;; RJF.
;; rl():=load("c:/lisp/decimalfloat.lisp");

;;;  problems to be expected:  

;;; 1. Need to search every other file for treatment
;;; of bigfloats that expects that fpprec is power-of-2. There are a few
;;; instances I found on quick glance. I tried sin(bf), it went into loop. 

;;; 2. I started to write the code so that the decimal point is on
;;; the right of the mantissa. That is  1/2  is  5 X 10^(-1).
;;; I suppose I could do this, for decimal precision 4, do 5000 x 10^-(-4).

;;; This latter version has the property that precision can be
;;; deduced from the length of the fraction part, except for 0x10^0.
;;; The previous binary version put the binary point on the left.
;;;   so in binary precision 10,  the code would be 1024 x 2^(-11)
;;;  we store the precision in the header, nevertheless, for speed.
;;;;;;;;;;;;;;


;; some decimal functions

(defun $decimalfptrim(B)			; B is a bigfloat
  (cons (car B)(decimalfptrim (cdr B))))

;; Need a  better version of decimalfptrim to carefully round instead of truncate
#+ignore
(defun decimalfptrim( f)  ;;  f is (significand exponent)
  (let ((s (car f))
	(e (cadr f)))
    (cond((eq s 0)(list 0 0))
	 (t
	  ;; first remove trailing zeros
	  (while (= (mod s 10) 0)(setf s (truncate s 10))(incf e))
	  ;; there may still be more digits than allowed by fpprec
	  
	  (let ((adjust (- (flatsize10 s) fpprec)))
	    (if (<= adjust 0) (list s e)
	      (list (truncate s (expt 10 adjust)) (+ e  adjust))))
	  ))))

;; the more careful version

(defun decimalfptrim( f)  ;;  f is (significand exponent)
  (let ((s (car f))
	(e (cadr f)))
    (cond((eq s 0)(list 0 0))
	 (t
	  ;; first remove trailing zeros
	  (while (= (mod s 10) 0)(setf s (truncate s 10))(incf e))
	  ;; there may still be more digits than allowed by fpprec
	  (let ((adjust (- (flatsize10 s) fpprec)))
;;;	    (format t"~%in decimalfptrim, s=~s, adjust=~s~%" s adjust)
	    (if (<= adjust 0) (list s e)
	      (decimalround s e adjust)	    ))
	  ))))

(defun decimalround (s e ad) ; we want to remove ad digits from the right
  (let*((guard 0)	     ; this is the guard digit   
	(sign (signum s)))
					;significand with one extra (guard) digit
					;sticky has the rest
    (multiple-value-bind 
	(swithguard sticky) 
	(truncate s (expt 10 (1- ad)))
      ;; now have last digit of s, the one to drop off,  at left end of swithguard
      (setf guard (mod swithguard 10))
;;;      (format t "~% in decimalround swg=~s guard=~s sign=~s~%" swithguard guard sign) 
      ;; because (mod 13 10) is 3, (mod -13 10) is 7...
      ;;  round 13b0 to 1b0,  round -13 to -1b0  in 1 digit..
    
      ;; because (mod 18 10) is 8, (mod -18 10) is 2
      ;;  round 18 to 2  round -18 to -2
    
      ;; now consider round-to-even
      ;; because (mod 15 10) is 5   and 1 is odd, round to 2.
      ;; because (mod -15 10) is 5  and 1 is odd, round to -2.

      (decimalfptrim
       (list (cond ((= sign 1) 
		    (cond ((< guard 5) (truncate swithguard 10)) ; no sweat
			  ;; in case we have 96 b0 we round up and get 10b0. too many digits
			  ;; so in the case of xxxx99X, or worse 9999.999X we need to re-trim.
			  ((or (> guard 5)
			       (and (= guard 5)(> sticky 0)))
			   (1+ (truncate swithguard 10))) ;no sweat again.
			  (t
			   ;; guard digit is 5 and sticky is 0. (cannot be negative if sign is 1)
			   ;; so find the even digit at or above this number
			   (let* ((snoguard(truncate swithguard 10))
				  (checkthis (mod snoguard 10)))
			     (if (evenp checkthis)(list snoguard (+  e ad)) ; this is already even
			       (1+ snoguard)))))) 
		   (t		
		  ;;; i.e. sign=-1, because sign=0 can't happen
		    (cond ((> guard 5) (truncate swithguard 10))
			  ((or (< guard 5)
			       (and (= guard 5)(< sticky 0)))
			   (1- (truncate swithguard 10)))
			  ;; guard digit is 5 and sticky is 0 (cannot be pos if sign is -1)
			  (t    (let* ((snoguard(truncate swithguard 10))
				       (checkthis (mod snoguard 10)))
				  (if (evenp checkthis)snoguard ; this is already even
				    (1- snoguard)))))))
	     (+ e ad))))))
	  
(defun flatsize10(r) (if (= r 0) 1 
		       (ceiling (* #.(log 2)(integer-length (abs r))) #.(log 10 ))))

(defvar *m)

;; *DECFP = T if the computation is being done in decimal radix.  NIL implies
;; base 2.  
;;;originally, Decimal radix was used only during output.
;;; now -- in this file it is used all the time...  1/2016/RJf




(declare-top (special *cancelled $float $bfloat $ratprint $ratepsilon $domain $m1pbranch))

;; Representation of a Bigfloat:  ((BIGFLOAT SIMP precision) mantissa exponent)
;;  NOW -- number of decimal digits 1/2016
;;		precision = (integer-length mantissa)
;;  NOW        The actual number represented is (* fraction (^ 10 exponent)).

(defun fpprec1 (assign-var q)
  (declare (ignore assign-var))
  (if (or (not (fixnump q)) (< q 1))
      (merror (intl:gettext "fpprec: value must be a positive integer; found: ~M") q))
  (setq fpprec q  ;; now the same as $fpprec.  1/2016
	))


(defun dim-bigfloat (form result)
  (let (($lispdisp nil))
    (dimension-atom (maknam (fpformat form)) result)))

;; Converts the bigfloat L to list of digits including |.| and the
;; exponent marker |b|. The number of significant digits is controlled
;; by $fpprintprec.

#+ignore
(defun fpformat (l) (append (explodec (cadr l))'(|b|) (explodec (caddr l))))

;; uses L marker instead of b.   123b0 prints 1.23L2

(defun fpformat (l) 
  (cond ((= (cadr l) 0) '(|0| |.| |L| |0|))
	((> (cadr l) 0)
	 (let ((M(explodec (format nil "~s"(cadr l)))))
	   (append (cons (car M)(cons '|.| (cdr M))) '(|L|)(explodec(+ (caddr l)(length M) -1)))))


	(t  (let ((M(explodec (format nil "~s"(cadr l)))))
	      (append (cons '|-|(cons (car M)(cons '|.| (cdr M))))'(|L|)(explodec(+ (caddr l)(length M) -1)))))))


;; Tells you if you have a bigfloat object.  BUT, if it is a bigfloat,
;; it will normalize it by making the precision of the bigfloat match
;; the current precision setting in fpprec.  And it will also convert
;; bogus zeroes (mantissa is zero, but exponent is not) to a true
;; zero.

(defun bigfloatp (x)
  ;; A bigfloat object looks like '((bigfloat simp <prec>) <mantissa> <exp>)
  ;; with decimal, we don't change the significand unless it has too many digits.
  (if (and (consp x)(eq (caar x) 'bigfloat))
      ($decimalfptrim x) nil))

(defun bigfloat2rat (x) ;; decimal version
  (let ((k (* (cadr x)(expt 10 (caddr x)))))
    (cons (numerator k)(denominator k)))) ; do simplification later

;; Convert a floating point number into a bigfloat.
;; Probably we don't do this a lot, so we just do it
;; with a short program.
(defun floattofp (x)
  (cdr ($bfloat ($rationalize x))))

;; Convert a bigfloat into a floating point number.
(defmfun fp2flo (l)
  (coerce  (* (cadr l)(expt 10 (caddr l))) 'double-float))

(defun bcons (s)
  `((bigfloat simp ,fpprec) . ,(decimalfptrim s)))

(defmfun $bfloat (x)
  (let (y)
    (cond ((bigfloatp x) x)
	  ((or (numberp x)
	       (member x '($%e $%pi $%gamma) :test #'eq))
	   (bcons (intofp x)))
	  ((or (atom x) (member 'array (cdar x) :test #'eq))
	   (if (eq x '$%phi)
	       ($bfloat '((mtimes simp)
			  ((rat simp) 1 2)
			  ((mplus simp) 1 ((mexpt simp) 5 ((rat simp) 1 2)))))
	       x))
	  ((eq (caar x) 'mexpt)
	   (if (equal (cadr x) '$%e)
	       (*fpexp ($bfloat (caddr x)))
	       (exptbigfloat ($bfloat (cadr x)) (caddr x))))
	  ((eq (caar x) 'mncexpt)
	   (list '(mncexpt) ($bfloat (cadr x)) (caddr x)))
	  ((eq (caar x) 'rat)
	   (ratbigfloat (cdr x)))
	  ((setq y (safe-get (caar x) 'floatprog))
	   (funcall y (mapcar #'$bfloat (cdr x))))
	  ((or (trigp (caar x)) (arcp (caar x)) (eq (caar x) '$entier))
	   (setq y ($bfloat (cadr x)))
	   (if ($bfloatp y)
	       (cond ((eq (caar x) '$entier) ($entier y))
		     ((arcp (caar x))
		      (setq y ($bfloat (logarc (caar x) y)))
		      (if (free y '$%i)
			  y (let ($ratprint) (fparcsimp ($rectform y)))))
		     ((member (caar x) '(%cot %sec %csc) :test #'eq)
		      (invertbigfloat
		       ($bfloat (list (ncons (safe-get (caar x) 'recip)) y))))
		     (t ($bfloat (exponentialize (caar x) y))))
	       (subst0 (list (ncons (caar x)) y) x)))
	  (t (recur-apply #'$bfloat x)))))

(defprop mplus addbigfloat floatprog)
(defprop mtimes timesbigfloat floatprog)
(defprop %sin sinbigfloat floatprog)
(defprop %cos cosbigfloat floatprog)
(defprop rat ratbigfloat floatprog)
(defprop %atan atanbigfloat floatprog)
(defprop %tan tanbigfloat floatprog)
(defprop %log logbigfloat floatprog)
(defprop mabs mabsbigfloat floatprog)

(defmfun addbigfloat (h)
  (prog (fans tst r nfans)
     (setq fans (setq tst bigfloatzero) nfans 0)
     (do ((l h (cdr l)))
	 ((null l))
       (cond ((setq r (bigfloatp (car l)))
	      (setq fans (bcons (fpplus (cdr r) (cdr fans)))))
	     (t (setq nfans (list '(mplus) (car l) nfans)))))
     (return (cond ((equal nfans 0) fans)
		   ((equal fans tst) nfans)
		   (t (simplify (list '(mplus) fans nfans)))))))


;; convert a rat to a float. 1/2015

;;  for n/d rational number   divide(n*10^fpprec,d) to get [s,r] such that  n/d= s/10^fprec+r.
;;  e.g.  321/4 .. fpprec=10, we get                         [802500000000, 0]
;; so 321/4  is  802500000000 x 10^-10 or 80.25.
;; better would be 8025 x10^-2 .
;;trimming trailing decimal zeros could be added..


(defun ratbigfloat(r)  ; r is ((rat) n d)
  (let* ((n (car r))
	 (d (cadr r))
	 (expon (+ (flatsize10 d)fpprec))
	 (shifter (expt 10 expon )))
    ;; add expon zeros to the numerator, divide by d (actually, truncate) then mult  by 10^ expon
    (list `(bigfloat simp ,fpprec) 
	(truncate (* n shifter) d)
	  (- expon))))
      

;; Returns a list of a significand  (ie. mantissa) and an exponent.
(defun intofp (l)
  (cond ((not (atom l)) ($bfloat l))
	((floatp l) (floattofp l))
	((equal 0 l) '(0 0))
	((eq l '$%pi) (fppi))
	((eq l '$%e) (fpe))
	((eq l '$%gamma) (fpgamma))

	(t (list l 0))
	))
;; we have to redefine decimalfptrim do fpround function. For now,
(defun fpround(l) l)

(defun fpquotient (a b) ; decimal version just a hack
  (cond ((equal (car b) 0)
	 (merror (intl:gettext "pquotient: attempted quotient by zero.")))
	((equal (car a) 0) '(0 0))
	(t (decimalfptrim (list (truncate (/ (* (expt 10 fpprec)(car a))
					  (car b)))
					   (- (cadr a)(cadr b) fpprec ))))))


(defun fpdifference (a b)
  (fpplus a (fpminus b)))

(defun fpminus (x)
  (if (equal (car x) 0)
      x
      (list (- (car x)) (cadr x))))

(defun fpplus(a b)
  ;; totally without rounding
  (let ((exp(- (cadr a) (cadr b))))

    (cond ((= exp 0)(list (+ (car a)(car b)) (cadr a)))
	  ((> exp 0)	
	   (list (+ (* (expt 10 exp) (car a))
		    (car b))
			   (cadr b)))
	  
	  (t  (list (+ (* (car b)(expt 10 (- exp))) 
		       (car a))
		    (cadr a))))))


(defun fptimes* (a b)
  ;; totally without rounding
  (list (* (car a)(car b))
	(+ (cadr a)(cadr b))))

;; Don't use the symbol BASE since it is SPECIAL.

(defun fpintexpt (int nn fixprec)	;INT is integer
  (setq fixprec (truncate fixprec (1- (integer-length int)))) ;NN is pos
  (let ((bas (intofp (expt int (min nn fixprec)))))
    (if (> nn fixprec)
	(fptimes* (intofp (expt int (rem nn fixprec)))
		  (fpexpt bas (quotient nn fixprec)))
	bas)))

;; NN is positive or negative integer
(defun fpone()'(1 0))



(defun timesbigfloat (h)
  (prog (fans r nfans)
     (setq fans (bcons (fpone)) nfans 1)
     (do ((l h (cdr l)))
	 ((null l))
       (if (setq r (bigfloatp (car l)))
	   (setq fans (bcons (fptimes* (cdr r) (cdr fans))))
	   (setq nfans (list '(mtimes) (car l) nfans))))
     (return (if (equal nfans 1)
		 fans
		 (simplify (list '(mtimes) fans nfans))))))

(defun invertbigfloat (a)
  (bcons (fpquotient '(1 0)(cdr a))))


(defun *fpexp (a)
  (fpend (let ((fpprec (+ 2. fpprec))) ; was 8
           (if ($bfloatp a)
               (fpexp (cdr (bigfloatp a)))
	     (list '(mexpt) '$%e a)))))


(defun *fpsin (a fl)
  (fpend (let ((fpprec (+ 2. fpprec))) ;was 8
	   (cond (($bfloatp a) (fpsin (cdr ($bfloat a)) fl))
		 (fl (list '(%sin) a))
		 (t (list '(%cos) a))))))

;;; taken out of nparse.lisp. much shortened

(defun make-number (data)
  (setq data (nreverse data))
  ;; Maxima really wants to read in any number as a flonum
  ;; (except when we have a bigfloat, of course!).  So convert exponent
  ;; markers to the flonum-exponent-marker.
  (let ((marker (car (nth 3 data))))
    (unless (eql marker flonum-exponent-marker)
      (when (member marker '(#\E #\F #\S #\D  #+cmu #\W))
        (setf (nth 3 data) (list flonum-exponent-marker)))))
  (if (not (member (car (nth 3 data)) '(#\B #\L) :test #'equal))
      (readlist (apply #'append data))
      (let*
	   ((*read-base* 10.)
	    (int-part (readlist (or (first data) '(#\0))))
	    (frac-part (readlist (or (third data) '(#\0))))
	    (frac-len (length (third data)))
	    (exp-sign (first (fifth data)))
	    (exp (readlist (sixth data))))
      	
	(let* ((bigmant
		  (+ frac-part (* int-part (expt 10 frac-len))))
		 (bigexp (- (if (char= exp-sign #\-)
					 (- exp)
			      exp) frac-len))		 )

	  (bcons (list bigmant bigexp)))
	
	)))

(defun fpe1 ()
  (let ((e (compe (+ fpprec 2))))               ;; compute additional digits 
    (decimalfptrim (list (car e) (cadr e))) ))
;; need the following stuff from float  for $%e, $%pi etc.  Here's the %e part..

(let ((table (make-hash-table)))
  (defun fpe ()
    (let ((value (gethash fpprec table)))
      (if value
	  value
	  (setf (gethash fpprec table) (cdr (fpe1))))))
  (defun fpe-table ()
    table)
  (defun clear_fpe_table ()
    (clrhash table)))
(clear_fpe_table)			; if it is there

(defun fpexp1 (x)
  (let ((term (fpone)) (ans (fpone)) 
	(oans '(0 0)))
     (do ((n 1 (1+ n)))
	 ((equal ans oans) ans)
       (setq term (decimalfptrim(fpquotient (fptimes* x term) (intofp n))))
       (setq oans ans)
       (setq ans (fpplus ans term))
       (format t "~% exp oans= ~s, ans= ~s" oans ans))
    ))

(eval-when
    #+gcl (load eval)
    #-gcl (:load-toplevel :execute)
    (fpprec1 nil $fpprec))		; Set up user's precision