source:
azure_iot_hub/trunk/musl-1.1.18/src/math/remquo.c
Last change on this file was 389, checked in by , 5 years ago | |
---|---|
|
|
File size: 1.4 KB |
Rev | Line | |
---|---|---|
[388] | 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.