external/oct/qd-dd.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 ;;; double-double float routines needed for quad-double.
  29 ;;;
  30 ;;; Not needed for CMUCL.
  31 ;;;
  32 ;;; These routines were taken directly from CMUCL.
  33
  34 (declaim (inline quick-two-sum))
  35 (defun quick-two-sum (a b)
  36   "Computes fl(a+b) and err(a+b), assuming |a| >= |b|"
  37   (declare (double-float a b))
  38   (let* ((s (+ a b))
  39          (e (- b (- s a))))
  40     (values s e)))
  41
  42 (declaim (inline two-sum))
  43 (defun two-sum (a b)
  44   "Computes fl(a+b) and err(a+b)"
  45   (declare (double-float a b))
  46   (let* ((s (+ a b))
  47          (v (- s a))
  48          (e (+ (- a (- s v))
  49                (- b v))))
  50     (values s e)))
  51
  52 (declaim (inline two-prod))
  53 (declaim (inline split))
  54 ;; This algorithm is the version given by Yozo Hida.  It has problems
  55 ;; with overflow because we multiply by 1+2^27.
  56 ;;
  57 ;; But be very careful about replacing this with a new algorithm.  The
  58 ;; values computed here are very important to get the rounding right.
  59 ;; If you change this, the rounding may be different, which will
  60 ;; affect other parts of the algorithm.
  61 ;;
  62 ;; I (rtoy) tried a different algorithm that split the number in two
  63 ;; as described, but without overflow.  However, that caused
  64 ;; -9.4294948327242751340284975915175w0/1w14 to return a value that
  65 ;; wasn't really close to -9.4294948327242751340284975915175w-14.
  66 ;;
  67 ;; This also means we can't print numbers like 1w308 with the current
  68 ;; printing algorithm, or even divide 1w308 by 10.
  69 #+nil
  70 (defun split (a)
  71   "Split the double-float number a into a-hi and a-lo such that a =
  72   a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
  73   a-lo contains the lower 26 bits."
  74   (declare (double-float a))
  75   (let* ((tmp (* a (+ 1 (expt 2 27))))
  76          (a-hi (- tmp (- tmp a)))
  77          (a-lo (- a a-hi)))
  78     (values a-hi a-lo)))
  79
  80
  81 (defun split (a)
  82   "Split the double-float number a into a-hi and a-lo such that a =
  83   a-hi + a-lo and a-hi contains the upper 26 significant bits of a and
  84   a-lo contains the lower 26 bits."
  85   (declare (double-float a)
  86            (optimize (speed 3)))
  87   ;; This splits the number a into 2 halves of 26 bits each, but the
  88   ;; halves are, I think, supposed to be properly rounded in an IEEE
  89   ;; fashion.
  90   ;;
  91   ;; For numbers that are very large, we use a different algorithm.
  92   ;; For smaller numbers, we can use the original algorithm of Yozo
  93   ;; Hida.
  94   (if (> (abs a) (scale-float 1d0 (- 1023 27)))
  95       ;; I've tested this algorithm against Yozo's method for 1
  96       ;; billion randomly generated double-floats between 2^(-995) and
  97       ;; 2^996, and identical results are obtained.  For numbers that
  98       ;; are very small, this algorithm produces different numbers
  99       ;; because of underflow.  For very large numbers, we, of course
 100       ;; produce different results because Yozo's method causes
 101       ;; overflow.
 102       (let* ((tmp (* a (+ 1 (scale-float 1d0 -27))))
 103              (as (* a (scale-float 1d0 -27)))
 104              (a-hi (* (- tmp (- tmp as)) (expt 2 27)))
 105              (a-lo (- a a-hi)))
 106         (values a-hi a-lo))
 107       ;; Yozo's algorithm.
 108       (let* ((tmp (* a (+ 1 (expt 2 27))))
 109              (a-hi (- tmp (- tmp a)))
 110              (a-lo (- a a-hi)))
 111         (values a-hi a-lo))))
 112
 113
 114 (defun two-prod (a b)
 115   "Compute fl(a*b) and err(a*b)"
 116   (declare (double-float a b))
 117   (let ((p (* a b)))
 118     (multiple-value-bind (a-hi a-lo)
 119         (split a)
 120       (multiple-value-bind (b-hi b-lo)
 121           (split b)
 122         (let ((e (+ (+ (- (* a-hi b-hi) p)
 123                        (* a-hi b-lo)
 124                        (* a-lo b-hi))
 125                     (* a-lo b-lo))))
 126           (values p e))))))
 127
 128 (declaim (inline two-sqr))
 129 (defun two-sqr (a)
 130   "Compute fl(a*a) and err(a*b).  This is a more efficient
 131   implementation of two-prod"
 132   (declare (double-float a))
 133   (let ((q (* a a)))
 134     (multiple-value-bind (a-hi a-lo)
 135         (split a)
 136       (values q (+ (+ (- (* a-hi a-hi) q)
 137                       (* 2 a-hi a-lo))
 138                    (* a-lo a-lo))))))