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 52typedef uint64_t man_t; 53#else 54typedef uint32_t man_t; 55#endif 56 57long double 58hypotl(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 137long double 138hypotl(long double x, long double y) 139{ 140 return hypot(x, y); 141} 142#endif 143