source: trunk/sys/libm/e_remainder.c @ 74

Last change on this file since 74 was 1, checked in by alain, 8 years ago

First import

File size: 1.9 KB
Line 
1
2/* @(#)e_remainder.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/* __ieee754_remainder(x,p)
15 * Return :                 
16 *      returns  x REM p  =  x - [x/p]*p as if in infinite
17 *      precise arithmetic, where [x/p] is the (infinite bit)
18 *      integer nearest x/p (in half way case choose the even one).
19 * Method :
20 *      Based on fmod() return x-[x/p]chopped*p exactlp.
21 */
22
23#include <libm/fdlibm.h>
24
25#ifdef __STDC__
26static const double zero = 0.0, one  = 1.0;
27#else
28static double zero = 0.0, one  = 1.0;
29#endif
30
31
32#ifdef __STDC__
33        double __ieee754_remainder(double x, double p)
34#else
35        double __ieee754_remainder(x,p)
36        double x,p;
37#endif
38{
39        int hx,hp,n0,n1;
40        unsigned sx,lx,lp;
41        double p_half;
42
43        n0 = ((*(int*)&one)>>29)^1;     /* index of high word */
44        n1 = 1-n0;                      /* index of low word */
45        hx = *( n0 + (int*)&x);         /* high word of x */
46        lx = *( n1 + (int*)&x);         /* low  word of x */
47        hp = *( n0 + (int*)&p);         /* high word of p */
48        lp = *( n1 + (int*)&p);         /* low  word of p */
49        sx = hx&0x80000000;
50        hp &= 0x7fffffff;
51        hx &= 0x7fffffff;
52
53    /* purge off exception values */
54        if((hp|lp)==0) return (x*p)/(x*p);      /* p = 0 */
55        if((hx>=0x7ff00000)||                   /* x not finite */
56          ((hp>=0x7ff00000)&&                   /* p is NaN */
57          (((hp-0x7ff00000)|lp)!=0)))
58            return (x*p)/(x*p);
59
60
61        if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);  /* now x < 2p */
62        if (((hx-hp)|(lx-lp))==0) return zero*x;
63        x  = fabs(x);
64        p  = fabs(p);
65        if (hp<0x00200000) {
66            if(x+x>p) {
67                x-=p;
68                if(x+x>=p) x -= p;
69            }
70        } else {
71            p_half = 0.5*p;
72            if(x>p_half) {
73                x-=p;
74                if(x>=p_half) x -= p;
75            }
76        }
77        *(n0+(int*)&x) ^= sx;
78        return x;
79}
Note: See TracBrowser for help on using the repository browser.