;; -*- mode: common-lisp; package: nanpack; -*-

#|
 NaNPack: Retrospective Floating-Point Diagnostics using Not-A-Numbers.
 reference: Steven Lumetta's CS283 term project report
 for CS283, Fall, 1992.
 Univ. of Calif. at Berkeley

 Authors: Steven Lumetta,  Richard Fateman, Dec. 92, Jan. 93

 some features needed are specific to Allegro CL 4.1,
 some features needed may be specific to Sun Microsystems SPARC

 Future generalizations: 
 We could duplicate this facility for double-NaN operations.
 See how to package up other information in NaNs 
 following the template in make-initialized-array.
|#



(eval-when (compile eval load)
  (require "ieee")
  (or (find-package :nanpack) (def-package :nanpack))
  (use-package :ieee)
  (in-package :nanpack)
  (export '(add subtract divide real-sqrt get-nan-info)))


;; This declaim should speed things up when run compiled, as well as
;; disable some safety considerations that may have the effect 
;; of disabling our diagnostics.

(declaim (optimize speed (debug 0)(safety 0)))

;; initial size of NaN-date array.  This array will be extended as needed.

(defparameter *NaN-data-size* 100)


(defvar *NaN-data* (make-array *NaN-data-size* :adjustable t
			       :fill-pointer 0))

(setf (fill-pointer *NaN-data*) 0)	; reloading this file resets array
(vector-push-extend "Data entered as NaN" *NaN-data* *NaN-data-size*)

;; *string-stream* is used to collect stack-frame information from the
;; debugger so that it can be inserted as the "value" of a NaN.

(defvar *string-stream* (make-string-output-stream))

(get-output-stream-string *string-stream*);;clear string-stream


#|
Functions that we don't want to see reported on the stack as
user functions with respect to NaN errors are indicated by setting
their *dont-report-on-stack* property.
|#

(defun dont-report-on-stack (fn) (setf (get fn '*dont-report-on-stack*) t))

#|
Put-NaN-tag is a macro to set the 22-bit tag of a NaN float--this must be
a macro to insure that the number itself gets modified and not a copy
of the number-- we assume that the exponent bits indicate NaN
and that the least significant mantissa bit is set.
|#

(defmacro put-NaN-tag (n tag)
  (declare  (single-float n))
  `(sys::.inv-memref ,n -6 4 :unsigned-long 
		     (logior #x7F800001 (ash ,tag 1))))


;;; THIS IS THE MAIN FUNCTION

(defun make-NaN (datum)
  "Make a single-nan with datum for data value"
  (let ((return-val (new-float 1.0)))
    (put-nan-tag return-val
		 (vector-push-extend
		  datum			;usually but not nec. a string
		 *NaN-data* *NaN-data-size*))
    return-val))

(dont-report-on-stack 'make-NaN)

(defun where-am-i()
  "return a string explaining what's on the stack now {Allegro CL specific}"
  (let ((frame-ref (debug::newest-frame
		    sys:*current-stack-group*
		    :visible-only-p t)))
    (do ((frame-ref 
	  (debug::next-older-frame frame-ref)
	  (debug::next-older-frame frame-ref)))
	((and (debug::frame-visible-p frame-ref)
	      (eq (debug::frame-type frame-ref) :function)
	      (not (get (debug::frame-name frame-ref)
			'*dont-report-on-stack*)))
	 (let ((*terminal-io* *string-stream*))
	   (tpl::zoom-print-frame frame-ref nil nil)
	   (get-output-stream-string *string-stream*)
	   )))))

(dont-report-on-stack 'where-am-i)

;; assure that we create a new float at run-time
(defun new-float (x) (+ x 1.0))       

(defun get-NaN-tag (n)
  "return the 22-bit tag of a single- float NaN."
  (declare  (single-float n))
  (ash (logand (sys:memref n -6 4 :unsigned-long) #x003FFFFF) -1))

;; The function nan-p returns non-nil if n is a NaN.
;; another useful definition would be--
;; (eql #x7F800001 (logand #x7F800001 (sys:memref n -6 4 :unsigned-long)))

(defun NaN-p (n)
  (declare (single-float n))
  (excl::nan-p n) )


(defun get-NaN-info (n)
  "given a NaN, returns stored info in it;
   given an integer N, gives info on the Nth NaN"
  (let ((tag (if (typep  n 'single-float) (get-NaN-tag n) n)))
    (if (>= tag (fill-pointer *NaN-data*))
	"This NaN appears to originate outside NaNPack."
      (format nil "#~4,'0d: ~a" tag (aref *NaN-data* tag)))))

;; two utility functions for use by arithmetic programs

(defun nanrept2(a b name)
  "report generation for a 2-arg nan-result"
  (concatenate 'string
     (if (NaN-p a)
	      (if (NaN-p b)
		  (format  nil
		   "~a called with NaN#~4,'0d and NaN#~4,'0d"
		   name (get-NaN-tag a) (get-NaN-tag b))
		
		(format nil
			"~a called with NaN#~4,'0d and ~a"
			name (get-NaN-tag a) b))
       
       (if (NaN-p b)
	   (format nil
		   "~a called with ~a and NaN#~4,'0d"
		   name a (get-NaN-tag b))
	 (format nil
		 "~a called with ~a and ~a" name a b)))
     (format nil " from ~a" (where-am-i))))

(dont-report-on-stack 'nanrept2)

(defun nanrept1(a name)
  "report generation for a 1-arg nan-result"
  (concatenate 'string
     (if (NaN-p a)
	 (format  nil  "~a called with NaN#~4,'0d"  name (get-NaN-tag a))
       (format nil "~a called with ~a" name a))
     (format nil " from ~a" (where-am-i))))

(dont-report-on-stack 'nanrept1)

;; these are examples of functions implementing basic arithmetic
;; with nan-trapping and exception reporting.

#+ignore
(defun divide (a b)
  ;;minimal version...
  "nanpack version of single-float division"
  (declare  (single-float a b))
  (let ((value  (/ a b)))
    (if (NaN-p value)
	(make-NaNf  (nanrept2 a b "Divide"))
      value)))

;; a more elaborate version, reporting exception flags as well.

(defun divide (a b)
  "nanpack version of single-float division"
  (declare  (inline /)  (single-float a b))
  
  (let* (
	 ;;save old EX flags set them to 0
	 (save_EX (ieee::swapEX 0))	
	 
	 ;; evaluate the division ( must compile this...)
	 (value (the single-float 
		  (/ (the single-float a) (the single-float b))))
	 
	 ;; analyze new settings of EX flags, restore old ones
	 (new_EX (ieee::analyzefp (ieee::swapEX save_EX))))

    (if (NaN-p value)
	(make-NaN 
	 (concatenate 'string
	   
	   (if new_EX
	       (format nil "Exception~p ~{~a~^ ~} in " (length new_EX) new_EX)
	     "")
	   
	   (nanrept2 a b "Divide")))
      
      ;; normal return value
      value)))

(dont-report-on-stack 'divide)


(defun mult (a b)
  "nanpack version of single-float multiply"
  (declare  (inline *)  (single-float a b))
  
  (let ((value (the single-float 
		 (* (the single-float a) (the single-float b)))))
    
    (if (NaN-p value)
	(make-NaN (nanrept2 a b "Mult"))
      value)))

(dont-report-on-stack 'mult)

(defun add (a b)
  (declare (inline +)	(single-float a b))
    (let ((value (the single-float 
		 (+ (the single-float a) (the single-float b)))))
    
    (if (NaN-p value)
	(make-NaN (nanrept2 a b "Add"))
      value)))

(dont-report-on-stack 'add)

(defun subtract (a b)
  (declare (optimize (speed 3)(safety 0)) (inline -)
	   (single-float a b))
      (let ((value (the single-float 
		 (- (the single-float a) (the single-float b)))))
    
    (if (NaN-p value)
	(make-NaN (nanrept2 a b "Subtract"))
      value)))

(dont-report-on-stack 'subtract)


(defun real-sqrt (a)
  "nanpack single-float square root"
  (declare (inline sqrt)  (single-float a))

  (if (or (eq 'excl::*nan-single* (NaN-p a))(< a 0.0))
      (make-NaN (nanrept1 a "Real-sqrt"))
    (the single-float (sqrt (the single-float a)))))

(dont-report-on-stack 'real-sqrt)


;;;;;;;;;;; some examples of use ;;;;;;;;;


(defun make-initialized-array (n name)
  "stuff into an array of length n, single-float NaNs with names Name[i]"
  (let ((a (make-array n :element-type 'single-float)))
    (dotimes (i n)			;i:= 0 to n-1
      (setf (aref a i)
	(make-nan (format nil "uninitialized ~a[~a]" name i))))
    a))

(defun make-initialized-2D-array (n m name)
  "stuff into an array of size n m, single-float NaNs with names Name[i,j]"
  (let ((a (make-array (list n m) :element-type 'single-float)))
    (dotimes (i n)
      (dotimes (j m)
	(setf (aref a i j)
	  (make-nan 
	   (format nil "uninitialized ~a[~a,~a]" name i j) ))))
    a))

;; a secant method root-finder that demonstrates the 
;; usefulness of retrospective diagnosis of Nans..

(defvar find-root-debug t)

(defun find-root (func x0 x1 tolerance &optional max)
  "return iterate x[n] such that |f(x[n])|< tolerance or n>max"
  (do ((oldval (funcall func x0) val)
       (val (funcall func x1) (funcall func x1))
       (count 0 (1+ count)))
      ((or (and (numberp max) (>= count max))
	   (nan-p val)
	   (< (abs val) tolerance))
       (values x1 val)) ;x1 is the final iterate, val is f(x1).
    (let ((x2 (subtract x1 (mult val
			    (divide (subtract x1 x0)
				    (subtract val oldval))))))
      ;; perhaps print info on the iteration
      (if find-root-debug
	  (format t "X0: ~a, X1: ~a, f(X0): ~a, f(X1): ~a  ->  X2: ~a~%"
		  x0 x1 oldval val x2))
      
      (setf x0 x1 x1 x2))))


(defun test2 (x)
  (add (real-sqrt (mult (subtract x 1.0) (subtract x 4.0)))
       (mult (subtract x 7.0) (subtract x 7.0))))

#| Example

USER(38): (setf a (multiple-value-bind (trash answer) 
	     (find-root 'test2 6.0 6.5 0.0001)
	     answer))

;; somewhat reformatted printout...

X0:  6.000, X1:  6.500, f(X0):  4.162, f(X1):  3.958  ->  X2: 16.193
X0:  6.500, X1: 16.193, f(X0):  3.958, f(X1): 98.117  ->  X2:  6.093
X0: 16.193, X1:  6.093, f(X0): 98.117, f(X1):  4.088  ->  X2:  5.653
X0:  6.093, X1:  5.653, f(X0):  4.088, f(X1):  4.587  ->  X2:  9.689
X0:  5.653, X1:  9.689, f(X0):  4.587, f(X1): 14.258  ->  X2:  3.740
#.EXCL::*NAN-SINGLE*

USER(39): (get-nan-info a)

"#0006: Add called with NaN#0005 and 10.629862 from    (TEST2 3.739653)"

USER(40): (get-nan-info 5)

"#0005: Real-sqrt called with -0.7132602 from    (TEST2 3.739653)"

;;note, the argument 3.73... will not be visible if test2 is compiled.
|#