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

;; NaNPack: Retrospective floating-point Diagnostics using NotANumbers
;; reference: Steven Lumetta's CS283 term project report
;; for CS283, Fall, 1992,
;; Univ. of Calif. at Berkeley

;; Authors: Steven Lumetta  Richard Fateman, December, 1992

;; some features are specific to Allegro CL 4.1,
;; some features are specific to Sun Microsystems SPARC

;; Generalizations: 
;; We could extend to double-NaN operations.
;; More or less information might be reported for given operations.

(require "ieee")

(eval-when (compile eval load)
  (or (find-package :nanpack) (make-package :nanpack))
  (use-package :ieee))

;; we do not export anything at the moment..

;; This proclaim should
;; speed things up on compiling; disable some safety considerations
;; that may have the effect of disabling our diagnostics.

(proclaim '(optimize speed (debug 0)(safety 0)))

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

(defparameter *NaN-data-size* 100)

;; *string-stream* is used to collect 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*)

(defvar *NaN-data* (make-array *NaN-data-size* :adjustable t
			       :element-type 'string
			       :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*)

;;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 *NaN-skip* property

(defun add-skip-fn (fn)
  (setf (get fn '*NaN-skip*) t))

(defun rem-skip-fn (fn)
  (setf (get fn '*NaN-skip*) nil))

;;; THIS IS THE MAIN FUNCTION

(defun make-NaN ()
  (let ((frame-ref (debug::newest-frame
		    sys:*current-stack-group*
		    :visible-only-p t))
	(return-val (new-float 1.0)))
    (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) '*NaN-skip*)))
	 (let ((*terminal-io* *string-stream*))
	   (tpl::zoom-print-frame frame-ref nil nil)
	   (put-nan-tag return-val 
			(vector-push-extend
			 (get-output-stream-string *string-stream*)
			 *NaN-data* *NaN-data-size*))
	      return-val)))))

(add-skip-fn 'make-NaN)

;; 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))


; Next 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--I will also 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))))

; The function nan-p returns true if a single-float is a NaN.

(defun NaN-p (n)
  (declare (single-float n))
  ;; (eql #x7F800001 (logand #x7F800001 (sys:memref n -6 4 :unsigned-long)))
  (excl::nan-p n) ;use built-in
  )

(defun get-NaN-info (n)
  (let ((tag (if (single-float-p n) (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)))))

(defmacro nanrept(a b name)
  "report generation for a 2-arg nan-result"
  `(if (NaN-p ,a)
	      (if (NaN-p ,b)
		  (format 
		   *string-stream*
		   "~s called with NaN#~4,'0d and NaN#~4,'0d from:~%"
		   ,name (get-NaN-tag ,a) (get-NaN-tag ,b))
		(format *string-stream*
			"~s called with NaN#~4,'0d and ~a from:~%"
			,name (get-NaN-tag ,a) ,b))
	    (if (NaN-p ,b)
		(format *string-stream*
			"~s called with ~a and NaN#~4,'0d from:~%"
			,name ,a (get-NaN-tag ,b))
	      (format *string-stream* 
		      "~s called with ~a and ~a from:~%" ,name ,a ,b))))

;; 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)
	(progn  (nanrept a b "Divide")  (make-NaN))
      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_EX (ieee::swapEX 0))
	 (value (the single-float 
		  (/ (the single-float a) (the single-float b))))
	 (new_EX (ieee::analyzefp (ieee::swapEX save_EX))))
    (if (NaN-p value)
	(progn
	  (if new_EX
	      (format *string-stream* "Exceptions ~{~a~^ ~} in~%" new_EX))
	  (nanrept a b "Divide")
	  (make-NaN))
      value)))




(add-skip-fn 'divide)

(defun mult (a b)
  "nanpack version of single-float multiply"
  (declare  (inline *)	   (single-float a b))
  (let* ((save_EX (ieee::swapEX 0))
	 (value (the single-float 
		  (* (the single-float a) (the single-float b))))
	 (new_EX (ieee::analyzefp (ieee::swapEX save_EX))))
    (if (NaN-p value)
	(progn
	  (if new_EX
	      (format *string-stream* "Exceptions ~{~a~^ ~} in~%" new_EX))
	  (nanrept a b "Mult")
	  (make-NaN))
      value)))

(add-skip-fn 'mult)

(defun add (a b)
  (declare (inline +)	   (single-float a b))
  (let* ((save_EX (ieee::swapEX 0))
	 (value (the single-float 
		  (+ (the single-float a) (the single-float b))))
	 (new_EX (ieee::analyzefp (ieee::swapEX save_EX))))
    (if (NaN-p value)
	(progn
	  (if new_EX
	      (format *string-stream* "Exceptions ~{~a~^ ~} in~%" new_EX))
	  (narept a b "Add")
	  (make-NaN))
      value)))

(add-skip-fn 'add)

(defun sub (a b)
  (declare (optimize (speed 3)(safety 0)) (inline -)
	   (single-float a b))
  (let* ((save_EX (ieee::swapEX 0))
	 (value (the single-float 
		  (- (the single-float a) (the single-float b))))
	 (new_EX (ieee::analyzefp (ieee::swapEX save_EX))))
    (if (NaN-p value)
	(progn
	  (if new_EX
	      (format *string-stream* "Exceptions ~{~a~^ ~} in~%" new_EX))
	  (nanrept a b "Sub")
	  (make-NaN))
      value)))

(add-skip-fn 'sub)


(defun root (a)
  (declare (optimize (speed 3)(safety 0)) (inline sqrt)
	   (single-float a))
  (if (< a 0.0)
      (progn
	(format *string-stream* 
		"Exceptions DOMAIN in~%Root called with ~a from:~%" a)
	(make-NaN))
    (let* ((save_EX (ieee::swapEX 0))
	   (value (the single-float (sqrt (the single-float a))))
	   (new_EX (ieee::analyzefp (ieee::swapEX save_EX))))
      (if (NaN-p value)
	  (progn
	    (if new_EX
		(format *string-stream* "Exceptions ~{~a~^ ~} in~%" new_EX))
	    (if (NaN-p a)
		(format *string-stream* "Root called with NaN#~4,'0d from:~%"
			(get-NaN-tag a))
	      (format *string-stream* "Root called with ~a from:~%" a))
	    (make-NaN))
	value))))

(add-skip-fn 'root)


;;;;;;;;;;; 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)
      (format *string-stream* "~s[~a] from:" name i)
      (setf (aref a i) (make-nan)))
    a))

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

(defun find-root (func x0 x1 tolerance &optional 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))
    (let ((x2 (sub x1 (mult val
			    (divide (sub x1 x0)
				    (sub val oldval))))))
      (format t "X0: ~a, X1: ~a, f(X0): ~a, f(X1): ~a  ->  X2: ~a~%"
	      x0 x1 oldval val x2)
      (setf x0 x1)
      (setf x1 x2))))


(defun test2 (x)
  (add (root (mult (sub x 1.0) (sub x 4.0)))
       (mult (sub x 7.0) (sub x 7.0))))


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

;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)

;"#0003: Exceptions INVALID in
;Add called with NaN#0002 and 10.629862 from:
;   (TEST2 3.739653)"

;USER(40): (get-nan-info 2)

;
;"#0002: Exceptions DOMAIN in
;Root called with -0.7132602 from:
;   (TEST2 3.739653)"