Elliott Hughes | a0ee078 | 2013-01-30 19:06:37 -0800 | [diff] [blame] | 1 | /* From: @(#)e_hypot.c 1.3 95/01/18 */ |
| 2 | /* |
| 3 | * ==================================================== |
| 4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 5 | * |
| 6 | * Developed at SunSoft, a Sun Microsystems, Inc. business. |
| 7 | * Permission to use, copy, modify, and distribute this |
| 8 | * software is freely granted, provided that this notice |
| 9 | * is preserved. |
| 10 | * ==================================================== |
| 11 | */ |
| 12 | |
| 13 | #include <sys/cdefs.h> |
| 14 | __FBSDID("$FreeBSD$"); |
| 15 | |
| 16 | /* long double version of hypot(). See e_hypot.c for most comments. */ |
| 17 | |
| 18 | #include <float.h> |
| 19 | |
| 20 | #include "fpmath.h" |
| 21 | #include "math.h" |
| 22 | #include "math_private.h" |
| 23 | |
| 24 | #define GET_LDBL_MAN(h, l, v) do { \ |
| 25 | union IEEEl2bits uv; \ |
| 26 | \ |
| 27 | uv.e = v; \ |
| 28 | h = uv.bits.manh; \ |
| 29 | l = uv.bits.manl; \ |
| 30 | } while (0) |
| 31 | |
| 32 | #undef GET_HIGH_WORD |
| 33 | #define GET_HIGH_WORD(i, v) GET_LDBL_EXPSIGN(i, v) |
| 34 | #undef SET_HIGH_WORD |
| 35 | #define SET_HIGH_WORD(v, i) SET_LDBL_EXPSIGN(v, i) |
| 36 | |
| 37 | #define DESW(exp) (exp) /* delta expsign word */ |
| 38 | #define ESW(exp) (MAX_EXP - 1 + (exp)) /* expsign word */ |
| 39 | #define MANT_DIG LDBL_MANT_DIG |
| 40 | #define MAX_EXP LDBL_MAX_EXP |
| 41 | |
| 42 | #if LDBL_MANL_SIZE > 32 |
| 43 | typedef uint64_t man_t; |
| 44 | #else |
| 45 | typedef uint32_t man_t; |
| 46 | #endif |
| 47 | |
| 48 | long double |
| 49 | hypotl(long double x, long double y) |
| 50 | { |
| 51 | long double a=x,b=y,t1,t2,y1,y2,w; |
| 52 | int32_t j,k,ha,hb; |
| 53 | |
| 54 | GET_HIGH_WORD(ha,x); |
| 55 | ha &= 0x7fff; |
| 56 | GET_HIGH_WORD(hb,y); |
| 57 | hb &= 0x7fff; |
| 58 | if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} |
| 59 | a = fabsl(a); |
| 60 | b = fabsl(b); |
| 61 | if((ha-hb)>DESW(MANT_DIG+7)) {return a+b;} /* x/y > 2**(MANT_DIG+7) */ |
| 62 | k=0; |
| 63 | if(ha > ESW(MAX_EXP/2-12)) { /* a>2**(MAX_EXP/2-12) */ |
| 64 | if(ha >= ESW(MAX_EXP)) { /* Inf or NaN */ |
| 65 | man_t manh, manl; |
| 66 | /* Use original arg order iff result is NaN; quieten sNaNs. */ |
| 67 | w = fabsl(x+0.0)-fabsl(y+0.0); |
| 68 | GET_LDBL_MAN(manh,manl,a); |
| 69 | if (manh == LDBL_NBIT && manl == 0) w = a; |
| 70 | GET_LDBL_MAN(manh,manl,b); |
| 71 | if (hb >= ESW(MAX_EXP) && manh == LDBL_NBIT && manl == 0) w = b; |
| 72 | return w; |
| 73 | } |
| 74 | /* scale a and b by 2**-(MAX_EXP/2+88) */ |
| 75 | ha -= DESW(MAX_EXP/2+88); hb -= DESW(MAX_EXP/2+88); |
| 76 | k += MAX_EXP/2+88; |
| 77 | SET_HIGH_WORD(a,ha); |
| 78 | SET_HIGH_WORD(b,hb); |
| 79 | } |
| 80 | if(hb < ESW(-(MAX_EXP/2-12))) { /* b < 2**-(MAX_EXP/2-12) */ |
| 81 | if(hb <= 0) { /* subnormal b or 0 */ |
| 82 | man_t manh, manl; |
| 83 | GET_LDBL_MAN(manh,manl,b); |
| 84 | if((manh|manl)==0) return a; |
| 85 | t1=0; |
| 86 | SET_HIGH_WORD(t1,ESW(MAX_EXP-2)); /* t1=2^(MAX_EXP-2) */ |
| 87 | b *= t1; |
| 88 | a *= t1; |
| 89 | k -= MAX_EXP-2; |
| 90 | } else { /* scale a and b by 2^(MAX_EXP/2+88) */ |
| 91 | ha += DESW(MAX_EXP/2+88); |
| 92 | hb += DESW(MAX_EXP/2+88); |
| 93 | k -= MAX_EXP/2+88; |
| 94 | SET_HIGH_WORD(a,ha); |
| 95 | SET_HIGH_WORD(b,hb); |
| 96 | } |
| 97 | } |
| 98 | /* medium size a and b */ |
| 99 | w = a-b; |
| 100 | if (w>b) { |
| 101 | t1 = a; |
| 102 | union IEEEl2bits uv; |
| 103 | uv.e = t1; uv.bits.manl = 0; t1 = uv.e; |
| 104 | t2 = a-t1; |
| 105 | w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1))); |
| 106 | } else { |
| 107 | a = a+a; |
| 108 | y1 = b; |
| 109 | union IEEEl2bits uv; |
| 110 | uv.e = y1; uv.bits.manl = 0; y1 = uv.e; |
| 111 | y2 = b - y1; |
| 112 | t1 = a; |
| 113 | uv.e = t1; uv.bits.manl = 0; t1 = uv.e; |
| 114 | t2 = a - t1; |
| 115 | w = sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b))); |
| 116 | } |
| 117 | if(k!=0) { |
| 118 | u_int32_t high; |
| 119 | t1 = 1.0; |
| 120 | GET_HIGH_WORD(high,t1); |
| 121 | SET_HIGH_WORD(t1,high+DESW(k)); |
| 122 | return t1*w; |
| 123 | } else return w; |
| 124 | } |