1*05a0b428SJohn Marino /* $OpenBSD: s_scalbnl.c,v 1.3 2013/11/12 18:28:02 martynas Exp $ */ 2*05a0b428SJohn Marino /* @(#)s_scalbn.c 5.1 93/09/24 */ 3*05a0b428SJohn Marino /* 4*05a0b428SJohn Marino * ==================================================== 5*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6*05a0b428SJohn Marino * 7*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business. 8*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this 9*05a0b428SJohn Marino * software is freely granted, provided that this notice 10*05a0b428SJohn Marino * is preserved. 11*05a0b428SJohn Marino * ==================================================== 12*05a0b428SJohn Marino */ 13*05a0b428SJohn Marino 14*05a0b428SJohn Marino /* 15*05a0b428SJohn Marino * scalbnl (long double x, int n) 16*05a0b428SJohn Marino * scalbnl(x,n) returns x* 2**n computed by exponent 17*05a0b428SJohn Marino * manipulation rather than by actually performing an 18*05a0b428SJohn Marino * exponentiation or a multiplication. 19*05a0b428SJohn Marino */ 20*05a0b428SJohn Marino 21*05a0b428SJohn Marino /* 22*05a0b428SJohn Marino * We assume that a long double has a 15-bit exponent. On systems 23*05a0b428SJohn Marino * where long double is the same as double, scalbnl() is an alias 24*05a0b428SJohn Marino * for scalbn(), so we don't use this routine. 25*05a0b428SJohn Marino */ 26*05a0b428SJohn Marino 27*05a0b428SJohn Marino #include <sys/types.h> 28*05a0b428SJohn Marino #include <machine/ieee.h> 29*05a0b428SJohn Marino #include <float.h> 30*05a0b428SJohn Marino #include <math.h> 31*05a0b428SJohn Marino 32*05a0b428SJohn Marino #if LDBL_MAX_EXP != 0x4000 33*05a0b428SJohn Marino #error "Unsupported long double format" 34*05a0b428SJohn Marino #endif 35*05a0b428SJohn Marino 36*05a0b428SJohn Marino static const long double 37*05a0b428SJohn Marino huge = 0x1p16000L, 38*05a0b428SJohn Marino tiny = 0x1p-16000L; 39*05a0b428SJohn Marino 40*05a0b428SJohn Marino long double 41*05a0b428SJohn Marino scalbnl (long double x, int n) 42*05a0b428SJohn Marino { 43*05a0b428SJohn Marino union { 44*05a0b428SJohn Marino long double e; 45*05a0b428SJohn Marino struct ieee_ext bits; 46*05a0b428SJohn Marino } u; 47*05a0b428SJohn Marino int k; 48*05a0b428SJohn Marino u.e = x; 49*05a0b428SJohn Marino k = u.bits.ext_exp; /* extract exponent */ 50*05a0b428SJohn Marino if (k==0) { /* 0 or subnormal x */ 51*05a0b428SJohn Marino if ((u.bits.ext_frach 52*05a0b428SJohn Marino #ifdef EXT_FRACHMBITS 53*05a0b428SJohn Marino | u.bits.ext_frachm 54*05a0b428SJohn Marino #endif /* EXT_FRACHMBITS */ 55*05a0b428SJohn Marino #ifdef EXT_FRACLMBITS 56*05a0b428SJohn Marino | u.bits.ext_fraclm 57*05a0b428SJohn Marino #endif /* EXT_FRACLMBITS */ 58*05a0b428SJohn Marino | u.bits.ext_fracl)==0) return x; /* +-0 */ 59*05a0b428SJohn Marino u.e *= 0x1p+128; 60*05a0b428SJohn Marino k = u.bits.ext_exp - 128; 61*05a0b428SJohn Marino if (n< -50000) return tiny*x; /*underflow*/ 62*05a0b428SJohn Marino } 63*05a0b428SJohn Marino if (k==0x7fff) return x+x; /* NaN or Inf */ 64*05a0b428SJohn Marino k = k+n; 65*05a0b428SJohn Marino if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow */ 66*05a0b428SJohn Marino if (k > 0) /* normal result */ 67*05a0b428SJohn Marino {u.bits.ext_exp = k; return u.e;} 68*05a0b428SJohn Marino if (k <= -128) { 69*05a0b428SJohn Marino if (n > 50000) /* in case integer overflow in n+k */ 70*05a0b428SJohn Marino return huge*copysign(huge,x); /*overflow*/ 71*05a0b428SJohn Marino else return tiny*copysign(tiny,x); /*underflow*/ 72*05a0b428SJohn Marino } 73*05a0b428SJohn Marino k += 128; /* subnormal result */ 74*05a0b428SJohn Marino u.bits.ext_exp = k; 75*05a0b428SJohn Marino return u.e*0x1p-128; 76*05a0b428SJohn Marino } 77*05a0b428SJohn Marino 78*05a0b428SJohn Marino long double 79*05a0b428SJohn Marino ldexpl(long double x, int n) 80*05a0b428SJohn Marino { 81*05a0b428SJohn Marino return scalbnl(x, n); 82*05a0b428SJohn Marino } 83