1 ;;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: System; Base: 10 -*-
3 ;;;; evcl - 12 - Number - float64-scale
4 ;;; arch/generic/lisp/math/gen-math-f32-scale.lisp
6 ;;; This file is part of Evita Common Lisp.
8 ;;; Copyright (C) 1996-2007 by Project Vogue.
9 ;;; Written by Yoshifumi "VOGUE" INOUE. (yosi@msn.com)
11 ;;; @(#)$Id: //proj/evcl3/mainline/arch/generic/lisp/libm/float64/gen-float64-scale.lisp#1 $
14 ;;; This file contains implementation of float64-scale.
16 (in-package #:strict-math)
19 * From fdlibm (http://www.netlib.org/fdlibm/)
20 /* @(#)s_scalbn.c 5.1 93/09/24 */
22 * ====================================================
23 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
29 * ====================================================
34 <<scalbn>>, <<scalbnf>>---scale by power of two
42 double scalbn(double <[x]>, int <[y]>);
43 float scalbnf(float <[x]>, int <[y]>);
47 double scalbn(<[x]>,<[y]>)
50 float scalbnf(<[x]>,<[y]>)
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
61 <[x]> times 2 to the power <[n]>.
64 Neither <<scalbn>> nor <<scalbnf>> is required by ANSI C or by the System V
65 Interface Definition (Issue 2).
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.
76 (defun scale-float64 (x n)
77 (declare (values double-float))
78 (declare (type double-float x))
79 (declare (type fixnum n))
81 (two54 #+nil 1.80143985094819840000e+16
82 #.(encode-float64 #x43500000 #x00000000) )
83 (twom54 #+nil 5.55111512312578270212e-17
84 #.(encode-float64 #x3C900000 #x00000000) )
86 #.(encode-float64 #x7E37E43C #x8800759B) )
88 #.(encode-float64 #x01A56E1F #xC2F8F359) )
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)
97 (multiple-value-setq (hx lx) (decode-float64 x))
98 (setq k (- (ash (logand hx #x7ff00000) -20) 54))
100 (when (< n -50000) (return (* tiny x))) )
102 (when (eql k #x7ff) (return (+ x x))) ; NaN or Inf
108 (return (* huge (float64-sign huge x))) )
110 (when (> k 0) ; normal result
111 (return (encode-float64
112 (logior (logand hx #x800fffff) (ash k 20)) lx )))
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)) ) ) ) ) )