;; -*- mode: common-lisp; package:ieee; -*-
;; Author: Richard Fateman  21 Sept, 1992
;; CHANGED 7/17/97 for Solaris, Allegro CL 4.3,
;; 
;; So far as I can tell this is SUN Solaris SPARC specific,
;; though other Sun platforms or other UNIX mfgs may now be 
;; more similar than in the past.
;; (convergent evolution?)
;;
;; where to get info? /usr/include/sys/ieeefp.h
;; or                 /usr/include/sys/fpu/ieee.h
;; or                 /usr/include/signal.h, siginfo.h

;; but really, most of what we need is access via C library functions
;; as documented in fpgetround(3C) manual page.

;; This package makes the status and mode flags of the IEEE 754
;; standard for binary floating-point arithmetic available.
;; re-written for Solaris, starting 7/17/97 RJF


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

(eval-when (compile eval load)
  (in-package :ieee))

(defvar *exceptions* 0 "ieee 754 accumulated exceptions")
(defvar *direction*  0 "ieee 754 rounding direction")
(defvar *precision*  0 "ieee 754 rounding precision")
#|
     The  following  floating-point  exception  masks  are  OR-ed
     together to form mask.
             FP_X_INV      /* invalid operation exception */
             FP_X_OFL      /* overflow exception */
             FP_X_UFL      /* underflow exception */
             FP_X_DZ       /* divide-by-zero exception */
             FP_X_IMP      /* imprecise (loss of precision) */
|#

(eval-when (compile eval load)
  
  ;;  enum fp_trap_enable_type 
  (defconstant inexact 0)  ;; in c, fp_trap_inexact. etc
  (defconstant division 1)
  (defconstant underflow 2)
  (defconstant overflow 3)
  (defconstant invalid 4)
  (defconstant domain 5)		;is this used in Solaris ?
  
  ;;enum fp_direction_type
  (defconstant nearest 0)
  (defconstant tozero 1)
  (defconstant negative 2)
  (defconstant positive 3)
  
  ;;enum fp_precision_type   
  (defconstant single 0)
  (defconstant precision_3 1)		;is this used?
  (defconstant double 2)
  (defconstant extended 3)
	   )

;; First see if ieee entry points are loaded, and load them if necessary.
;; these are documented in the C Library Functions part of the manual
;; under fpgetround(3C).

#+ignore  ;;;
(unless (ff:get-entry-point "_fpgetmask")
  (load "" :unreferenced-lib-names 

	'(
  	  "_fpsetmask"
	  "_fpgetmask"
  	  "_fpgetround"
    	  "_fpsetround"
      	  "_fpsetsticky"
       	  "_fpgetsticky"
	  "_fp_rounding_direction"
	  "_swap_rounding_direction"
  	  "_swap_accrued_exceptions"

	  "_log")
	;; ignored? not needed?  :system-libraries '("m")
	))


;; This next form sets up the correspondence between the
;; lisp and C worlds.

(ff:defforeign 'raw_swapRD :entry-point "_swap_rounding_direction"
	       :arguments  '(integer)
	       :arg-checking nil
	       :callback nil
	       :call-direct t
	       :return-type :integer)

(ff:defforeign 'fpgetmask :entry-point "_fpgetmask"
	       :arguments  '()
	       :arg-checking nil
	       :callback nil
	       :call-direct t
	       :return-type :integer)

(ff:defforeign 'fpsetmask :entry-point "_fpsetmask"
	       :arguments  '(integer)
	       :arg-checking nil
	       :callback nil
	       :call-direct t
	       :return-type :integer)


;; when rounding from extended, to what target?
#+ignore (ff:defforeign 'raw_swapRP :entry-point "__swapRP"
	       :arguments  '(integer)
	       :arg-checking nil
	       :callback nil
	       :call-direct t
	       :return-type :integer)

(ff:defforeign 'raw_swapEX :entry-point "_swap_accrued_exceptions"
	       :arguments  '(integer)
	       :arg-checking nil
	       :callback nil
	       :call-direct t
	       :return-type :integer)

;; each one of these routines below sets the global (known to lisp)
;; data, and returns as a value the previous setting of the data.

