;; -*- Mode:Common-Lisp;Package:mma; Base:10 -*-
(in-package :mma)
;;; Common Lisp Bigfloat Package, Part I: Basic arithmetic only.

;;  copyright (c) 1990    Richard Fateman, UC Berkeley
;;  copyright (c) 2011    Richard Fateman


;;; Bigfloats are stored as structures: precision, fraction and exponent are
;;; integers, with an implied binary point to the left of the fraction. 
;;; The fraction is normalized (except for zero).

;;; We use the common lisp data structure type discrimination mechanism
;;; for tagging, printing, etc. bigfloats.

;;; Another file, bfelem provides sin/cos/log/exp of bigfloats.
;;; bfappend.lisp makes it all accessible from MockMMA.

;;; much of this file and bfelem  (and more) can be done more efficiently
;;; by the mpfr library, should that be available in the host lisp.
;;; see http;//www.cs.berkeley.edu/~fateman/generic for examples.

;;; 
(provide 'bf)

(defvar bigfloat-bin-prec 10 
  "The number of bits in the fraction of a new bigfloat")

(defvar bigfloat-chop nil
"If T makes bigfloat arithmetic faster and less accurate by disabling rounding")  

(defvar bigfloat-print-trunc t 
  "If T means drop trailing 0's in the fraction when printing bigfloats")

(defvar bigfloat-printprec 0
  "If >0, then restrict printing to that many digits when printing bigfloats")

(proclaim 
  '(special bigfloatone bigfloatzero))

;; declare the structure of bigfloats

(defstruct (bigfloat 
	    (:constructor bcons2
			  (fraction exponent
				    &optional (precision bigfloat-bin-prec)))
	    (:print-function bigfloatprintfunction))
  (fraction  0                  :type integer :read-only t)
  (exponent  0                  :type integer :read-only t)
  (precision bigfloat-bin-prec  :type integer :read-only t)
  )

;; test for equality of two bigfloat

(defun bigfloat-eql(x y)
  (and (eql (bigfloat-fraction x)(bigfloat-fraction y))
       (eql (bigfloat-exponent x)(bigfloat-exponent y))))

;; how many bits does it take to represent x decimal digits?
(defun bits-in-decimal(x)
  (ceiling (* x #.(log 10 2))))

;; initialization routines
;; Any time bigfloat precision is changed, one of these routine should
;; be called:  either bigfloat-init-dec or bigfloat-init-bin.

;; bigfloat-init-dec is
;; given: a certain number of decimal digits to be carried (approximately)

(defun bigfloat-init-dec (q)
  (bigfloat-init-bin (bits-in-decimal q)))

;; bigfloat-init-bin is
;; given: a certain number of bits to be carried in the fraction

(defun bigfloat-init-bin (q)  
  (setq bigfloat-bin-prec q
	bigfloatone (intofp 1) ;;initialize one
	bigfloatzero (bcons2 0 0)) ;; initialize zero
  q)

(defun intlen(x)(integer-length (abs x)))
  
;; into-bigfloat is called by the reader to make a bigfloat.
;; The two arguments are integers representing frac * 10^exp.
;; The input routines should convert 1.23*10^5 to 123 * 10^(-2),
;; and then call into-bigfloat on 123 and -2.
;; the return value is an accurately converted number.

(defun into-bigfloat (frac exp)
  ;; log[2](10) is 3.321928 ...
   ;; jack up precision if it is needed for good conversion
   (let ((bigfloat-bin-prec 
	  (max bigfloat-bin-prec
	       (+ (ceiling 3.321193 (intlen exp))
		  (intlen frac)))))
     ;; compute frac*10^exp
     (bigfloat-* (intofp frac) (bigfloat-expt(intofp 10)exp))))

(defun bigfloat-convert(x) ;; takes whatever x is, and returns a bigfloat.
  (typecase x
	     (bigfloat
	      (cond((= bigfloat-bin-prec (bigfloat-precision x)) x)
	       (t (bfnormal 
		   (bigfloat-fraction x)
		   (bigfloat-exponent x)))))
	     (integer (intofp x)) ;bignum or fixnum
	     	 ;; this is the common-lisp rational number type
	     (rational
	      (bigfloat-/ ;; can't lose on under/overflow
	       (intofp (numerator x))
	       (intofp (denominator x))))
	     (float ;single, double or other float
	      (bigfloat-convert (rational x)))
	     (otherwise x))) ;; if it can't be converted, return x? probably not useful

	     
;; integer to bigfloat conversion

(defun intofp(l)(if (equal 0 l) bigfloatzero (bfnormal l (intlen l))))

;; bfnormal "normalizes" a floating point fraction.
;; the proposed number has fraction frac (with binary bit at the left)
;; and exponent exp.
;; But frac may have the wrong number of bits.  To be
;; normalized, (intlen frac) must be = to bigfloat-bin-prec
;; unless frac =0.

(defun bfnormal(frac exp &optional (il (intlen frac)))
  (cond ((= il bigfloat-bin-prec) (bcons2 frac exp))
	(bigfloat-chop (bfnormal-chop frac exp il))
	(t(bfnormal-round frac exp il))))

(defun bfnormal-chop (frac exp il);
  (if (= frac 0) bigfloatzero
    (bcons2 (ash frac (- bigfloat-bin-prec il))
	    exp)))

(defun bfnormal-round (frac exp il);
  (if (= frac 0) bigfloatzero
    (let
	((s (- bigfloat-bin-prec  il)))
      (cond ((> s 0) ;; not enough bits, so shift left what you have
	     (bcons2 (ash frac s) exp))
	    ;; s= 0 taken care of in bfnormal.
	    (t  ;; that is, (< s 0) ;; too many bits. Must round.
	     (let ((r (round frac (ash 1 (- s)))))
	       (if (= (intlen r) bigfloat-bin-prec)
		   (bcons2 r exp)  ;; usual case
		 ;; Take care of rare case when rounding up bumps the length.
		 ;; Analogy in decimal: round 999 x 10^0 to 2 digits.
		 ;; You get 100 x 10^3.
		 (bcons2 (ash r -1) (1+ exp)))))))))

;; add two bigfloats of any precision; converts each to global precision.
;; c = a+b.

(defun bigfloat-+ (a b)
  (setq a (bigfloat-convert a))
  (setq b (bigfloat-convert b))
  (let* ((fa (bigfloat-fraction a))
	 (fb (bigfloat-fraction b))
	 (ea (bigfloat-exponent a))
	 (eb (bigfloat-exponent b))
	 fc
	 (dif (- ea eb)))
    (cond ((> dif     bigfloat-bin-prec) a) ;b is insignificant
	  ((> (- dif) bigfloat-bin-prec) b) ;a is insignificant
	  ((>= dif 0)
	   (setq fc (+ (ash fa dif) fb))
	   (bfnormal 
	    fc
	    (+ eb (- (intlen fc) bigfloat-bin-prec))))
	  (t (setq fc (+ (ash fb (- dif)) fa))
	     (bfnormal 
	      fc
	      (+ ea (- (intlen fc) bigfloat-bin-prec)))))))

;; bigfloat-* multiplies two bigfloats of arbitrary (perhaps different)
;; precisions, and return a bigfloat of global (bigfloat-bin-prec)
;; precision.

(defun bigfloat-* (a b) 
  (let* ((f (* (bigfloat-fraction a)(bigfloat-fraction b)))
	 (il (intlen f))
	 (exp   (+ (bigfloat-exponent a)
		   (bigfloat-exponent b)
		   (- (bigfloat-precision a))
		   (- (bigfloat-precision b))
		   il)))
    (bfnormal f exp il)))

;; bigfloat-/ divides two bigfloats of arbitrary (perhaps different)
;; precisions, and return a bigfloat of global (bigfloat-bin-prec)
;; precision.

(defun bigfloat-/ (a b) 
  ;; division by zero will be noticed in the "round" below
  (let* ((f (round (ash (bigfloat-fraction a)
			(+ bigfloat-bin-prec (bigfloat-precision b)))
		   (bigfloat-fraction b)))
	 (il (intlen f))
	 (exp   (+ (bigfloat-exponent a)
		   (- (bigfloat-exponent b))
		   (- (bigfloat-precision a))
		   (- bigfloat-bin-prec)
		   il)))
    (bfnormal f exp il)))

;; compute 1/x for bigfloat x

(defun bigfloat-inv(x)  (bigfloat-/ (bfone) x))

;; coerce-bigfloat 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-bigfloat(x typ) 
  (cond ;; convert bigfloat to bigfloat of global precision
        ((eq typ 'bigfloat) 
	 (bigfloat-convert x))
	((eq typ 'ratio)
	 ;; convert a bigfloat to a rational fraction. This is exact in 
	 ;; the sense that for any represented bigfloat x in bigfloat-bin-prec,
	 ;; (bigfloat-convert( bigfloattorat x)) = x
	 (* (bigfloat-fraction x)
	    (expt 2 (- (bigfloat-exponent x) (bigfloat-precision x)))))
	((eq typ 'double-float)
	 ;; convert a bigfloat into a floating point double when it is
	 ;; simple to do.  Won't work for very high precision, since
	 ;; the fraction part will overflow.
	 (* (coerce (bigfloat-fraction x) 'double-float)
	    (expt 2.0d0 (-  (bigfloat-exponent x)(bigfloat-precision x)))))
	((member typ '(float single-float) :test #'eq)
	 (coerce (coerce-bigfloat x 'double-float) 'single-float))
	((member typ '(integer bignum fixnum))
	 (round (coerce-bigfloat x 'ratio)))
	(t (error "Can't coerce ~s to type ~s~%" x typ))))

;; how many decimal digits in an integer x?
;; slow, but accurate.  Used in i/o conversions.

(defun decimalsin(x)  (length (format nil "~s" x)))

;; various useful predicates

(defun bigfloat-> (a b) (bigfloat-posp (bigfloat-diff  a b))) 
(defun bigfloat-< (a b) (bigfloat-posp (bigfloat-diff  b a))) 
(defun bigfloat-posp (x)(> (bigfloat-fraction x) 0))
(defun bigfloat-zerop(x)(= (bigfloat-fraction x) 0))

(defun bigfloat-min(&rest args)
  (let (min)
    (cond ((null args) (error "bigfloat-min with no arguments")))
    (setq min (car args))
    (do ((a (cdr args)  (cdr a)))
	((null a) min)
	(cond ((bigfloat-< (car a) min)
	       (setq min (car a)))))))

(defun bigfloat-max(&rest args)
  (let (max)
    (cond ((null args) (error "bigfloat-max with no arguments")))
    (setq max (car args))
    (do ((a (cdr args)  (cdr a)))
	((null a) max)
	(cond ((bigfloat-> (car a) max)
	       (setq max (car a)))))))

;; return an existing 1 or make a new one.

(defun bfone nil 
       (cond
	((= bigfloat-bin-prec (bigfloat-precision bigfloatone)) 
	 bigfloatone)
	(t (setq bigfloatone (intofp 1)))))

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

(defun bigfloat-- (x &optional (y nil))
  (cond (y (bigfloat-+ x (bigfloat-- y)));; compute x-y
	((bigfloat-zerop x) x)
	(t (bcons2 
	    (- (bigfloat-fraction x)) (bigfloat-exponent 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 bigfloat-expt (p nn) 
  (cond ((eql nn 0) (bfone))
	((eql nn 1) p)
	((< nn 0) (bigfloat-inv (bigfloat-expt p (- nn))))
	(t 
	 (do ((n (floor nn 2)(floor n 2))
	      (s (cond ((oddp nn) p) (t (bfone)))))
	     ((zerop n) s)
	     (setq p (bigfloat-* p p))
	     (and (oddp n) (setq s (bigfloat-* s p))) )
	   )))

;; compute square-root of bigfloat x

(defun bigfloat-sqrt(x)(bigfloat-root x 2))

;; compute a^(1/n)  see Fitch, SIGSAM Bull Nov 74

(defun bigfloat-root (a n) 
  (let* ((ofprec bigfloat-bin-prec) 
	 (bigfloat-bin-prec (+ bigfloat-bin-prec 2))  ;assumes a>0 n>=2
	 (bk (bigfloat-expt (intofp 2)
			    (1+ (truncate (bigfloat-exponent
					   (setq a (bigfloat-convert a)))
					  n)))))
    (do ((x bk
	    (bigfloat-diff
	     x (setq bk (bigfloat-/ (bigfloat-diff
				     x 
				     (bigfloat-/ a (bigfloat-expt x n1)))
				    n))))
	 (n1 (1- n))
	 (n (intofp n)))
	((or (bigfloat-zerop bk)
	     (> (- (bigfloat-exponent x)
		   (bigfloat-exponent bk)) ofprec)) (setq a x))))
  (bigfloat-convert a)) ;return to previous precision
 
(defun bigfloat-abs (x) 
       (cond ((>= (bigfloat-fraction x) 0) x)
	     (t (bcons2 (- (bigfloat-fraction x)) (bigfloat-exponent x))))) 

;; integer part of a bigfloat

(defun bigfloat-intpart (f) 
 (let ((m (- (bigfloat-precision f) (bigfloat-exponent f))))
   (if (> m 0) (truncate (bigfloat-fraction f) (ash 1 m))
     (* (bigfloat-fraction f) (ash 1 (- m))))))


;;set up user's precision

;;(bigfloat-init-dec 3);; 3 decimal digits  for debugging

;; some tests, some utilities

(defun bc (x &optional
	     (bigfloat-bin-prec bigfloat-bin-prec )) (bigfloat-convert x))

;; this definition helps for tracing functions so you can read everything..
;; at least while you are using bignums that are exactly representable

(defvar bigfloat-print-trace nil "t if you want to see internals of bigfloats")

(defun bigfloatprintfunction (x s pl)  ;pl, print-level, is not used.
  (cond (bigfloat-print-trace
	 (format s "[~s*2^~s]=~a"
		 (bigfloat-fraction x)
		 (bigfloat-exponent x)
;;		 (coerce-bigfloat x 'double-float) ; test
		 (bigfloat-print x nil)
		 ))
	(t (bigfloat-print x s))))

;;(bigfloat-print x) normally prints the bigfloat to the standard output.
;;  it looks like d.dd*10^ee  where the number of digits, d, is specified
;;  by the width w, or if w is not provided, by the number of decimal
;;  digits assumed plausible by the binary precision.
;;  ee, the exponent, is computed as appropriate.

;;(bigfloat-print x stream) prints to a stream, and if stream=nil,
;; it returns a string.

;;(bigfloat-print x stream w) prints to a stream (if stream=t, standard
;; output; and also uses a width for the fraction part of w digits.

;; The total field width will be at most  w + 4 + width of exponent.
;; optionally there will be a negative sign taking one space, also.

;; this version of bigfloat-print prints d.ddd..d *10^n
;; except if n=0, when it prints d.ddd..d  (no *10^n).  Also
;; trailing zeros in the fraction are dropped off if bigfloat-print-trunc
;; is t, and 

(defun bigfloat-print (x &optional (stream t) (w bigfloat-printprec)&aux sign)
  (flet 
   ;; utility functions
   ((t0(s); change "12300" to "123"
       (if bigfloat-print-trunc (string-right-trim "0" s) s))
    (putindot(s) ; change "123" to "1.23"
	     (concatenate 'string (subseq s 0 1) "." (subseq s 1))))
   (cond((bigfloat-zerop x) (format stream "0.0"))
	(t (multiple-value-bind (frac exp)
				(bf-format x w)
				(setq sign (< frac 0))
				(setq frac (format nil "~s" (abs frac)))
				(format stream "~:[~;-~]~a~:[*10^~s~;~]" 
					sign
					(t0(putindot frac))
					(zerop 
					 (setq exp(+ (length frac) -1 exp)))
					exp))))))


;; bf-format returns a possible pair of decimal fraction and exponent
;; values for printing x.
;; In particular, two integer values f and e are returned, and the
;; number x is closely approximated by f x 10^e.

(defun bf-format( x w)
  (cond ((bigfloat-zerop x) (values 0 0))
	(t
	 (let*((f (abs(bigfloat-fraction x)))
	       (e (bigfloat-exponent x))
	       (p (bigfloat-precision x))
	       (bigfloat-bin-prec p)   ;use the precision of the number
	       (width (if (> w 0) w (truncate (* bigfloat-bin-prec 0.30103))))
	       ;; to convert .f  x 2^e to decimal
	       ;; compute the log[10] of the number, approximately
	       (d (floor
		   (+
		   ;; log of the exponent
		    (* (+ e (- p) (intlen f))
		      #.(/ 1.0d0(log 10.0d0 2)))
		   ;; log of the fraction
		   ;; we can do this in double-float since it is
		   ;; unreasonable to have an exponent exceeding +/- 10^307.
		   ;; Numbers of the size (2 ^(10^307)) won't
		   ;; print.
		   (log (coerce (/ f (expt 2 (intlen  f)))
				'double-float)
			10))))
	       (s (if (>= d 0)
		      (bigfloat-/ x ;;note: x still has its sign
			      ;;this is cheap 
			      ;;(bigfloat-expt (intofp 10) d)
			      (intofp (expt 10 d));; this is accurate
			      )
		    (bigfloat-* x (intofp (expt 10 (- d))))))
	       ;; now  1 < |s| <= 10 ? , x = s * 10^d
	       (left(round (/(* (expt 10 (1- width))(bigfloat-fraction s))
				    (expt 2 (- (bigfloat-precision s)
					       (bigfloat-exponent s)))))) )
	   (cond ((> (decimalsin left) width)
		  (setq d (1+ d) left (round left 10))))	
	   ;; this could be set up as a string to return..
	   ;; The two items of interest that could be returned are
	   ;; the fraction and exponent

	   (values left (- d (1- width)))))))