;; -*- Mode:Common-Lisp;Package:mma; Base:10 -*-
;; Mathematica(tm)-like evaluator / bigfloat appendix
;; copyright (c) 1990 Richard J. Fateman; pieces by Tak Yan, Derek Lai
;;
(proclaim
 '(special Precision BigfloatChop BigfloatPrintTruncate 
	   BigfloatPrintPrecision funnyvars))

(require 'math-eval "eval")
(require "bfelem")

;; There is an erroneous treatment of the scope of these variables.
;; We use global bindings when other, local, bindings are probably
;; appropriate.  This can be fixed later, when
;; we implement Set for local variables.  Right now we 
;; are using local bindings only in pattern matching.
;; -- not done as of 2/22/91 (RJF)

(setf (gethash 'Precision funnyvars)
      #'(lambda(x)(bigfloat-init (the integer x))
	  ;; this next line is like Set, except ignoring funnyvars
	  (setf (gethash 'Precision  (gethash 'Precision symtab)) x)))

;; BigfloatChop=0 is off  or   BigfloatChop=<anything else> is the choice

(setf (gethash 'BigfloatChop funnyvars)
      #'(lambda(x)(setq bigfloat-chop (cond((=(the integer x) 0)nil)
					   (t nil)))
	  (setf (gethash 'BigfloatChop  
			 (gethash 'BigfloatChop symtab)) x)))

(setf (gethash 'BigfloatPrintTruncate funnyvars)
      #'(lambda(x)(setq bigfloat-print-trunc (cond((=(the integer x) 0)nil)
						  (t t)))
	  (setf (gethash 'BigfloatPrintTruncate
			 (gethash 'BigfloatPrintTruncate symtab)) x) ))

(setf (gethash 'BigfloatPrintPrecision funnyvars)
      #'(lambda(x)(setq bigfloat-chop (the integer x))
	  (setf (gethash 'BigfloatPrintPrecision 
			 (gethash 'BigfloatPrintPrecision symtab)) x)))

(bigfloat-init 3)

(mapc #'chash '(Precision BigfloatChop BigfloatPrintTruncate BigfloatPrintPrecision))

(Set 'Precision 3)
(Set 'BigfloatChop 0) ;off
(Set 'BigfloatPrintTruncate 1) ;on
(Set 'BigfloatPrintPrecision 3) ;low precision

(require 'math-eval "eval")

#+ignore
(setf (gethash 'Precision funnyvars) #'(lambda(h)(format t "Precision set to ~s" h)))
  
;;; This appendix assumes you've read in eval.lisp already,
;;; These are examples of how to write Lisp functions for lmath evaluation.
;;; Look especially at Sin, for a simple 1-argument function.

;;; substitute an appropriate operator for the built-in
;;;      operators like * + sin cos etc.

(defun Real(a b)
  (bigfloat-+  (intofp a) (bigfloat-/ (intofp b)
				      (intofp(expt 10 (decimalsin b))))))

(defun decimalsin(x)
  (ceiling (* (integer-length x) 0.30102999566398114d0)))

#+ignore (defun Times (&rest x)(cond ((null x)bigfloatone)
			    ((null (cdr x))(car x))
			    (t (let ((s (apply #'Times (cdr x)))
				     (r (car x)))
				 (cond ((and (bigfloat-p r)
					     (bigfloat-p s))
					(bigfloat-*  r s))
				       (t (ulist 'Times  r s)))))))


;; This is the place where we have to determine the nature of 
;; result for the product of any member of the set of types below
;; by any other member of the set
;;  { 
;; (a) CL number
;;    subcases include integer, floats (single, double, extended)
;;                     ratio, 
;;                     complex (integer), complex (float) etc.
;; (b) bigfloat (our defstruct).
;;      subcase: complex bigfloat
;; (c) rational form  (that is, our defstruct, rat)
;;     plausible subcases might be various rings or fields such as
;;         Z[x1,...,xn], Z(x1,...xn), Q[x1,...,xn], Z[x1,...x(n-1)] (xn)
;;         algebraic number or algebraic function fields, etc.
;; (d) Element of a finite field (e.g. 10 mod 23)
;; (e) Interval (endpoints any manifest scalar CL number or bigfloat)
;;     subcases interior or exterior, open or closed.
;; (f) matrix 
;;     subcases: sparse, dense, different element types
;;     vector (cross or dot or by scalar)
;; (g) formal Sum (e.g. Sum[f[#]&,{#,low,high}],

;; (h) Series object (e.g. 1+x + O[x]^2)
;; (i) Operator (e.g. D[x]* x^2 => 2*x, or D*D [x] for 2nd derivative)
;; (j) Sets
;; (h) other symbolic form (e.g. z, Pi, 3*x, or x+y^2)}

;; Rules:  (Note: there seem to be at least 11x11 cases, and maybe far more
;;          if you consider subcases and alternatives not listed.
;;          This is why Scratchpad II seems like a good idea.)

;; (1) The product of any two types in the same category is treated
;;     by a program that depends on the type.  For example, the product
;;     of any two CL types (elements in (a)) is computed by using 
;;     the CL coercions.
;;     (This is not so good because an overflowing flt.pt. product
;;     can still be represented by a large ratio. Also, we might differ
;;     with CL on its treatment of 1/0 and 0/0.)
;;     The product of (for example) two formal sums can be dealt with by
;;     a program that need not know about other stuff.

;; (2) The product of any element in (b) by an element in (a) is in (b).
;;     (Bigfloats rule over other CL numbers)

;; (3) The product of any element in (c) by an element of (a) or (b)
;;     requires either adoption of a new coefficient domain for the
;;     rationals or conversions from rational to "symbol" (h).
;;     Do we want /rational/  0.7500 x^2  ? For (c) by elements of (d)-(h),
;;     encapsulate them as rational kernels.

;; (4) Finite field arithmetic isn't handled at all because we aren't
;;     representing those objects yet.


;; (n) (Stopgap) elements in (d)-(h) are treated indistinguishably.

;; Assume that, after collection, we have at most one element in each
;; category.  We then combine them if appropriate.  That's pretty bad.

;; Another strategy:  if you have disparate rational forms, disparate
;; bigfloat forms, etc., etc. then you proclaim at any one time
;; a preferred form.  e.g. order of variables, type of coefficient,
;; bigfloat precision.  Convert to that global form.  This is already
;; the strategy in bigfloats.  Can it be generalized?
 
(defun Times (&rest x &aux (nums 1) oths)
  (dolist (h x   ;; iterate with h running over all args
	     ;; two sub-results are formed. The numerical coefficient is
	     ;; accumulated in  nums.  The symbolic stuff, if any, is
	     ;; accumulated in  oths.
	     
	   ;;resultform is computed as the product of nums and oths
	   (cond((eql 1 nums)
		 (if (null oths) 1 (ucons 'Times (uniq(nreverse oths)))))
		((null oths) nums)
		(t (ucons 'Times (ucons nums (uniq(nreverse oths)))))))
	  
	 ;; body
	 (typecase h
		   ((and number (not complex))  ;; would need bigfloat-complex
		    (setq nums (cond ((bigfloat-p nums)
				      (bigfloat-* (bigfloat-convert h)
						  nums))
					    (t(* nums h))))) ;; collect CL numbers
		   (bigfloat (setq nums
				   (if (eql nums 1) h
				       (bigfloat-* h (bigfloat-convert nums)))))
		   (rat 	   ; if you find a rat, break out, apply Rule 3!
		    (return-from Times
				 (reduce-rat #'rat* 
					     (into-rat (car x))
					     (cdr x))))
		   (t (push h oths)))))

;; reprogrammed along the model for Times
(defun Plus  (&rest x &aux (nums 0) oths)
  (dolist (h x   ;; iterate with h running over all args
	     ;; two sub-results are formed. The numerical coefficient is
	     ;; accumulated in  nums.  The symbolic stuff, if any, is
	     ;; accumulated in  oths.
	     
	   ;;resultform is computed as the product of nums and oths
	   (cond((eql 0 nums)
		 (if (null oths) 0 (ucons 'Plus (uniq(nreverse oths)))))
		((null oths) nums)
		(t (ucons 'Plus (ucons nums (uniq(nreverse oths)))))))
	  
	 ;; body
	 (typecase h
		   ((and number (not complex))
		    (setq nums (cond ((bigfloat-p nums)
				      (bigfloat-+ (bigfloat-convert h)
						  nums))
					    (t(+ nums h))))) ;; collect CL numbers
		   (bigfloat (setq nums
				   (if (eql nums 0) h
				       (bigfloat-+ h (bigfloat-convert nums)))))
		   (rat 	   ; if you find a rat, break out, apply Rule 3!
		    (return-from Plus
				 (reduce-rat #'rat+ 
					     (into-rat (car x))
					     (cdr x))))
		   (t (push h oths)))))

  ;;; this will be a bigfloat if b is bigfloat, e is integer
(defun Power (b e)

		    (cond((and (bigfloat-p b)(integerp e))
			  (bigfloat-expt b e))
			 ;; what if b and e are both bigfloats
			 ((and (numberp b)(numberp e))(expt b e))
			 (t (ulist 'Power b e))))

;; these programs deal fully only with non-complex bigfloats.
;; this should probably be fixed at some point.  The kind of
;; changes are hinted at in the definition of Log.

(defun Log(x)
  (typecase x
	    (bigfloat (if (bigfloat-posp x)
			  (bigfloat-log x) ;; x must positive 
			(ulist 'Log x))) ;; really should produce complex
	    (number (log x))
	    ;; here we could have some checks on special arguments
	    (t (ulist 'Log x))))


(defun Sin(x)
  (typecase x
	    (bigfloat (bigfloat-sin x))
	    ;; here, should check for 0 to make sure it's exact.
	    
	    (number (sin x))
	    ;; here, check for complex/bigfloat
	    ;; here, should check for familiar multiples of Pi,
	    ;; here, could also check sin(ArcSin), sin(ArcCos), 
	    ;; ArcTan special formulas
	    ;; perhaps Sin[x+y]-> ...
	    ;; perhaps Sin[I*x] -> Sinh[x]
	    ;; distribute over compound objects e.g. Sin[Matrix...]
	    ;; 
	    (t (ulist 'Sin x))))
(defun Cos(x)
  (typecase x
	    (bigfloat (bigfloat-cos x))
	    (number (cos x))
	    (t (ulist 'Cos x))))

(defun Tan(x)
  (typecase x
	    (bigfloat (bigfloat-tan x))
	    (number (tan x))
	    (t (ulist 'Tan x))))
(defun Exp(x)
  (typecase x
	    (bigfloat (bigfloat-exp x))
	    (number (exp x))
	    (t (ulist 'Exp  x))))



;; Pi, to precision  Precision.

(defun Pi()(fppi))