12116Sjkh/* @(#)s_scalbn.c 5.1 93/09/24 */ 22116Sjkh/* 32116Sjkh * ==================================================== 42116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 52116Sjkh * 62116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 72116Sjkh * Permission to use, copy, modify, and distribute this 88870Srgrimes * software is freely granted, provided that this notice 92116Sjkh * is preserved. 102116Sjkh * ==================================================== 112116Sjkh */ 122116Sjkh 13324006Sdim#include <sys/cdefs.h> 14324006Sdim__FBSDID("$FreeBSD: stable/11/lib/msun/src/s_scalbn.c 324006 2017-09-26 09:01:56Z dim $"); 152116Sjkh 168870Srgrimes/* 172116Sjkh * scalbn (double x, int n) 188870Srgrimes * scalbn(x,n) returns x* 2**n computed by exponent 198870Srgrimes * manipulation rather than by actually performing an 202116Sjkh * exponentiation or a multiplication. 212116Sjkh */ 222116Sjkh 23143220Sdas#include <float.h> 24143220Sdas 252116Sjkh#include "math.h" 262116Sjkh#include "math_private.h" 272116Sjkh 282116Sjkhstatic const double 292116Sjkhtwo54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ 302116Sjkhtwom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ 312116Sjkhhuge = 1.0e+300, 322116Sjkhtiny = 1.0e-300; 332116Sjkh 3497413Salfreddouble 35117912Speterscalbn (double x, int n) 362116Sjkh{ 372116Sjkh int32_t k,hx,lx; 382116Sjkh EXTRACT_WORDS(hx,lx,x); 392116Sjkh k = (hx&0x7ff00000)>>20; /* extract exponent */ 402116Sjkh if (k==0) { /* 0 or subnormal x */ 412116Sjkh if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ 428870Srgrimes x *= two54; 432116Sjkh GET_HIGH_WORD(hx,x); 448870Srgrimes k = ((hx&0x7ff00000)>>20) - 54; 452116Sjkh if (n< -50000) return tiny*x; /*underflow*/ 462116Sjkh } 472116Sjkh if (k==0x7ff) return x+x; /* NaN or Inf */ 488870Srgrimes k = k+n; 492116Sjkh if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ 502116Sjkh if (k > 0) /* normal result */ 512116Sjkh {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} 52324006Sdim if (k <= -54) { 532116Sjkh if (n > 50000) /* in case integer overflow in n+k */ 542116Sjkh return huge*copysign(huge,x); /*overflow*/ 55324006Sdim else 56324006Sdim return tiny*copysign(tiny,x); /*underflow*/ 57324006Sdim } 582116Sjkh k += 54; /* subnormal result */ 592116Sjkh SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); 602116Sjkh return x*twom54; 612116Sjkh} 62143220Sdas 63143220Sdas#if (LDBL_MANT_DIG == 53) 64143264Sdas__weak_reference(scalbn, ldexpl); 65143264Sdas__weak_reference(scalbn, scalbnl); 66143220Sdas#endif 67