Last change
on this file since 445 was 444, checked in by coas-nagasima, 4 years ago |
muslのソースコードを追加
|
-
Property svn:eol-style
set to
native
-
Property svn:mime-type
set to
text/x-csrc;charset=UTF-8
|
File size:
1.3 KB
|
Rev | Line | |
---|
[444] | 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.