Return to e_remainder.c CVS log

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

Annotation of src/lib/libm/src/e_remainder.c, Revision 1.4

1.1       jtc         1: /* @(#)e_remainder.c 5.1 93/09/24 */
                      2: /*
                      3:  * ====================================================
                      4:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
                      5:  *
                      6:  * Developed at SunPro, a Sun Microsystems, Inc. business.
                      7:  * Permission to use, copy, modify, and distribute this
                      8:  * software is freely granted, provided that this notice
                      9:  * is preserved.
                     10:  * ====================================================
                     11:  */
1.3       jtc        12:
                     13: #ifndef lint
1.4     ! jtc        14: static char rcsid[] = "$Id: e_remainder.c,v 1.3 1994/02/18 02:25:44 jtc Exp $";
1.3       jtc        15: #endif
1.1       jtc        16:
                     17: /* __ieee754_remainder(x,p)
                     18:  * Return :
                     19:  *     returns  x REM p  =  x - [x/p]*p as if in infinite
                     20:  *     precise arithmetic, where [x/p] is the (infinite bit)
                     21:  *     integer nearest x/p (in half way case choose the even one).
                     22:  * Method :
                     23:  *     Based on fmod() return x-[x/p]chopped*p exactlp.
                     24:  */
                     25:
1.2       jtc        26: #include <math.h>
1.4     ! jtc        27: #include <machine/endian.h>
        !            28:
        !            29: #if BYTE_ORDER == LITTLE_ENDIAN
        !            30: #define n0     1
        !            31: #define n1     0
        !            32: #else
        !            33: #define n0     0
        !            34: #define n1     1
        !            35: #endif
1.1       jtc        36:
                     37: #ifdef __STDC__
1.4     ! jtc        38: static const double zero = 0.0;
1.1       jtc        39: #else
1.4     ! jtc        40: static double zero = 0.0;
1.1       jtc        41: #endif
                     42:
                     43:
                     44: #ifdef __STDC__
                     45:        double __ieee754_remainder(double x, double p)
                     46: #else
                     47:        double __ieee754_remainder(x,p)
                     48:        double x,p;
                     49: #endif
                     50: {
1.4     ! jtc        51:        int hx,hp;
1.1       jtc        52:        unsigned sx,lx,lp;
                     53:        double p_half;
                     54:
                     55:        hx = *( n0 + (int*)&x);         /* high word of x */
                     56:        lx = *( n1 + (int*)&x);         /* low  word of x */
                     57:        hp = *( n0 + (int*)&p);         /* high word of p */
                     58:        lp = *( n1 + (int*)&p);         /* low  word of p */
                     59:        sx = hx&0x80000000;
                     60:        hp &= 0x7fffffff;
                     61:        hx &= 0x7fffffff;
                     62:
                     63:     /* purge off exception values */
                     64:        if((hp|lp)==0) return (x*p)/(x*p);      /* p = 0 */
                     65:        if((hx>=0x7ff00000)||                   /* x not finite */
                     66:          ((hp>=0x7ff00000)&&                   /* p is NaN */
                     67:          (((hp-0x7ff00000)|lp)!=0)))
                     68:            return (x*p)/(x*p);
                     69:
                     70:
                     71:        if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);  /* now x < 2p */
                     72:        if (((hx-hp)|(lx-lp))==0) return zero*x;
                     73:        x  = fabs(x);
                     74:        p  = fabs(p);
                     75:        if (hp<0x00200000) {
                     76:            if(x+x>p) {
                     77:                x-=p;
                     78:                if(x+x>=p) x -= p;
                     79:            }
                     80:        } else {
                     81:            p_half = 0.5*p;
                     82:            if(x>p_half) {
                     83:                x-=p;
                     84:                if(x>=p_half) x -= p;
                     85:            }
                     86:        }
                     87:        *(n0+(int*)&x) ^= sx;
                     88:        return x;
                     89: }