external/oct/qd-rep.lisp

   1 ;;;; -*- Mode: lisp -*-
   2 ;;;;
   3 ;;;; Copyright (c) 2007 Raymond Toy
   4 ;;;;
   5 ;;;; Permission is hereby granted, free of charge, to any person
   6 ;;;; obtaining a copy of this software and associated documentation
   7 ;;;; files (the "Software"), to deal in the Software without
   8 ;;;; restriction, including without limitation the rights to use,
   9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
  10 ;;;; copies of the Software, and to permit persons to whom the
  11 ;;;; Software is furnished to do so, subject to the following
  12 ;;;; conditions:
  13 ;;;;
  14 ;;;; The above copyright notice and this permission notice shall be
  15 ;;;; included in all copies or substantial portions of the Software.
  16 ;;;;
  17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
  25
  26 (in-package #:qdi)
  27
  28 ;;; This file contains the actual representation of a %quad-double
  29 ;;; number.  The only real requirement for a %quad-double number is an
  30 ;;; object that can hold four double-float values.
  31 ;;;
  32 ;;; This object is created by %MAKE-QD-D.  The four double-float
  33 ;;; elements of a %quad-double are accessed via QD-0, QD-1, QD-2, and
  34 ;;; QD-3.  A convenience function, QD-PARTS, is also provided to
  35 ;;; return all four values at once.
  36
  37 ;; All of the following functions should be inline to reduce consing.
  38 (declaim (inline
  39           qd-0 qd-1 qd-2 qd-3
  40           %make-qd-d
  41           qd-parts))
  42 #+cmu
  43 (progn
  44 ;; For CMUCL (at least recent enough versions that support
  45 ;; double-double-float), we can use a (complex double-double-float) to
  46 ;; hold our 4 double-float values.  This has a nice advantage: Much of
  47 ;; the arithmetic can be done without consing, provided the key
  48 ;; functions are inline.
  49 (deftype %quad-double ()
  50   '(complex double-double-float))
  51
  52 ;; QD-0, QD-1, QD-2, and QD-3 extract the various parts of a
  53 ;; quad-double.  QD-0 is the most significant part and QD-3 is the
  54 ;; least.
  55 (defun qd-0 (q)
  56   (declare (type %quad-double q)
  57            (optimize (speed 3)))
  58   (kernel:double-double-hi (realpart q)))
  59 (defun qd-1 (q)
  60   (declare (type %quad-double q)
  61            (optimize (speed 3)))
  62   (kernel:double-double-lo (realpart q)))
  63 (defun qd-2 (q)
  64   (declare (type %quad-double q)
  65            (optimize (speed 3)))
  66   (kernel:double-double-hi (imagpart q)))
  67 (defun qd-3 (q)
  68   (declare (type %quad-double q)
  69            (optimize (speed 3)))
  70   (kernel:double-double-lo (imagpart q)))
  71
  72 (eval-when (:compile-toplevel :load-toplevel :execute)
  73 (defun %make-qd-d (a0 a1 a2 a3)
  74   "Make a %quad-double from 4 double-floats, exactly using the given
  75   values.  No check is made to see if the values make sense.  A0 is
  76   the most significant part and A3 is the least.
  77 "
  78   (declare (double-float a0 a1
  79                          a2 a3))
  80   (complex (kernel:%make-double-double-float a0 a1)
  81            (kernel:%make-double-double-float a2 a3)))
  82 )
  83
  84
  85 (defun qd-parts (qd)
  86   "Extract the four doubles comprising a quad-double and return them
  87   as multiple values.  The most significant double is the first value."
  88   (declare (type %quad-double qd))
  89   (let ((re (realpart qd))
  90         (im (imagpart qd)))
  91     (values (kernel:double-double-hi re)
  92             (kernel:double-double-lo re)
  93             (kernel:double-double-hi im)
  94             (kernel:double-double-lo im))))
  95
  96 ) ; end progn
  97
  98 #-cmu
  99 (progn
 100 ;; For Lisp's without a double-double-float type, I think the best we
 101 ;; can do is a simple-array of four double-floats.  Even with
 102 ;; inlining, I think there will lots of consing when working with this
 103 ;; type.
 104 ;;
 105 ;; A defstruct would also work but I think a simple-array is the
 106 ;; simplest and smallest representation.
 107 (deftype %quad-double ()
 108   '(simple-array double-float (4)))
 109
 110 (defun qd-0 (q)
 111   (declare (type %quad-double q)
 112            (optimize (speed 3)))
 113   (aref q 0))
 114
 115 (defun qd-1 (q)
 116   (declare (type %quad-double q)
 117            (optimize (speed 3)))
 118   (aref q 1))
 119
 120 (defun qd-2 (q)
 121   (declare (type %quad-double q)
 122            (optimize (speed 3)))
 123   (aref q 2))
 124
 125 (defun qd-3 (q)
 126   (declare (type %quad-double q)
 127            (optimize (speed 3)))
 128   (aref q 3))
 129
 130 (eval-when (:compile-toplevel :load-toplevel :execute)
 131 (defun %make-qd-d (a0 a1 a2 a3)
 132   "Make a %quad-double from 4 double-floats, exactly using the given
 133   values.  No check is made to see if the values make sense.  A0 is
 134   the most significant part and A3 is the least.
 135 "
 136   (declare (double-float a0 a1
 137                          a2 a3)
 138            (optimize (speed 3)))
 139   (let ((a (make-array 4 :element-type 'double-float)))
 140     (setf (aref a 0) a0)
 141     (setf (aref a 1) a1)
 142     (setf (aref a 2) a2)
 143     (setf (aref a 3) a3)
 144     a))
 145 )
 146
 147 (defun qd-parts (qd)
 148   "Extract the four doubles comprising a quad-double and return them
 149   as multiple values.  The most significant double is the first value."
 150   (declare (type %quad-double qd))
 151   (values (aref qd 0)
 152           (aref qd 1)
 153           (aref qd 2)
 154           (aref qd 3)))
 155
 156 ) ; end progn