00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
#if defined(LIBM_SCCS) && !defined(lint)
00014
static char rcsid[] =
"$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
00015
#endif
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
#include "math.h"
00050
#include "mathP.h"
00051
00052
#ifdef __STDC__
00053
double __ieee754_hypot(
double x,
double y)
00054 #
else
00055 double __ieee754_hypot(x,y)
00056 double x, y;
00057 #endif
00058 {
00059
double a=x,b=y,t1,t2,
y1,y2,w;
00060 int32_t j,k,ha,hb;
00061
00062
GET_HIGH_WORD(ha,x);
00063 ha &= 0x7fffffff;
00064
GET_HIGH_WORD(hb,y);
00065 hb &= 0x7fffffff;
00066
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;}
else {a=x;b=y;}
00067
SET_HIGH_WORD(a,ha);
00068
SET_HIGH_WORD(b,hb);
00069
if((ha-hb)>0x3c00000) {
return a+b;}
00070 k=0;
00071
if(ha > 0x5f300000) {
00072
if(ha >= 0x7ff00000) {
00073 u_int32_t low;
00074 w = a+b;
00075
GET_LOW_WORD(low,a);
00076
if(((ha&0xfffff)|low)==0) w = a;
00077
GET_LOW_WORD(low,b);
00078
if(((hb^0x7ff00000)|low)==0) w = b;
00079
return w;
00080 }
00081
00082 ha -= 0x25800000; hb -= 0x25800000; k += 600;
00083
SET_HIGH_WORD(a,ha);
00084
SET_HIGH_WORD(b,hb);
00085 }
00086
if(hb < 0x20b00000) {
00087
if(hb <= 0x000fffff) {
00088 u_int32_t low;
00089
GET_LOW_WORD(low,b);
00090
if((hb|low)==0)
return a;
00091 t1=0;
00092
SET_HIGH_WORD(t1,0x7fd00000);
00093 b *= t1;
00094 a *= t1;
00095 k -= 1022;
00096 }
else {
00097 ha += 0x25800000;
00098 hb += 0x25800000;
00099 k -= 600;
00100
SET_HIGH_WORD(a,ha);
00101
SET_HIGH_WORD(b,hb);
00102 }
00103 }
00104
00105 w = a-b;
00106
if (w>b) {
00107 t1 = 0;
00108
SET_HIGH_WORD(t1,ha);
00109 t2 = a-t1;
00110 w =
__ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
00111 }
else {
00112 a = a+a;
00113
y1 = 0;
00114
SET_HIGH_WORD(
y1,hb);
00115 y2 = b -
y1;
00116 t1 = 0;
00117
SET_HIGH_WORD(t1,ha+0x00100000);
00118 t2 = a - t1;
00119 w =
__ieee754_sqrt(t1*
y1-(w*(-w)-(t1*y2+t2*b)));
00120 }
00121
if(k!=0) {
00122 u_int32_t high;
00123 t1 = 1.0;
00124
GET_HIGH_WORD(high,t1);
00125
SET_HIGH_WORD(t1,high+(k<<20));
00126
return t1*w;
00127 }
else return w;
00128 }