source:
EcnlProtoTool/trunk/musl-1.1.18/src/math/remquo.c@
444
Last change on this file since 444 was 444, checked in by , 4 years ago | |
---|---|
|
|
File size: 1.4 KB |
Line | |
---|---|
1 | #include <math.h> |
2 | #include <stdint.h> |
3 | |
4 | double remquo(double x, double y, int *quo) |
5 | { |
6 | union {double f; uint64_t i;} ux = {x}, uy = {y}; |
7 | int ex = ux.i>>52 & 0x7ff; |
8 | int ey = uy.i>>52 & 0x7ff; |
9 | int sx = ux.i>>63; |
10 | int sy = uy.i>>63; |
11 | uint32_t q; |
12 | uint64_t i; |
13 | uint64_t uxi = ux.i; |
14 | |
15 | *quo = 0; |
16 | if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) |
17 | return (x*y)/(x*y); |
18 | if (ux.i<<1 == 0) |
19 | return x; |
20 | |
21 | /* normalize x and y */ |
22 | if (!ex) { |
23 | for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1); |
24 | uxi <<= -ex + 1; |
25 | } else { |
26 | uxi &= -1ULL >> 12; |
27 | uxi |= 1ULL << 52; |
28 | } |
29 | if (!ey) { |
30 | for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1); |
31 | uy.i <<= -ey + 1; |
32 | } else { |
33 | uy.i &= -1ULL >> 12; |
34 | uy.i |= 1ULL << 52; |
35 | } |
36 | |
37 | q = 0; |
38 | if (ex < ey) { |
39 | if (ex+1 == ey) |
40 | goto end; |
41 | return x; |
42 | } |
43 | |
44 | /* x mod y */ |
45 | for (; ex > ey; ex--) { |
46 | i = uxi - uy.i; |
47 | if (i >> 63 == 0) { |
48 | uxi = i; |
49 | q++; |
50 | } |
51 | uxi <<= 1; |
52 | q <<= 1; |
53 | } |
54 | i = uxi - uy.i; |
55 | if (i >> 63 == 0) { |
56 | uxi = i; |
57 | q++; |
58 | } |
59 | if (uxi == 0) |
60 | ex = -60; |
61 | else |
62 | for (; uxi>>52 == 0; uxi <<= 1, ex--); |
63 | end: |
64 | /* scale result and decide between |x| and |x|-|y| */ |
65 | if (ex > 0) { |
66 | uxi -= 1ULL << 52; |
67 | uxi |= (uint64_t)ex << 52; |
68 | } else { |
69 | uxi >>= -ex + 1; |
70 | } |
71 | ux.i = uxi; |
72 | x = ux.f; |
73 | if (sy) |
74 | y = -y; |
75 | if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) { |
76 | x -= y; |
77 | q++; |
78 | } |
79 | q &= 0x7fffffff; |
80 | *quo = sx^sy ? -(int)q : (int)q; |
81 | return sx ? -x : x; |
82 | } |
Note:
See TracBrowser
for help on using the repository browser.