Return to s_modf.c CVS log

Up to [cvs.NetBSD.org] / src / lib / libm / src

Annotation of src/lib/libm/src/s_modf.c, Revision 1.2

1.1       jtc         1:
                      2: /* @(#)s_modf.c 5.1 93/09/24 */
                      3: /*
                      4:  * ====================================================
                      5:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
                      6:  *
                      7:  * Developed at SunPro, a Sun Microsystems, Inc. business.
                      8:  * Permission to use, copy, modify, and distribute this
                      9:  * software is freely granted, provided that this notice
                     10:  * is preserved.
                     11:  * ====================================================
                     12:  */
                     13:
                     14: /*
                     15:  * modf(double x, double *iptr)
                     16:  * return fraction part of x, and return x's integral part in *iptr.
                     17:  * Method:
                     18:  *     Bit twiddling.
                     19:  *
                     20:  * Exception:
                     21:  *     No exception.
                     22:  */
                     23:
1.2     ! jtc        24: #include <math.h>
1.1       jtc        25:
                     26: #ifdef __STDC__
                     27: static const double one = 1.0;
                     28: #else
                     29: static double one = 1.0;
                     30: #endif
                     31:
                     32: #ifdef __STDC__
                     33:        double modf(double x, double *iptr)
                     34: #else
                     35:        double modf(x, iptr)
                     36:        double x,*iptr;
                     37: #endif
                     38: {
                     39:        int i0,i1,n0,n1,j0;
                     40:        unsigned i;
                     41:        n0 = (*((int *)&one)>>29)^1;    /* high word index */
                     42:        n1 = 1-n0;
                     43:        i0 =  *(n0+(int*)&x);           /* high x */
                     44:        i1 =  *(n1+(int*)&x);           /* low  x */
                     45:        j0 = ((i0>>20)&0x7ff)-0x3ff;    /* exponent of x */
                     46:        if(j0<20) {                     /* integer part in high x */
                     47:            if(j0<0) {                  /* |x|<1 */
                     48:                *(n0+(int*)iptr) = i0&0x80000000;
                     49:                *(n1+(int*)iptr) = 0;           /* *iptr = +-0 */
                     50:                return x;
                     51:            } else {
                     52:                i = (0x000fffff)>>j0;
                     53:                if(((i0&i)|i1)==0) {            /* x is integral */
                     54:                    *iptr = x;
                     55:                    *(n0+(int*)&x) &= 0x80000000;
                     56:                    *(n1+(int*)&x)  = 0;        /* return +-0 */
                     57:                    return x;
                     58:                } else {
                     59:                    *(n0+(int*)iptr) = i0&(~i);
                     60:                    *(n1+(int*)iptr) = 0;
                     61:                    return x - *iptr;
                     62:                }
                     63:            }
                     64:        } else if (j0>51) {             /* no fraction part */
                     65:            *iptr = x*one;
                     66:            *(n0+(int*)&x) &= 0x80000000;
                     67:            *(n1+(int*)&x)  = 0;        /* return +-0 */
                     68:            return x;
                     69:        } else {                        /* fraction part in low x */
                     70:            i = ((unsigned)(0xffffffff))>>(j0-20);
                     71:            if((i1&i)==0) {             /* x is integral */
                     72:                *iptr = x;
                     73:                *(n0+(int*)&x) &= 0x80000000;
                     74:                *(n1+(int*)&x)  = 0;    /* return +-0 */
                     75:                return x;
                     76:            } else {
                     77:                *(n0+(int*)iptr) = i0;
                     78:                *(n1+(int*)iptr) = i1&(~i);
                     79:                return x - *iptr;
                     80:            }
                     81:        }
                     82: }