lsfloat.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
   3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   4 ;;; Since 1991, ANSI was finally finished.  Edited for ANSI Common Lisp.
   5
   6 ;;;; lsfloat -- Floating point specs and transcendental functions
   7 ;;;;
   8 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   9 ;;;; unrestricted use.
  10 ;;;;
  11 ;;;; Common Lisp allows for four different floating point types that need
  12 ;;;; not be distinct. For statistical work, the type I prefer to use is
  13 ;;;; the one that is closest to a C double. This type is named stat-float.
  14 ;;;; By setting the variable *read-default-float-format* to this type, you
  15 ;;;; insure that data entered as floating point data is read in with this
  16 ;;;; type. The problem arises with data read as integers that is passed to
  17 ;;;; a transcendental, like sqrt. Floating point contagion rules say these
  18 ;;;; integers are to be converted to type single-float. Unless single-float
  19 ;;;; is equivalent to C double, as it is in Mac CL and KCL, this is not
  20 ;;;; what I want. Hence this file redefines the transcendentals to first
  21 ;;;; coerce their arguments to stat-float before applying the built-in
  22 ;;;; functions.
  23 ;;;;
  24 ;;;; No actual modifications to the transcendentals are needed if
  25 ;;;; single-float is the same as stat-float. The fearure
  26 ;;;; :stat-float-is-double-float is used to indicate this.
  27 ;;;;
  28 ;;;; KCL NOTE:
  29 ;;;; In (A)KCL the type short-float corresponds to C float and the types
  30 ;;;; single-float, double-float and long-float correspond to C double.
  31 ;;;; But in the implementation of the transcendentals (A)KCL coerces
  32 ;;;; rationals to short-float, not single-float. CLtL1 is a little vague
  33 ;;;; on this (it talks about "single precision") but CLtL2 clarifies that
  34 ;;;; rationals should produce single-float results. So (A)KCL is wrong, at
  35 ;;;; least relative to the clarification in CLtL2. I therefore decided
  36 ;;;; to fix (A)KCL in files c/num_sfun.c and lsp/numlib.lsp. If these
  37 ;;;; fixes are applied, the feature :stat-float-is-double-float should be
  38 ;;;; defined.
  39 ;;;;
  40
  41 ;;;;
  42 ;;;; Package Setup
  43 ;;;;
  44
  45 (in-package #:lisp-stat-basics)
  46
  47 ;;(export '(+stat-float-typing+ +stat-cfloat-typing+ +stat-float-template+
  48 ;;          machine-epsilon))
  49
  50 ;; This should technically be conditionalized to the Lisp
  51 ;; implementation, i.e.
  52 #+sbcl(pushnew :stat-float-is-double-float *features*)
  53 #+cmu(pushnew :stat-float-is-double-float *features*)
  54 #+clisp(pushnew :stat-float-is-double-float *features*)
  55 ;; etc... the above need to be verified!
  56
  57 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  58 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  59 ;;;;
  60 ;;;;  Constants determining default floating point type for statistical
  61 ;;;;  operations. This type generally corresponds to a C double.
  62 ;;;;
  63 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  64 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  65
  66 ;; Why do we do this instead of defparameter?  mostly to allow for doc
  67 ;; defs as well as to minimize dereferencing mistakes.
  68 (defmacro define-constant (name value &optional doc)
  69   `(defparameter ,name (if (boundp ',name) (symbol-value ',name) ,value)
  70     ,@(when doc (list doc))))
  71
  72 ;(define-constant +stat-float-typing+ 'long-float)
  73 ;(define-constant +stat-cfloat-typing+ '(complex long-float))
  74 ;(define-constant +stat-float-template+ 0.d0)
  75
  76 (defparameter +stat-float-typing+ 'long-float)
  77 (defparameter +stat-cfloat-typing+ '(complex long-float))
  78 (defparameter +stat-float-template+ 0.d0)
  79
  80 (deftype stat-float () +stat-float-typing+)
  81 (deftype stat-cfloat () +stat-cfloat-typing+)
  82
  83 (defparameter machine-epsilon
  84   (do ((epsilon (float 1.0 +stat-float-template+) (/ epsilon 2.0)))
  85       ((= (+ 1.0 (/ epsilon 2.0)) 1.0) epsilon)))
  86
  87 (defmacro declare-double (&rest vars) `(declare (long-float ,@vars)))
  88
  89 (setf *read-default-float-format* +stat-float-typing+)
  90
  91 ;;; FIX-BASE-DOC adds note about modification to the end of the
  92 ;;; documentation string argument.
  93 (defmacro fix-base-doc (doc)
  94   `(format nil
  95     "~a~%Modified to coerce arguments(s) to stat-float or stat-cfloat."
  96     ,doc))
  97
  98 ;;; FIX-BASE-FUN-DOC fixes documentation of SYM and installs in BASE-SYM
  99 (defun fix-base-fun-doc (sym base-sym)
 100   (let ((doc (documentation sym 'function)))
 101     (if doc (setf (documentation base-sym 'function) (fix-base-doc doc)))))
 102
 103 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 104 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 105 ;;;;
 106 ;;;;  Functions and Macros for modifying functions to coerce to standard
 107 ;;;;  floating point type.
 108 ;;;;
 109 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 110 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 111
 112 ;;; FFIX - coerces its arguments to standard real or complex floating
 113 ;;; point number if needed.
 114 #-stat-float-is-double-float
 115 (defmacro ffix (x)
 116   `(if (complexp ,x)
 117     (coerce ,x ',+stat-cfloat-typing+)
 118     (float ,x ',+stat-float-template+)))
 119
 120 #+stat-float-is-double-float
 121 (defmacro ffix (x) x)
 122
 123 ;;; MAKEDOUBLE coerces its argument to the standard floating point
 124 ;;; type specified by +stat-float-template+
 125
 126 (defun makedouble (x) (float x +stat-float-template+))
 127
 128 #+stat-float-is-double-float
 129 (eval-when (compile)
 130   (proclaim '(function makedouble (t) long-float)))
 131
 132 ;;;
 133 ;;; MAKE-BASE-TRANS-FUN Macro for re-defining one argument transcendental
 134 ;;; functions
 135 ;;;
 136 #-stat-float-is-double-float
 137 (defmacro make-base-trans-fun (sym)
 138   (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
 139          (doc (documentation sym 'function))
 140          (doc-list (if doc (list (fix-base-doc doc)))))
 141   `(defun ,base-sym (x)
 142      ,@doc-list
 143      (declare (inline ,sym coerce float))
 144      (,sym (ffix x)))))
 145
 146 #+stat-float-is-double-float
 147 (defmacro make-base-trans-fun (sym)
 148   (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
 149          (doc (documentation sym 'function)))
 150     `(progn (setf (symbol-function ',base-sym) (symbol-function ',sym))
 151             (if ,doc (setf (documentation ',base-sym 'function) ,doc)))))
 152
 153 ;;;
 154 ;;; MAKE-BASE-TRANS-FUN-2 Macro for re-defining transcendental functions
 155 ;;; with an optional second argument
 156 ;;;
 157 #-stat-float-is-double-float
 158 (defmacro make-base-trans-fun-2 (sym)
 159   (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
 160          (doc (documentation sym 'function))
 161          (doc-list (if doc (list (fix-base-doc doc)))))
 162   `(defun ,base-sym (x &optional y)
 163      ,@doc-list
 164      (declare (inline ,sym coerce float))
 165      (if y (,sym (ffix x) (ffix y)) (,sym (ffix x))))))
 166
 167 #+stat-float-is-double-float
 168 (defmacro make-base-trans-fun-2 (sym)
 169   (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
 170          (doc (documentation sym 'function)))
 171     `(progn (setf (symbol-function ',base-sym) (symbol-function ',sym))
 172             (if ,doc (setf (documentation ',base-sym 'function) ,doc)))))
 173
 174 ;;; FIX-BASE-DOC adds note about modification to documentation string argument
 175 (defmacro fix-base-doc (doc)
 176   `(format nil
 177            "~a~%Modified to coerce arguments(s) to stat-float or stat-cfloat."
 178            ,doc))
 179
 180 ;;; FIX-BASE-FUN-DOC fixes documentation of SYM and installs in BASE-SYM
 181 (defun fix-base-fun-doc (sym base-sym)
 182   (let ((doc (documentation sym 'function)))
 183     (if doc (setf (documentation base-sym 'function) (fix-base-doc doc)))))
 184
 185 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 186 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 187 ;;;;
 188 ;;;;  Modified base functions to coerce to standard floating point type
 189 ;;;;
 190 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 191 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 192
 193 ;;; BASE-FLOAT
 194 #-stat-float-is-double-float
 195 (progn
 196   (defun base-float (x &optional (y +stat-float-template+)) (float x y))
 197   (fix-base-fun-doc 'float 'base-float))
 198
 199 #+stat-float-is-double-float
 200 (make-base-trans-fun float)
 201
 202 ;;; BASE-EXPT
 203 #-stat-float-is-double-float
 204 (progn
 205   (defun base-expt (x y)
 206     (declare (inline expt coerce float integerp))
 207     (if (integerp y) (expt x y) (expt (ffix x) (ffix y))))
 208
 209   (fix-base-fun-doc 'expt 'base-expt))
 210
 211 #+stat-float-is-double-float
 212 (make-base-trans-fun expt)
 213
 214 ;;; Others
 215 (make-base-trans-fun-2 log)
 216 (make-base-trans-fun exp)
 217 (make-base-trans-fun sqrt)
 218 (make-base-trans-fun sin)
 219 (make-base-trans-fun cos)
 220 (make-base-trans-fun tan)
 221 (make-base-trans-fun asin)
 222 (make-base-trans-fun acos)
 223 (make-base-trans-fun-2 atan)
 224 (make-base-trans-fun sinh)
 225 (make-base-trans-fun cosh)
 226 (make-base-trans-fun tanh)
 227 (make-base-trans-fun asinh)
 228 (make-base-trans-fun acosh)
 229 (make-base-trans-fun atanh)
 230 (make-base-trans-fun abs)
 231 (make-base-trans-fun phase)
 232 (make-base-trans-fun-2 ffloor)
 233 (make-base-trans-fun-2 fceiling)
 234 (make-base-trans-fun-2 ftruncate)
 235 (make-base-trans-fun-2 fround)
 236 (make-base-trans-fun signum)
 237 (make-base-trans-fun cis)