(defun swapRD (x)(setq *direction* x) (raw_swapRD x))
(defun swapRP (x)(setq *precision* x) (raw_swapRP x))
(defun swapEX (x)(setq *exceptions* x) (raw_swapEX x))

;; getEX() finds out the current exceptions vector and stores it in the
;; lisp variable exceptions, but also restores it to the system vector.

(defun getEX ()(raw_swapEX(setq *exceptions* (raw_swapEX 0))) *exceptions*)

;; (clearEX) changes the exceptions vector to clear all accumulated bits.

(defun clearEX ()(swapEX 0))

#|
 (ieee:tryfp x)  evaluates a single s-expression  x and returns
 two values: The result of evaluating x, and a list
 of the names of flags that were set.

 implementation note:
 we use ignore-errors below, but all we really
 want to do is mask off the software-created lisp
 signals (i.e. "division-by-zero") which are NOT produced
 by the floating point arithmetic, but by programs that
 check for divisors being zero, and put us in a debug loop.
 We do better using the "raw" division instructions as
 made possible by ieee::single-/ and ieee::double-/ below.

;; This version is not for multiprocessing ....
|#

(defmacro stryfp (x) ;simple tryfp  (not multiprocessing)
  "2 values: x evaluated, and flags that were set"
  `(let(oldexceptions)
     (setq oldexceptions(swapEX 0)); save old exceptions
     (values 
      (ignore-errors ,x)		;evaluate x
      (analyzefp(swapEX oldexceptions)))))


;; ieee::analyzefp converts the returned exception flag (an integer)
;; into a list of exception names.  The mapping we use below
;; is particular to the SPARC architecture. The 80x86 (for example)
;; is different.

  
(defun analyzefp(n)
  "convert exception integer to a list of atoms"
  (declare (fixnum n)) ; n<=31 in usage
  (let (ilist 
	(tlist '(inexact division underflow overflow invalid domain)))
    ;; for SPARC these are the flag meanings (yours may be different)
    ;;	fp_inexact	= 0,
    ;;	fp_division	= 1,
    ;;	fp_underflow	= 2,
    ;;	fp_overflow	= 3,
    ;;	fp_invalid	= 4,
    ;;  fp_domain       = 5
    (mapc #'(lambda (h)
	      (when (oddp n)
		(setf ilist (cons h ilist)))
	      (setf n (ash n -1)))
	      tlist)
	ilist))

;;The next two programs must be compiled to have the right side-effects
;;in case of overflow or underflow.
;;They make a single (or double) division happen without first checking 
;;for division by 0.0.  Note that the "normal" division in Allegro does 
;;not set the zero-divide flag because it checks in the software before
;;dividing! 

(defun single-/ (x y)
  (declare (optimize (speed 3)(safety 0)) (inline /))
  (the single-float (/ (the single-float x)(the single-float y))))

(defun double-/ (x y) 
  (declare (optimize (speed 3)(safety 0))
	   (inline /))
  (/ (the double-float x)(the double-float y)))




;; Just for the fun of it, here's an entry to Sun-OS's c-library
;; routine for double-precision real positive log.
;; compare (clog 0.0d0) to (log 0.0d0)
;; compare (clog -1.0d0) to (log -1.0d0)

(ff:defforeign 'clog :entry-point "_log" :arguments '(double-float)
  :return-type :double-float)


#|
.......... A more lisp-ish interface to the IEEE float flags......
 Inquire what the current rounding direction is:
 by simply typing (ieee::direction).

for SPARC,  the result 0=nearest, 1=tozero, 2=positive, 3=negative

 You can change the rounding direction by typing, e.g.
 (setf (ieee::direction) 2)   ;etc for positive


 similarly for (ieee::precision)  the current rounding precision.

 You can determine the condition of the underflow flag by
 (ieee::underflow)   ; t or nil

 you can change the underflow bit(simulating an underflow) 
 (setf (ieee::underflow) nil)  ;clears it
 (setf (ieee::underflow) t)    ;simulates an underflow

 you can clear all the exception flags by
 (ieee::clearEX)

 you can look at all the exception flags by
 (ieee::getEX) 
 (ieee::analyzfp (ieee::getEX))  tells you which are set "in English")

 the global values (in this package) of *direction*, *precision*,
 and *exceptions*
 are not reliable indicators of the ieee modes or status unless they
 have been updated by the access functions defined here..

 you can save and restore the exception bits in one swoop rather
 than individually by using (setq oldbits (swapEX newbits))

|#






(defun direction() 
;;for SPARC,  the result 0=nearest, 1=tozero, 2=positive, 3=negative
  (swapRD *direction*) ;returns old direction
  *direction*)



(defsetf direction ()(val)
;;possible values are constants: nearest, tozero, negative, positive
  `(swapRD ,val))
		
