Go to most recent revision | Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
3362 | Serge | 1 | |
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 | */ |
||
12 | |||
13 | |||
14 | * Return : |
||
15 | * returns x REM p = x - [x/p]*p as if in infinite |
||
16 | * precise arithmetic, where [x/p] is the (infinite bit) |
||
17 | * integer nearest x/p (in half way case choose the even one). |
||
18 | * Method : |
||
19 | * Based on fmod() return x-[x/p]chopped*p exactlp. |
||
20 | */ |
||
21 | |||
22 | |||
23 | |||
24 | |||
25 | |||
26 | |||
27 | static const double zero = 0.0; |
||
28 | #else |
||
29 | static double zero = 0.0; |
||
30 | #endif |
||
31 | |||
32 | |||
33 | |||
34 | double __ieee754_remainder(double x, double p) |
||
35 | #else |
||
36 | double __ieee754_remainder(x,p) |
||
37 | double x,p; |
||
38 | #endif |
||
39 | { |
||
40 | __int32_t hx,hp; |
||
41 | __uint32_t sx,lx,lp; |
||
42 | double p_half; |
||
43 | |||
44 | |||
45 | EXTRACT_WORDS(hp,lp,p); |
||
46 | sx = hx&0x80000000; |
||
47 | hp &= 0x7fffffff; |
||
48 | hx &= 0x7fffffff; |
||
49 | |||
50 | |||
51 | if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */ |
||
52 | if((hx>=0x7ff00000)|| /* x not finite */ |
||
53 | ((hp>=0x7ff00000)&& /* p is NaN */ |
||
54 | (((hp-0x7ff00000)|lp)!=0))) |
||
55 | return (x*p)/(x*p); |
||
56 | |||
57 | |||
58 | |||
59 | if (((hx-hp)|(lx-lp))==0) return zero*x; |
||
60 | x = fabs(x); |
||
61 | p = fabs(p); |
||
62 | if (hp<0x00200000) { |
||
63 | if(x+x>p) { |
||
64 | x-=p; |
||
65 | if(x+x>=p) x -= p; |
||
66 | } |
||
67 | } else { |
||
68 | p_half = 0.5*p; |
||
69 | if(x>p_half) { |
||
70 | x-=p; |
||
71 | if(x>=p_half) x -= p; |
||
72 | } |
||
73 | } |
||
74 | GET_HIGH_WORD(hx,x); |
||
75 | SET_HIGH_WORD(x,hx^sx); |
||
76 | return x; |
||
77 | } |
||
78 | |||
79 | |||
80 |