r9709: update policy

[kmrcl.git] / math.lisp

diff --git a/math.lisp b/math.lisp
index be89825d27890bc0e6268a751fe363832662ad18..fcbab3c4baa2e922b0eba049af0bd04079a19739 100644 (file)
--- a/math.lisp
+++ b/math.lisp
@@ -34,3 +34,65 @@
   `(if (stringp ,obj)
       (parse-integer ,obj)
      ,obj))
+
+(defun histogram (v n-bins &key min max)
+  (declare (fixnum n-bins))
+  (when (listp v)
+    (setq v (coerce v 'vector)))
+  (when (zerop (length v))
+    (return-from histogram (values nil nil nil)) )
+  (let ((n (length v))
+       (bins (make-array n-bins :element-type 'integer :initial-element 0))
+       found-min found-max)
+    (declare (fixnum n))
+    (unless (and min max)
+      (setq found-min (aref v 0)
+           found-max (aref v 0))
+      (loop for i fixnum from 1 to (1- n)
+         do
+           (let ((x (aref v i)))
+             (cond
+              ((> x found-max)
+               (setq found-max x))
+              ((< x found-min)
+               (setq found-min x)))))
+      (unless min
+       (setq min found-min))
+      (unless max
+       (setq max found-max)))
+    (let ((width (/ (- max min) n-bins)))
+      (setq width (+ width (* double-float-epsilon width)))
+      (dotimes (i n)
+       (let ((bin (nth-value 0 (truncate (- (aref v i) min) width))))
+         (declare (fixnum bin))
+         (when (and (not (minusp bin))
+                    (< bin n-bins))
+           (incf (aref bins bin))))))
+    (values bins min max)))
+             
+
+(defun fixnum-width ()
+  (nth-value 0 (truncate (+ (/ (log (1+ most-positive-fixnum)) (log 2)) .5))))
+
+(defun scaled-epsilon (float &optional (operation '+))
+  "Return the smallest number that would return a value different from
+  FLOAT if OPERATION were applied to FLOAT and this number.  OPERATION
+  should be either + or -, and defauls to +."
+  (multiple-value-bind (significand exponent)
+      (decode-float float)
+    (multiple-value-bind (1.0-significand 1.0-exponent)
+       (decode-float (float 1.0 float))
+      (if (and (eq operation '-)
+              (= significand 1.0-significand))
+         (scale-float (typecase float
+                        (short-float short-float-negative-epsilon)
+                        (single-float single-float-negative-epsilon)
+                        (double-float double-float-negative-epsilon)
+                        (long-float long-float-negative-epsilon))
+                      (- exponent 1.0-exponent))
+       (scale-float (typecase float
+                      (short-float short-float-epsilon)
+                      (single-float single-float-epsilon)
+                      (double-float double-float-epsilon)
+                      (long-float long-float-epsilon))
+                    (- exponent 1.0-exponent))))))