;; get/set rounding precision

(defun precision()
  (swapRP *precision*) ;returns old precision
  *precision*)

(defsetf precision ()(val) 
    ;;possible values are constants: extended, single, double, precision_3n
  `(swapRP ,val))

;; example: is there an underflow? (underflow) tells you

(defun underflow()
  (logbitp #.underflow (getEX)))

;; set underflow on or off (t or nil) by (setf (underflow) nil) etc.

(defsetf underflow ()(val)
  `(let ((%z ,val))
     (if %z (swapEX(logior #.(expt 2 underflow) *exceptions*))
       (swapEX (logandc1 #.(expt 2 underflow) *exceptions*)))
     %z))

(defun overflow()
  (logbitp #.overflow (getEX)))

(defsetf overflow ()(val)
  `(let ((%z ,val))
     (if %z (swapEX(logior #.(expt 2 overflow) *exceptions*))
       (swapEX (logandc1 #.(expt 2 overflow) *exceptions*)))
     %z))

(defun inexact()
  (logbitp #.inexact (getEX)))

(defsetf inexact ()(val)
  `(let ((%z ,val))
     (if %z (swapEX(logior #.(expt 2 inexact) *exceptions*))
       (swapEX (logandc1 #.(expt 2 inexact) *exceptions*)))
     %z))

(defun division()
  (logbitp #.division (getEX)))

(defsetf division ()(val)
  `(let ((%z ,val))
     (if %z (swapEX(logior #.(expt 2 division) *exceptions*))
       (swapEX (logandc1 #.(expt 2 division) *exceptions*)))
     %z))

(defun invalid()
  (logbitp #.invalid (getEX)))

(defsetf invalid ()(val)
  `(let ((%z ,val))
     (if %z (swapEX(logior #.(expt 2 invalid) *exceptions*))
       (swapEX (logandc1 #.(expt 2 invalid) *exceptions*)))
     %z))


#|
-----------example---------


(let (a b old-dir)
  (setq old-dir (swapRD 0))
  (setf (direction) 2)			;round negative
  (setf a (/ 1.0 3.0))
  (setf (direction) 3)			;round positive
  (setf b (/ 1.0 3.0))
  (swapRD old-dir)
  (list a b (- a b)))

;;result is
(0.3333333 0.33333334 2.9802322e-8)

;; for your edification:
 (* 0.25 single-float-epsilon) -->
 2.9802322e-8


;;; some additional data for testing
(defvar big most-positive-single-float)
(defvar inf (* 2.0 big))
(defvar small (/ 1.0 big))
(defvar nan (- inf inf))

(defun test3 (x y) ;; 
  (declare (optimize (speed 3)(safety 0))
	   (double-float x y)
	   (inline / * +))
  (/  (+ (* (the double-float x)(the double-float x))
	 (* (the double-float y)(the double-float y)))
      (+ (the double-float x)(the double-float y))))
|#

#|
 An aside:
 All we use here is __swapRD, RP, EX, but we also
 picked up some additional entrypoints.  Useful ones might include
 nextafter, copysign.  Some of the other math functions
 that you might be tempted to use, like log, use the UNIX
 "matherr" routine, and sometimes gives different results:
 For example, lisp uses complex numbers sometimes, and most of its
 math functions are generic over integers, rationals, single, double.
 Allegro CL's (log 0.0) gives -3.4028232e+38 which I think is wrong;
 it should be #.EXCL::*NEGATIVE-INFINITY-SINGLE*.  But if you
 use the log from the c library (see def'n of clog below)
 you get 3 things (you may not want them all.)

 a message: "log: SING error  " is printed on error-io.
 an answer: #.EXCL::*NEGATIVE-INFINITY-DOUBLE* is returned
 flags set: (DIVISION)

 also, clog is for double floats only.

 end of aside.

 example of use .. tryfp returns a list of two results: the answer
 and then the list of bits that were turned on.

 (in-package :ieee)
 (compile 'ieee::single-/) ;if it has not already been done
 (ieee::tryfp (list (ieee::single-/ 0.0 0.0) (expt 10.0 10000)))
   RESULTS--->
 (#.EXCL::*NAN-SINGLE* #.EXCL::*INFINITY-SINGLE*)
 (IEEE::INVALID IEEE::OVERFLOW IEEE::INEXACT)

 some other comments about lisp and floating point conditions...
 let g(x) = 1+ 3/(1/(x-3))  { a strange way of writing 3*x-8}

 this one gives a non-continuable error for (g 3.0)
 (defun g(x)(+ 1(/ 3.0d0 (/ 1.0d0 (- x 3.0d0)))) )
 this one gives the answer 1.0d0 (correct)
 it also leaves the division-by-zero flag on.
 (defun g(x)(+ 1(double-/ 3.0d0 (double-/ 1.0d0 (- x 3.0d0)))) )

|#

;;; multiprocessing can mess us up;;;;;;;

#|
Allegro CL 4.1 is a multiprocessing lisp system.  It is possible
to have multiple lightweight processes running, in which case it
is possible that there will be several processes computing floating
point results.  Under these circumstances it is necessary for each of
the processes to maintain its own copy of the hardware registers
that pertain to the floating-point status.  The system already keeps
track of the floating-point data registers on a per-process basis,
but apparently not the extra data needed here.

One possible way of dealing with this is to have, for a floating-point
process (or rather, each fp process), the following structure:
|#

;; for example

(use-package :mp)

(defun fp1-hook-fun () ;;should compile this for sure..
		 (swapRD *direction*)
		 (swapRP *precision*)
		 (swapEX *exceptions*))

;(compile 'fp1-hook-fun)

(defmacro tryfp (x)
  "returns 2 values: x evaluated, and flags that were set"
  `(let ((ended (list nil nil nil)) ; use for flag
	 (oldex (swapEX 0))		; old exeptions
	 (theprocess 
	  (make-process :name "fp1-proc"
			:resume-hook #'fp1-hook-fun
			:suspend-hook #'fp1-hook-fun)
	 ))
     (process-enable theprocess)
     (process-preset theprocess
      #'(lambda(ended *package*)
	  (let ((*direction* 0)(*precision* 0)(*exceptions* 0))
	    (funcall #'fp1-hook-fun)
	    (setf (cadr ended) (ignore-errors ,x))
	    (setf (caddr ended) (analyzefp(getEX)))
	    (funcall #'fp1-hook-fun)
	    (setf (car ended) t);; indicate we're done via changing shared val
	    ))
      ended
      *package*)
     (process-wait "waiting for fp1" #'(lambda() (car ended)))
     (swapEX oldex)
     (values (cadr ended)(caddr ended))
	 
     ))

#|
 This section contains non-working code for enabling, on unix, some
of the other stuff...  this is not multi-processing compatible in the
sense that trapping could become broken if there is an attempt at a
(lisp) process switch during the handling of a trap.

(ff:defforeign 'ieee_flags  ;; not used now..
    :arguments
  '(string string string (vector (unsigned-byte 32) 1))
  :return-type :integer) ;; wrong return type?

(ff:defforeign 'ieee_handler 
    :arguments
  '(string string integer) ;not really integer, but fake it?
  :return-type :integer)

(ff:defforeign 'csignal :entry-point "_signal"
    :arguments
  '(integer integer integer) ;;also faked
  :return-type :integer)


(defun foo (signal &optional ignore)
  (format t "~&; got signal ~d~%" signal) t)

 These 2 lines "should" set the handling of all arithmetic exceptions
 to the same kind of handling as a signal-2  (keyboard interrupt)
 in lisp

(ieee_handler "set" "all" (csignal 2 *out* 0))
(excl::add-signal-handler 8 'foo)	;sigfe

|#