Last change
on this file since 351 was 337, checked in by coas-nagasima, 6 years ago |
ASP3版ECNLを追加
|
-
Property svn:eol-style
set to
native
-
Property svn:mime-type
set to
text/x-csrc;charset=UTF-8
|
File size:
1.3 KB
|
Line | |
---|
1 | #include <math.h>
|
---|
2 | #include <stdint.h>
|
---|
3 | #include <float.h>
|
---|
4 |
|
---|
5 | #if FLT_EVAL_METHOD > 1U && LDBL_MANT_DIG == 64
|
---|
6 | #define SPLIT (0x1p32 + 1)
|
---|
7 | #else
|
---|
8 | #define SPLIT (0x1p27 + 1)
|
---|
9 | #endif
|
---|
10 |
|
---|
11 | static void sq(double_t *hi, double_t *lo, double x)
|
---|
12 | {
|
---|
13 | double_t xh, xl, xc;
|
---|
14 |
|
---|
15 | xc = (double_t)x*SPLIT;
|
---|
16 | xh = x - xc + xc;
|
---|
17 | xl = x - xh;
|
---|
18 | *hi = (double_t)x*x;
|
---|
19 | *lo = xh*xh - *hi + 2*xh*xl + xl*xl;
|
---|
20 | }
|
---|
21 |
|
---|
22 | double hypot(double x, double y)
|
---|
23 | {
|
---|
24 | union {double f; uint64_t i;} ux = {x}, uy = {y}, ut;
|
---|
25 | int ex, ey;
|
---|
26 | double_t hx, lx, hy, ly, z;
|
---|
27 |
|
---|
28 | /* arrange |x| >= |y| */
|
---|
29 | ux.i &= -1ULL>>1;
|
---|
30 | uy.i &= -1ULL>>1;
|
---|
31 | if (ux.i < uy.i) {
|
---|
32 | ut = ux;
|
---|
33 | ux = uy;
|
---|
34 | uy = ut;
|
---|
35 | }
|
---|
36 |
|
---|
37 | /* special cases */
|
---|
38 | ex = ux.i>>52;
|
---|
39 | ey = uy.i>>52;
|
---|
40 | x = ux.f;
|
---|
41 | y = uy.f;
|
---|
42 | /* note: hypot(inf,nan) == inf */
|
---|
43 | if (ey == 0x7ff)
|
---|
44 | return y;
|
---|
45 | if (ex == 0x7ff || uy.i == 0)
|
---|
46 | return x;
|
---|
47 | /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */
|
---|
48 | /* 64 difference is enough for ld80 double_t */
|
---|
49 | if (ex - ey > 64)
|
---|
50 | return x + y;
|
---|
51 |
|
---|
52 | /* precise sqrt argument in nearest rounding mode without overflow */
|
---|
53 | /* xh*xh must not overflow and xl*xl must not underflow in sq */
|
---|
54 | z = 1;
|
---|
55 | if (ex > 0x3ff+510) {
|
---|
56 | z = 0x1p700;
|
---|
57 | x *= 0x1p-700;
|
---|
58 | y *= 0x1p-700;
|
---|
59 | } else if (ey < 0x3ff-450) {
|
---|
60 | z = 0x1p-700;
|
---|
61 | x *= 0x1p700;
|
---|
62 | y *= 0x1p700;
|
---|
63 | }
|
---|
64 | sq(&hx, &lx, x);
|
---|
65 | sq(&hy, &ly, y);
|
---|
66 | return z*sqrt(ly+lx+hy+hx);
|
---|
67 | }
|
---|
Note:
See
TracBrowser
for help on using the repository browser.