;;-*- mode: common-lisp; -*-

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

(defun two-sum(a b)			;jrs's short paper theorem2
  (declare (double-float a b)) 
  (let*((x (+ a b))
	(bv (- x a))
	(ar (- a (- x bv)))
	(y (+ ar (- b bv))))
    (declare (double-float a b x y av ar bv br))     
    (values x y)))

(defun fast-two-sum(a b)  ;;assumes |a| > |b|, both double-floats.
  (declare (double-float a b) (optimize (speed 3)(safety 0)))
  (let ((x (+ a b)))
    (declare (double-float x))
    (values x (-  b  (- x a)))))



(defun numerge (e f)
  ;;non-destructive merge. (similar CL merge is destructive)
  (cond ((null e) f)((null f) e)
	((> (the double-float (car e)) 
	    (the double-float (car f)))
	 (cons (car e)(numerge (cdr e) f)))
	(t
	 (cons (car f)(numerge e (cdr f))))))
	
(defun linear-expansion-sum (e f)	;from jrs' description
  (let ((g (numerge e f))
	(bqi 0.0d0);;big Q[i], i=2 initially
	(lqi 0.0d0);;little q[i]
	(r  0.0d0)
	(htemp 0.0d0)
	(h nil))
    (declare (double-float bqi lqi r htemp) (optimize (speed 3)(safety 0)))
    (multiple-value-setq 
	(bqi lqi)
      (fast-two-sum (cadr g)(car g)))	;bq2,lq2 =twosum of g2, g1
    (setf g (cddr g))
    (do ((i g (cdr i)))
	((null i)			;if at end
	 (unless (= lqi 0.0d0) (push lqi h));tack last two pieces on
	 (unless (= bqi 0.0d0) (push bqi h))
	 h);; return the preferred ordering
      ;; this is the inside of the loop	  
      (multiple-value-setq   (r htemp)
	(fast-two-sum (car i)lqi))
      (multiple-value-setq   (bqi lqi)
	(two-sum bqi r))
      (unless (= htemp 0.0d0) (push htemp h)))))

    
