code.consxy.com Git - clnl/blob - resources/strictmath/otherconversions/gen-float64-scale.lisp

   1 ;;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: System; Base: 10 -*-
   2 ;;;;
   3 ;;;; evcl - 12 - Number - float64-scale
   4 ;;; arch/generic/lisp/math/gen-math-f32-scale.lisp
   5 ;;;
   6 ;;; This file is part of Evita Common Lisp.
   7 ;;;
   8 ;;; Copyright (C) 1996-2007 by Project Vogue.
   9 ;;; Written by Yoshifumi "VOGUE" INOUE. (yosi@msn.com)
  10 ;;;
  11 ;;; @(#)$Id: //proj/evcl3/mainline/arch/generic/lisp/libm/float64/gen-float64-scale.lisp#1 $
  12 ;;;
  13 ;;; Description:
  14 ;;;  This file contains implementation of float64-scale.
  15 ;
  16 (in-package #:strict-math)
  17
  18 #|
  19  * From fdlibm (http://www.netlib.org/fdlibm/)
  20 /* @(#)s_scalbn.c 5.1 93/09/24 */
  21 /*
  22  * ====================================================
  23  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  24  *
  25  * Developed at SunPro, a Sun Microsystems, Inc. business.
  26  * Permission to use, copy, modify, and distribute this
  27  * software is freely granted, provided that this notice
  28  * is preserved.
  29  * ====================================================
  30  */
  31
  32 /*
  33 FUNCTION
  34 <<scalbn>>, <<scalbnf>>---scale by power of two
  35 INDEX
  36         scalbn
  37 INDEX
  38         scalbnf
  39
  40 ANSI_SYNOPSIS
  41         #include <math.h>
  42         double scalbn(double <[x]>, int <[y]>);
  43         float scalbnf(float <[x]>, int <[y]>);
  44
  45 TRAD_SYNOPSIS
  46         #include <math.h>
  47         double scalbn(<[x]>,<[y]>)
  48         double <[x]>;
  49         int <[y]>;
  50         float scalbnf(<[x]>,<[y]>)
  51         float <[x]>;
  52         int <[y]>;
  53
  54 DESCRIPTION
  55 <<scalbn>> and <<scalbnf>> scale <[x]> by <[n]>, returning <[x]> times
  56 2 to the power <[n]>.  The result is computed by manipulating the
  57 exponent, rather than by actually performing an exponentiation or
  58 multiplication.
  59
  60 RETURNS
  61 <[x]> times 2 to the power <[n]>.
  62
  63 PORTABILITY
  64 Neither <<scalbn>> nor <<scalbnf>> is required by ANSI C or by the System V
  65 Interface Definition (Issue 2).
  66
  67 */
  68
  69 /*
  70  * scalbn (double x, int n)
  71  * scalbn(x,n) returns x* 2**n  computed by  exponent
  72  * manipulation rather than by actually performing an
  73  * exponentiation or a multiplication.
  74  */
  75 |#
  76 (defun scale-float64 (x n)
  77     (declare (values double-float))
  78     (declare (type double-float x))
  79     (declare (type fixnum n))
  80   (prog* (
  81     (two54  #+nil 1.80143985094819840000e+16
  82             #.(encode-float64 #x43500000 #x00000000) )
  83     (twom54 #+nil 5.55111512312578270212e-17
  84             #.(encode-float64 #x3C900000 #x00000000) )
  85     (huge   #+nil 1.0e+300
  86             #.(encode-float64 #x7E37E43C #x8800759B) )
  87     (tiny   #+nil 1.0e-300
  88             #.(encode-float64 #x01A56E1F #xC2F8F359) )
  89     )
  90     ;;
  91     (multiple-value-bind (hx lx) (decode-float64 x)
  92     (let ((k (ash (logand hx #x7ff00000) -20))) ; extract exponent
  93       (when (eql k 0)   ; 0 or subnormal x
  94         (when (eql (logior lx (logand hx #x7fffffff)) 0)
  95           (return x) )  ; +-0
  96         (setq x (* x two54))
  97         (multiple-value-setq (hx lx) (decode-float64 x))
  98         (setq k (- (ash (logand hx #x7ff00000) -20) 54))
  99         ;; underflow
 100         (when (< n -50000) (return (* tiny x))) )
 101
 102       (when (eql k #x7ff) (return (+ x x))) ; NaN or Inf
 103
 104       (incf k n)
 105
 106       (when (> k  #x7fe)
 107         ;; overflow
 108         (return (* huge (float64-sign huge x))) )
 109
 110       (when (> k  0)  ; normal result
 111         (return (encode-float64
 112             (logior (logand hx #x800fffff) (ash k 20)) lx )))
 113
 114       (when (<= k -54)
 115         (return  (if (> n  50000)    ; in case integer overflow in n+k
 116             (* huge (float64-sign huge x))     ; *overflow*
 117             (* tiny (float64-sign tiny x)) ))) ; *underflow*
 118       (incf k 54)   ; subnormal result
 119       (let ((x (encode-float64
 120                   (logior (logand hx #x800fffff) (ash k 20)) lx) ))
 121         (return (* x twom54)) ) ) ) ) )