/* $Header:scalb.c 12.0$ */
/* $ACIS:scalb.c 12.0$ */
/* $Source: /ibm/acis/usr/src/lib/libc/ca/gen/RCS/scalb.c,v $ */

#ifndef lint
static char *rcsid = "$Header:scalb.c 12.0$";
#endif

/* 
 * Copyright (c) 1985 Regents of the University of California.
 * 
 * Use and reproduction of this software are granted  in  accordance  with
 * the terms and conditions specified in  the  Berkeley  Software  License
 * Agreement (in particular, this entails acknowledgement of the programs'
 * source, and inclusion of this notice) with the additional understanding
 * that  all  recipients  should regard themselves as participants  in  an
 * ongoing  research  project and hence should  feel  obligated  to report
 * their  experiences (good or bad) with these elementary function  codes,
 * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
 */

/*
 *  scalb(x,N) 
 *              returns  x * (2**N), for integer values N.
 *  logb(x) 
 *              returns the unbiased exponent of x, a signed integer in 
 *              double precision, except that logb(0) is -INF, logb(INF) 
 *              is +INF, and logb(NAN) is that NAN.
 */

#include <errno.h>
extern errno;

static unsigned short msign=0x7fff, mexp =0x7ff0  ;
static short prep1=54, gap=4, bias=1023           ;
static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
extern double copysign();

#ifdef RTFL
double _scalb(x,N)
#else
double scalb(x,N)
#endif RTFL

double x; int N;
{
        int k;
        double scalb();

        unsigned short *px=(unsigned short *) &x;

        if( (k= *px & mexp ) != mexp ) {
            if( x == zero )  return(x); 
            if( N<-2100) 
		return((copysign(nunf,x))*nunf); 
	    else if(N>2100) {
		errno = ERANGE; return(novf*(copysign(2.0,x))); }
            if( k == 0 ) {
                 x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}

            if((k = (k>>gap)+ N) > 0 )
                if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
                else {					/* overflow */
			errno = ERANGE; x=(novf*(copysign(2.0,x))); }
            else
                if( k > -prep1 ) 
                                        /* gradual underflow */
                    {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
                else
		return((copysign(nunf,x))*nunf); 
            }
        return(x);
}

#define IEEE854
#if defined(IEEE854) && defined(IEEE754)
#undef IEEE754
#endif

#ifdef IEEE754
#ifdef RTFL
double _logb(x)
#else
double logb(x)
#endif RTFL

double x; 
{

        short *px=(short *) &x, k;

        if( (k= *px & mexp ) != mexp )
            if ( k != 0 )
                return ( (k>>gap) - bias );
            else if( x != zero)
                return ( -1022.0 );
            else        
                return(-1.0/zero);    
        else if (isnan(x))
            return(x);
        else
            {*px &= msign; return(x);}
}
#endif IEEE754

#ifdef IEEE854
#ifdef RTFL
double _logb(x)
#else
double logb(x)
#endif RTFL

double x; 

{
#define MSIGN	0x7fffffff
#define MEXP	0x7ff00000
#define GAP	20
#define BIAS	1023


	unsigned int *px, uhi, ulo, k;
	int i;
	
	px = (unsigned int *) &x;
	uhi = *px;
	ulo = *(px + 1);
	if ( (i = uhi & MEXP) != MEXP ) {
		if ( i != 0 )			/* normalized */
			return ( (i>>GAP) - BIAS );
		k = uhi & 0x000fffff;
		if ( (k) || (ulo) ) {		/* denormalized */
			if (k)			/* check hi 20 bit */
				i = (-1022-20);
			else			/* check lo 32 bit */
				{ k = ulo; i = (-1022-20-32); }
			while (k >>= 1) i++;
			return (i);
		} else				/* zero */
			return(-1.0/zero);	/* signal div by 0 */
	} else if (isnan(x))
		return(x);
	else					/* +INF or -INF */
		{*px &= MSIGN; return(x);}
}
#endif IEEE854
