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.2 KB
|
Line | |
---|
1 | #include "libm.h"
|
---|
2 |
|
---|
3 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
---|
4 | long double hypotl(long double x, long double y)
|
---|
5 | {
|
---|
6 | return hypot(x, y);
|
---|
7 | }
|
---|
8 | #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
|
---|
9 | #if LDBL_MANT_DIG == 64
|
---|
10 | #define SPLIT (0x1p32L+1)
|
---|
11 | #elif LDBL_MANT_DIG == 113
|
---|
12 | #define SPLIT (0x1p57L+1)
|
---|
13 | #endif
|
---|
14 |
|
---|
15 | static void sq(long double *hi, long double *lo, long double x)
|
---|
16 | {
|
---|
17 | long double xh, xl, xc;
|
---|
18 | xc = x*SPLIT;
|
---|
19 | xh = x - xc + xc;
|
---|
20 | xl = x - xh;
|
---|
21 | *hi = x*x;
|
---|
22 | *lo = xh*xh - *hi + 2*xh*xl + xl*xl;
|
---|
23 | }
|
---|
24 |
|
---|
25 | long double hypotl(long double x, long double y)
|
---|
26 | {
|
---|
27 | union ldshape ux = {x}, uy = {y};
|
---|
28 | int ex, ey;
|
---|
29 | long double hx, lx, hy, ly, z;
|
---|
30 |
|
---|
31 | ux.i.se &= 0x7fff;
|
---|
32 | uy.i.se &= 0x7fff;
|
---|
33 | if (ux.i.se < uy.i.se) {
|
---|
34 | ex = uy.i.se;
|
---|
35 | ey = ux.i.se;
|
---|
36 | x = uy.f;
|
---|
37 | y = ux.f;
|
---|
38 | } else {
|
---|
39 | ex = ux.i.se;
|
---|
40 | ey = uy.i.se;
|
---|
41 | x = ux.f;
|
---|
42 | y = uy.f;
|
---|
43 | }
|
---|
44 |
|
---|
45 | if (ex == 0x7fff && isinf(y))
|
---|
46 | return y;
|
---|
47 | if (ex == 0x7fff || y == 0)
|
---|
48 | return x;
|
---|
49 | if (ex - ey > LDBL_MANT_DIG)
|
---|
50 | return x + y;
|
---|
51 |
|
---|
52 | z = 1;
|
---|
53 | if (ex > 0x3fff+8000) {
|
---|
54 | z = 0x1p10000L;
|
---|
55 | x *= 0x1p-10000L;
|
---|
56 | y *= 0x1p-10000L;
|
---|
57 | } else if (ey < 0x3fff-8000) {
|
---|
58 | z = 0x1p-10000L;
|
---|
59 | x *= 0x1p10000L;
|
---|
60 | y *= 0x1p10000L;
|
---|
61 | }
|
---|
62 | sq(&hx, &lx, x);
|
---|
63 | sq(&hy, &ly, y);
|
---|
64 | return z*sqrtl(ly+lx+hy+hx);
|
---|
65 | }
|
---|
66 | #endif
|
---|
Note:
See
TracBrowser
for help on using the repository browser.