(defun split(a)
  (let* ((c (* a #.(+ 1 (expt 2 26)))) ;compile time constant
	(abig (- c a))
	(ahi (- c abig)))
    (values (- c abig) (- a ahi))))

(defun two-product (a b)
  (multiple-value-bind (ahi alo)
      (split a)
    (multiple-value-bind (bh blo)
	(split b)
      (let* ((x (* a b))
	     (err1 (- x (* ahi bhi)))
	     (err2 (- err1 (* alo bhi)))
	     (err3 (- err2 (* ahi blo))))
	(values x (- (* alo blo) err3))))))
	
(defun scale-expansion(e b)
  (multiple-value-bind
      (q2 h1) (two-product (car e) b)
    (do ((i (cdr e)(cdr i)))
	((null i) (push bq h) h)))	;return value
  ;;; this is wrong.  Need to re-think this from
  ;;; arrays to lists
	(multiple-value-setq (bti ti)
	  (two-product (car i) b))
	(multiple-value-setq (qti hi)
	  (two-sum qti t))
	(multiple-value-setq (qti hi)
	(fast-twosum bti qi)))
    
;; testing

(defun les(x y)(linear-expansion-sum x y))

(defun double(x) (setf x(mapcar #'double-float x))
       (les x x))

(defparameter t1 (list (expt 2 60) (expt 2 30) (expt  2 10) (expt 2 -40)))

;;(double t1) => (2.3058430113611794d+18 1.8189894035458565d-12)
;; correct answer via mathematica is 
2.305843011361179648000000000001818989403545856475830078125000d18

:::SYMBOLIC

(defun l-s(e f)
  (newsym "W" -1)
  (newsymvarsclear) 
  (let* ((sumlen(+ (length e)(length f)))
	 (res (linear-expansion-sum-s sumlen)))
	 
    `(lambda (e f);;
       (setf g (merge e f))
      ,(lda-let  `(let*,(append (list (list 'acount sumlen))
		     (genelts sumlen)
	    
	       (newsymvars))
	    		      
		      
		       ,@(append
		 (apply #'slist1 res)
	      '((values acount bstack))))
	      )
	      
	      )))


(defun genelts(s)
  (let ((i -1))
  (mapcar #'(lambda(r)(incf i)`(,r (elt g ,i)))
			      (pp (vpm 0 "G")  s))))
  
(defun les-s (n)(linear-expansion-sum-s n))

(defun linear-expansion-sum-s (sumlen)
;;symbolic.  sumlen = sum of lengths of e and f
  (let ((g ;;(numerge e f)
	 (pp (vpm 0 "G")  sumlen)
	 )
	
	(bqi nil)
	(lqi nil)
	(r  nil)
	(htemp nil)
	(h nil))
    
    (multiple-value-setq 
	(bqi lqi)
      (fast-two-sum-s (cadr g)(car g)))	;bq2,lq2 =twosum of g2, g1
    (setf g (cddr g))
    (do ((i g (cdr i)))
	((null i)			;if at end
	 (push bqi h)
	 (push lqi h)			;we can't tell if they're zero..

	 (reverse h)
	 );; return the preferred ordering.
      ;;this is the inside of the  loop
      (multiple-value-setq   (r htemp)
	(fast-two-sum-s (car i)lqi))
      (multiple-value-setq   (bqi lqi)
	(two-sum-s bqi r))
      (push htemp h))))


(defun fast-two-sum-s(a b)  ;;assumes |a| > |b|, both double-floats.
  (let ((x (+s a b)))
    (values x (newname(-s  b  (-s x a))))))

(defun -s( &rest k)`(- ,@ k))

(defun two-sum-s(a b)			;jrs's short paper theorem2
  (let*((x (r+s a b))
	(bv (r-s x a))
	(ar (r-s a (-s x bv)))
	(y (r+s ar (-s b bv))))
    (values x y)))

#+ignore
(defun rename (h)
  (cond ((atom h)h)
  (t (let ((n (makename)))
       (push `(,n ,h) vars) n))))

(defun r+s (g h)(newname(+s g h)))
(defun r*s (g h)(newname(*s g h)))
(defun r-s (g h)(newname(-s g h)))

(setf (get '- 'simp) #'(lambda(r)(cons '- r))) ;; don't simplify -..

(defun slist(&rest s)
  ;; generate code to produce a pair acount bstack
  ;; where acount = number of non-zero elements in s
  ;; bstack is the list of them
  (let ((ccode  (apply #'slist1 (reverse s))))
    `(let ((acount 0)
	   (bstack nil))
       ,@ccode
       (values acount bstack))))

(defun slist1(&rest s)
  (cond((null s) '())
       (t
	(cons
	 `(when (= ,(car s) 0.0)(push ,(car s) bstack) (incf acount))
	 (apply #'slist1 (cdr s))))))


;;;
;;;multi-term stuff

(defun +m (&rest z)
  (if (not(= (length z 2)))(error "fix +m to do other than 2 args ~s" z))
  `(,(intern(format nil "MT+~a&~a" (mtlength (car z))
		(mtlength(cadr z)))),@(mapcar #'mtbody z	)))
(defun mtlength(z)
  (cond((and (listp z)(eq (car z) 'mt)) (caddr z))
       (t 1)))
(defun mtbody(z)
  (cond((and (listp z)(eq (car z) 'mt)) (cadr z))
       (t z)))

(defun +m (&rest z)
  (if (not(= (length z 2)))(error "fix +m to do other than 2 args ~s" z))
  (let ((l0 (mtlength (elt z 0)))
	(l1 (mtlength (elt z 1))))
    `(MT (,(intern(format nil "MT+~a&~a" l0 l1)),@(mapcar #'mtbody z))
	 ,(+ l0 l1))))

(defun *m (&rest z)
  (if (not(= (length z 2)))(error "fix *m to do other than 2 args ~s" z))
  (let ((l0 (mtlength (elt z 0)))
	(l1 (mtlength (elt z 1))))
    `(MT (,(intern(format nil "MT*~a&~a" l0 l1)),@(mapcar #'mtbody z))
	 ,(+ l0 l1))))


(setf (get 'mt 'simp) #'(lambda(r)(simp (car r)))) ;; (mt a 4) => a


#|
;; adding a 4 and 2 element  list


(LAMBDA (E F)
  (SETF G (MERGE E F))
  (LET* ((ACOUNT 5)
         (G0 (ELT G -24))
         (G1 (ELT G -25))
         (G2 (ELT G -26))
         (G3 (ELT G -27))
         (G4 (ELT G -28))
         (W0 (- G0 (- (+ G0 G1) G1)))
         (W1 (- W0 (- (+ W0 G2) G2)))
         (W2 (+ G0 W0 G1 G2))
         (W3 (- W2 (+ G0 G1)))
         (W4 (- (+ G0 G1) (- W2 W3)))
         (W5 (+ W0 G2 W4 (* -1 W3)))
         (W6 (- W5 (- (+ G3 W5) G3)))
         (W7 (+ W2 G3 W5))
         (W8 (- W7 W2))
         (W9 (- W2 (- W7 W8)))
         (W10 (+ G3 W5 W9 (* -1 W8)))
         (W11 (- W10 (- (+ W10 G4) G4)))
         (W12 (+ W10 G4 W7))
         (W13 (- W12 W7))
         (W14 (- W7 (- W12 W13)))
         (W15 (+ W10 W14 G4 (* -1 W13))))
    (WHEN (= W15 0.0) (PUSH W15 BSTACK) (INCF ACOUNT))
    (WHEN (= W12 0.0) (PUSH W12 BSTACK) (INCF ACOUNT))
    (WHEN (= W11 0.0) (PUSH W11 BSTACK) (INCF ACOUNT))
    (WHEN (= W6 0.0) (PUSH W6 BSTACK) (INCF ACOUNT))
    (WHEN (= W1 0.0) (PUSH W1 BSTACK) (INCF ACOUNT))
    (VALUES ACOUNT BSTACK)))

|#