math.lisp

   1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*-
   2 ;;;; *************************************************************************
   3 ;;;; FILE IDENTIFICATION
   4 ;;;;
   5 ;;;; Name:          math.lisp
   6 ;;;; Purpose:       General purpose math functions
   7 ;;;; Programmer:    Kevin M. Rosenberg
   8 ;;;; Date Started:  Nov 2002
   9 ;;;;
  10 ;;;; This file, part of KMRCL, is Copyright (c) 2002 by Kevin M. Rosenberg
  11 ;;;;
  12 ;;;; KMRCL users are granted the rights to distribute and use this software
  13 ;;;; as governed by the terms of the Lisp Lesser GNU Public License
  14 ;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL.
  15 ;;;; *************************************************************************
  16
  17
  18 (in-package #:kmrcl)
  19
  20 (defun deriv (f dx)
  21   #'(lambda (x)
  22       (/ (- (funcall f (+ x dx)) (funcall f x))
  23          dx)))
  24
  25 (defun sin^ (x)
  26     (funcall (deriv #'sin 1d-8) x))
  27
  28 ;;; (sin^ pi)
  29
  30 (defmacro ensure-integer (obj)
  31   "Ensure object is an integer. If it is a string, then parse it"
  32   `(if (stringp ,obj)
  33       (parse-integer ,obj)
  34      ,obj))
  35
  36 (defun histogram (v n-bins &key min max)
  37   (declare (fixnum n-bins))
  38   (when (listp v)
  39     (setq v (coerce v 'vector)))
  40   (when (zerop (length v))
  41     (return-from histogram (values nil nil nil)) )
  42   (let ((n (length v))
  43         (bins (make-array n-bins :element-type 'integer :initial-element 0))
  44         found-min found-max)
  45     (declare (fixnum n))
  46     (unless (and min max)
  47       (setq found-min (aref v 0)
  48             found-max (aref v 0))
  49       (loop for i fixnum from 1 to (1- n)
  50           do
  51             (let ((x (aref v i)))
  52               (cond
  53                ((> x found-max)
  54                 (setq found-max x))
  55                ((< x found-min)
  56                 (setq found-min x)))))
  57       (unless min
  58         (setq min found-min))
  59       (unless max
  60         (setq max found-max)))
  61     (let ((width (/ (- max min) n-bins)))
  62       (setq width (+ width (* double-float-epsilon width)))
  63       (dotimes (i n)
  64         (let ((bin (nth-value 0 (truncate (- (aref v i) min) width))))
  65           (declare (fixnum bin))
  66           (when (and (not (minusp bin))
  67                      (< bin n-bins))
  68             (incf (aref bins bin))))))
  69     (values bins min max)))
  70
  71
  72 (defun fixnum-width ()
  73   (nth-value 0 (truncate (+ (/ (log (1+ most-positive-fixnum)) (log 2)) .5))))
  74
  75 (defun scaled-epsilon (float &optional (operation '+))
  76   "Return the smallest number that would return a value different from
  77   FLOAT if OPERATION were applied to FLOAT and this number.  OPERATION
  78   should be either + or -, and defauls to +."
  79   (multiple-value-bind (significand exponent)
  80       (decode-float float)
  81     (multiple-value-bind (1.0-significand 1.0-exponent)
  82         (decode-float (float 1.0 float))
  83       (if (and (eq operation '-)
  84                (= significand 1.0-significand))
  85           (scale-float (typecase float
  86                          (short-float short-float-negative-epsilon)
  87                          (single-float single-float-negative-epsilon)
  88                          (double-float double-float-negative-epsilon)
  89                          (long-float long-float-negative-epsilon))
  90                        (- exponent 1.0-exponent))
  91         (scale-float (typecase float
  92                        (short-float short-float-epsilon)
  93                        (single-float single-float-epsilon)
  94                        (double-float double-float-epsilon)
  95                        (long-float long-float-epsilon))
  96                      (- exponent 1.0-exponent))))))
  97
  98 (defun sinc (x)
  99   (if (zerop x)
 100       1d0
 101     (let ((x (coerce x 'double-float)))
 102       (/ (sin x) x))))
 103
 104
 105 (defun numbers-within-percentage (a b percent)
 106   "Determines if two numbers are equal within a percentage difference."
 107   (let ((abs-diff (* 0.01 percent 0.5 (+ (abs a) (abs b)))))
 108     (< (abs (- a b)) abs-diff)))