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