;;; -*- Mode: LISP; Syntax: Common-Lisp; Package: qd; Base: 10 -*-
;; Author: Richard Fateman, Jan, 2006
;; quad-double

;(load "qd.dll")
(defpackage :qd
  (:use :ga :cl)
    (:shadowing-import-from  
   :ga    
   "+" "-" "/" "*" "expt"		;... n-ary arith
   "=" "/=" ">" "<" "<=" ">="		;... n-ary comparisons
   "1-" "1+" "abs" "incf" "decf"
   "min" "max"
   asin log tan atanh acos sin sinh cosh cos sqrt exp atan
   "tocl" "re-intern" 
   )
  )

(in-package :qd)

(eval-when (compile load)
  (declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
  )

(eval-when (compile load eval)
  (ff:def-foreign-type cqd (cl::* :double))
  (ff:def-foreign-type cdd (cl::* :double))

  (ff:def-foreign-call
   (qd_mul  "c_qd_mul")
   ((op1 cqd)(op2 cqd) (target cqd))
   :returning :void 		         	     
   :arg-checking nil :call-direct t)
  
  (ff:def-foreign-call
   (qd_add  "c_qd_add")
   ((op1 cqd)(op2 cqd) (target cqd))
   :returning :void 		         	     
   :arg-checking nil :call-direct t)
 (ff:def-foreign-call
   (qd_sub  "c_qd_sub")
   ((op1 cqd)(op2 cqd) (target cqd))
   :returning :void 		         	     
   :arg-checking nil :call-direct t)
  
   (ff:def-foreign-call
   (qd_div  "c_qd_div")
   ((op1 cqd)(op2 cqd) (target cqd))
   :returning :void 		         	     
   :arg-checking nil :call-direct t)
  )

;; extend the list of operations with other entry points as
;; needed. There are 120 external entry points for qd and dd
;; routines. We can write macros to expand them into lisp
;; ff interfaces, based on number of arguments and some
;; minor string-hacking.

;; An timing-minded programmer might wish to access the
;; raw qd routines than use the overloaded +,*, etc. We
;; expect that the overload/discrimination mechanism may
;; be (percentagewise) more of an issue since the qd routines
;; may be much faster than (say) the arbitrary precision ones.
;; So far (2/3/06) this file does not have the overloaded
;; names put in there.

(defun init-qd() 
  (make-array 4 :element-type 'double-float :initial-element 0.0d0))

;;If we want to tag these guys we can do something like this:

(defstruct aqd :q)
;; Here q is the "guts" of the repr. An array of doubles.
;; The structure "aqd" provides a wrapper or tag so that
;; we can hang methods on it, like print-object, +, * etc.

(defmethod print-object ((a aqd) stream)(format stream "~a"
						 (qd2string a)))

(defun alloc-qd()  ;; provides an empty qd, with a tag so lisp knows it
    (make-aqd :q (init-qd)))

;; The next program converts a qd number to a lisp rational.
;; probably a "more efficient" way to do this; note denom is a power of 2

(defmethod qd2lisp((x aqd))
  (or (excl::nan-p (aref (aqd-q x) 0)) ; if qd contains a NaN, use it.
      (apply #'+  (map 'list #'rational (aqd-q x)))))

;; To convert from an ordinary lisp "real" number one of these should do:

(defmethod lisp2qd((x rational)) ;to encode a lisp rational, divide num/den
  (let ((ans (init-qd))
	(nu (aqd-q (lisp2qd (numerator x))))
	(de (aqd-q (lisp2qd (denominator x)))))
    (qd_div nu de ans)
    (make-aqd :q ans)))

(defmethod lisp2qd((x fixnum))  ;small integers are easy.
  (lisp2qd (coerce x 'double-float)))

(defmethod lisp2qd ((x double-float)) ;so are double-floats
    (let ((ans (init-qd)))
    (setf (aref ans 0) x)
    (make-aqd :q ans)))

(defmethod lisp2qd((x single-float)) ;so are single-floats
 (lisp2qd (coerce x 'double-float)))
  
(defmethod lisp2qd ((x integer)) 
  ;; integers that are shorter than a double-float fraction are easy.
  (if (< (cl::abs x) #.(cl::expt 2 53)) (lisp2qd  (coerce x 'double-float))
    ;; the rest of this code is for when x is a bignum.
    ;; We convert it, 52-bit section by 52-bit section
    (let ((ans (init-qd))
	  (s (signum x))
	  (shifter (aqd-q (lisp2qd 1))) ;initially, shift by 1
	  (newdig 0)
	  (shiftam(aqd-q (lisp2qd #.(expt 2 52))) ))
      
      (if (< s 0)(setf x (- x)))
      (loop while (> x 0) do
	    (setf newdig (logand x #.(cl::1- (expt 2 52))))
	    (setf newdig (aqd-q (lisp2qd newdig))) ; grab some bits
	    (qd_mul shifter newdig newdig) ;newdig now a qd
	    (setf x (ash x -52)) ;remove them from x
	    (qd_add ans newdig ans)
	    (qd_mul shifter shiftam shifter))
      (if (< s 0)(qd_mul ans (lisp2qd -1) ans))
      (make-aqd :q ans))))

;; should be faster ways of shifting and changing sign in the above.
  

;; Converting from qd to lisp to decimal (or other base) string we
;; look at the first n decimal digits of a fraction, also show exponent.
;; We only use base 10 in qd2string.

(defun decimalize(r n &optional (base 10))	; r rational number >=0
  (let ((expon  (if (= r 0) 0 (floor(cl::log (cl::abs r) base) ))))
       (values   (signum r)
	        (round (*(cl::abs r) (expt base (- n expon))))
		(cl::1+ expon))))

;; make a formatted output . Something like this..
;; this is  used by print-object.

(defparameter *qd-digits-to-show* 66)

(defmethod qd2string((x aqd))
  ;; check if x is a NaN
  (if (excl::nan-p (aref (aqd-q x) 0)) "NaN(qd)"
    (multiple-value-bind (s r e h)  ;h is extra variable
	(decimalize (qd2lisp x) *qd-digits-to-show* 10)
      (format nil "~a0.~aQ~s" (if (< s 0)"-" "") 
	      (string-right-trim "0"  
				 (subseq (setf h(format nil "~a" r)) 0  
					(cl::min (length h) *qd-digits-to-show*)) )
				 e))))
    

(defmacro defarithmetic (op pgm)
  (let ((two-arg
	 (intern (concatenate 'string "two-arg-" (symbol-name op))
		 :ga )) )
    `(progn
       ;; new defmethods for qd. .. note order of args different from gmp
       (defmethod ,two-arg ((arg1 aqd) (arg2 aqd))
	 (let* ((r (alloc-qd)) (in (aqd-q r)))
	   (,pgm (aqd-q arg1)(aqd-q arg2) in) r))
       
       (defmethod ,two-arg ((arg1 real) (arg2 aqd))
	 (let* ((r (lisp2qd arg1))(in (aqd-q r)))
	   (,pgm in (aqd-q arg2) in) r))
       
       (defmethod ,two-arg ((arg1 aqd) (arg2 real))
	 (let* ((r (lisp2qd arg2))(in (aqd-q r))) 
	   (,pgm (aqd-q arg1) in in) r))
       
       (compile ',two-arg)
       (compile ',op)
       ',op)))

(defarithmetic + qd_add)
(defarithmetic - qd_sub)
(defarithmetic * qd_mul)
(defarithmetic / qd_div)

;; do analogous stuff for other entry points,

(defmacro r (op);;
  (let ((fun-name (intern op :ga ))
	(string-entry (format nil "~a~s" "c_dd_" op))
	(lisp-name-for-entry (intern (format nil "~a~s" "dd_" op))))
		      
    `(progn
       (ff:def-foreign-call
	   (,lisp-name-for-entry
	    ,string-entry);; e.g. c_qd_sin
	((op1 cqd)(target cqd))
	:returning :void 		         	     
	:arg-checking nil :call-direct t)
	   
       (defmethod ,fun-name ((arg aqd))
	 (let* ((h (alloc-qd)) (in (aqd-q h)))
	   (,lisp-name-for-entry (aqd-q arg) in) h))
       (compile ',fun-name))))

(r sin) (r cos)  (r tan)
(r atan)(r asin) (r acos)
(r sinh)(r cosh) (r atanh)
(r exp)(r log)(r abs)