Last change
on this file was 352, checked in by coas-nagasima, 6 years ago |
arm向け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
|
Rev | Line | |
---|
[352] | 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.