Some random thoughts on future features

>>>>> "Barton" == Barton Willis <willisb@unk.edu> writes:

    Barton> When I wrote the orthopoly numerical stuff, I briefly thought about using 
    Barton> CLOS to do double-float,  bigfloat, and complex floating point math. 

Here is a very, very preliminary version that is very gross.  It
barely does anything, but it's enough to compute jacobi_sn for
bigfloat args to bigfloat accuracy:

(%i1) fpprec:50;

(%o1) 				      50
(%i2) jacobi_sn(x,1);

(%o2) 				    tanh(x)
(%i3) tanh(0.5b0);

(%o3) 	    4.6211715726000975850231848364367254873028928033011b-1
(%i4) jacobi_sn(0.5b0,0.99999999999999999999999999999b0),bfloat;

(%o4) 	    4.6211715726000975850231848364384478195723189546368b-1


At this point, I'm kind of soliciting comments on the names of the
classes and the names of the methods.

It's going to take some time to get all of this in place.  To be
really useful, it has to basically provide support for all of the
basic math operations and special functions that Lisp provides.

Ray

;; Bigfloat and complex bigfloat classes
;; A bigfloat is basically a maxima bigfloat and contains the mantissa
;; and exponent part.

(defclass bigfloat ()
  ((value :reader bfloat
	  :initarg :value)))

(defmethod make-bigfloat ((x real))
  (make-instance 'bigfloat :value (intofp x)))

(defmethod print-object ((x bigfloat) (stream t))
  (print-unreadable-object (x stream :type t :identity t)
    (format stream "~A" (bfloat x))))

(defmethod make-bigfloat ((x list))
  (make-instance 'bigfloat :value x))

(defclass complex-bfloat ()
  ((real :reader real-bfloat :initarg :real)
   (imag :reader imag-bfloat :initarg :imag)))

(defmethod make-complex-bigfloat ((x real) (y real))
  (make-instance 'complex-bfloat
		 :real (intofp x)
		 :imag (intofp y)))

(defmethod bf+ ((x real) (y real))
  (+ x y))

(defmethod bf- ((x real) (y real))
  (- x y))

(defmethod bf* ((x real) (y real))
  (* x y))

(defmethod bf/ ((x real) (y real))
  (/ x y))

(defmethod bf-sqrt ((x real))
  (sqrt x))

(defmethod bf-expt ((x real) (y integer))
  (expt x y))

(defmethod bf+ ((x bigfloat) (y bigfloat))
  (make-bigfloat (fpplus (bfloat x) (bfloat y))))

(defmethod bf- ((x bigfloat) (y bigfloat))
  (make-bigfloat (fpdifference (bfloat x) (bfloat y))))

(defmethod bf* ((x bigfloat) (y bigfloat))
  (make-bigfloat (fptimes* (bfloat x) (bfloat y))))

(defmethod bf/ ((x bigfloat) (y bigfloat))
  (make-bigfloat (fpquotient (bfloat x) (bfloat y))))

(defmethod bf+ ((x real) (y bigfloat))
  (make-bigfloat (fpplus (intofp x) (bfloat y))))

(defmethod bf- ((x real) (y bigfloat))
  (make-bigfloat (fpdifference (intofp x) (bfloat y))))

(defmethod bf- ((x bigfloat) (y real))
  (make-bigfloat (fpdifference (bfloat x) (intofp y))))

(defmethod bf* ((x real) (y bigfloat))
  (make-bigfloat (fptimes* (intofp x) (bfloat y))))

(defmethod bf/ ((x real) (y bigfloat))
  (make-bigfloat (fpquotient (intofp x) (bfloat y))))

(defmethod bf-sqrt ((x bigfloat))
  (make-bigfloat (fproot `((bigfloat simp ,fpprec) ,@(bfloat x)) 2)))

(defmethod bf-expt ((x bigfloat) (y integer))
  (make-bigfloat (fpexpt (bfloat x) y)))

;; An example of computing jacobi_sn for bigfloats.
;; Only works real args and for 0 < m < 1.

(defprop %jacobi_sn jacobi-sn-bigfloat floatprog)

(defmfun jacobi-sn-bigfloat (args)
  (labels ((descending-transform (u m)
	     (let* ((root-m1 (bf-sqrt (bf- 1 m)))
		    (root-mu (bf/ (bf- 1 root-m1)
				  (bf+ 1 root-m1)))
		    (mu (bf* root-mu root-mu))
		    (v (bf/ u (bf+ 1 root-mu))))
	       (values v mu root-mu)))
	   (sn-descending (u m)
	     (cond ((< (second (bfloat m)) (- fpprec))
		    (make-bigfloat (fpsin1 (bfloat u))))
		   (t
		    (multiple-value-bind (v mu root-mu)
			(descending-transform u m)
		      (let* ((new-sn (sn-descending v mu)))
			(bf/ (bf* (bf+ 1 root-mu)
				  new-sn)
			     (bf+ 1 (bf* root-mu
					 (bf* new-sn new-sn))))))))))
    (bcons (bfloat (sn-descending (make-bigfloat (rest (first args)))
				  (make-bigfloat (rest (second args))))))))