1*05a0b428SJohn Marino /* @(#)s_floor.c 5.1 93/09/24 */ 2*05a0b428SJohn Marino /* 3*05a0b428SJohn Marino * ==================================================== 4*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5*05a0b428SJohn Marino * 6*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business. 7*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this 8*05a0b428SJohn Marino * software is freely granted, provided that this notice 9*05a0b428SJohn Marino * is preserved. 10*05a0b428SJohn Marino * ==================================================== 11*05a0b428SJohn Marino */ 12*05a0b428SJohn Marino 13*05a0b428SJohn Marino /* 14*05a0b428SJohn Marino * floorl(x) 15*05a0b428SJohn Marino * Return x rounded toward -inf to integral value 16*05a0b428SJohn Marino * Method: 17*05a0b428SJohn Marino * Bit twiddling. 18*05a0b428SJohn Marino * Exception: 19*05a0b428SJohn Marino * Inexact flag raised if x not equal to floor(x). 20*05a0b428SJohn Marino */ 21*05a0b428SJohn Marino 22*05a0b428SJohn Marino #include <math.h> 23*05a0b428SJohn Marino 24*05a0b428SJohn Marino #include "math_private.h" 25*05a0b428SJohn Marino 26*05a0b428SJohn Marino static const long double huge = 1.0e4930L; 27*05a0b428SJohn Marino 28*05a0b428SJohn Marino long double 29*05a0b428SJohn Marino floorl(long double x) 30*05a0b428SJohn Marino { 31*05a0b428SJohn Marino int32_t i1,jj0; 32*05a0b428SJohn Marino u_int32_t i,j,se,i0,sx; 33*05a0b428SJohn Marino GET_LDOUBLE_WORDS(se,i0,i1,x); 34*05a0b428SJohn Marino sx = (se>>15)&1; 35*05a0b428SJohn Marino jj0 = (se&0x7fff)-0x3fff; 36*05a0b428SJohn Marino if(jj0<31) { 37*05a0b428SJohn Marino if(jj0<0) { /* raise inexact if x != 0 */ 38*05a0b428SJohn Marino if(huge+x>0.0) { 39*05a0b428SJohn Marino if(sx==0) 40*05a0b428SJohn Marino return 0.0L; 41*05a0b428SJohn Marino else if(((se&0x7fff)|i0|i1)!=0) 42*05a0b428SJohn Marino return -1.0L; 43*05a0b428SJohn Marino } 44*05a0b428SJohn Marino } else { 45*05a0b428SJohn Marino i = (0x7fffffff)>>jj0; 46*05a0b428SJohn Marino if(((i0&i)|i1)==0) return x; /* x is integral */ 47*05a0b428SJohn Marino if(huge+x>0.0) { /* raise inexact flag */ 48*05a0b428SJohn Marino if(sx) { 49*05a0b428SJohn Marino if (jj0>0 && (i0+(0x80000000>>jj0))>i0) 50*05a0b428SJohn Marino i0 += (0x80000000)>>jj0; 51*05a0b428SJohn Marino else 52*05a0b428SJohn Marino { 53*05a0b428SJohn Marino i = 0x7fffffff; 54*05a0b428SJohn Marino ++se; 55*05a0b428SJohn Marino } 56*05a0b428SJohn Marino } 57*05a0b428SJohn Marino i0 &= (~i); i1=0; 58*05a0b428SJohn Marino } 59*05a0b428SJohn Marino } 60*05a0b428SJohn Marino } else if (jj0>62) { 61*05a0b428SJohn Marino if(jj0==0x4000) return x+x; /* inf or NaN */ 62*05a0b428SJohn Marino else return x; /* x is integral */ 63*05a0b428SJohn Marino } else { 64*05a0b428SJohn Marino i = ((u_int32_t)(0xffffffff))>>(jj0-31); 65*05a0b428SJohn Marino if((i1&i)==0) return x; /* x is integral */ 66*05a0b428SJohn Marino if(huge+x>0.0) { /* raise inexact flag */ 67*05a0b428SJohn Marino if(sx) { 68*05a0b428SJohn Marino if(jj0==31) i0+=1; 69*05a0b428SJohn Marino else { 70*05a0b428SJohn Marino j = i1+(1<<(63-jj0)); 71*05a0b428SJohn Marino if(j<i1) i0 +=1 ; /* got a carry */ 72*05a0b428SJohn Marino i1=j; 73*05a0b428SJohn Marino } 74*05a0b428SJohn Marino } 75*05a0b428SJohn Marino i1 &= (~i); 76*05a0b428SJohn Marino } 77*05a0b428SJohn Marino } 78*05a0b428SJohn Marino SET_LDOUBLE_WORDS(x,se,i0,i1); 79*05a0b428SJohn Marino return x; 80*05a0b428SJohn Marino } 81