123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122 |
- #include <openlibm_math.h>
- #include "math_private.h"
- long double
- hypotl(long double x, long double y)
- {
- long double a,b,t1,t2,yy1,y2,w;
- int64_t j,k,ha,hb;
- GET_LDOUBLE_MSW64(ha,x);
- ha &= 0x7fffffffffffffffLL;
- GET_LDOUBLE_MSW64(hb,y);
- hb &= 0x7fffffffffffffffLL;
- if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
- SET_LDOUBLE_MSW64(a,ha);
- SET_LDOUBLE_MSW64(b,hb);
- if((ha-hb)>0x78000000000000LL) {return a+b;}
- k=0;
- if(ha > 0x5f3f000000000000LL) {
- if(ha >= 0x7fff000000000000LL) {
- u_int64_t low;
- w = a+b;
- GET_LDOUBLE_LSW64(low,a);
- if(((ha&0xffffffffffffLL)|low)==0) w = a;
- GET_LDOUBLE_LSW64(low,b);
- if(((hb^0x7fff000000000000LL)|low)==0) w = b;
- return w;
- }
-
- ha -= 0x2580000000000000LL;
- hb -= 0x2580000000000000LL; k += 9600;
- SET_LDOUBLE_MSW64(a,ha);
- SET_LDOUBLE_MSW64(b,hb);
- }
- if(hb < 0x20bf000000000000LL) {
- if(hb <= 0x0000ffffffffffffLL) {
- u_int64_t low;
- GET_LDOUBLE_LSW64(low,b);
- if((hb|low)==0) return a;
- t1=0;
- SET_LDOUBLE_MSW64(t1,0x7ffd000000000000LL);
- b *= t1;
- a *= t1;
- k -= 16382;
- } else {
- ha += 0x2580000000000000LL;
- hb += 0x2580000000000000LL;
- k -= 9600;
- SET_LDOUBLE_MSW64(a,ha);
- SET_LDOUBLE_MSW64(b,hb);
- }
- }
-
- w = a-b;
- if (w>b) {
- t1 = 0;
- SET_LDOUBLE_MSW64(t1,ha);
- t2 = a-t1;
- w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
- } else {
- a = a+a;
- yy1 = 0;
- SET_LDOUBLE_MSW64(yy1,hb);
- y2 = b - yy1;
- t1 = 0;
- SET_LDOUBLE_MSW64(t1,ha+0x0001000000000000LL);
- t2 = a - t1;
- w = sqrtl(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
- }
- if(k!=0) {
- u_int64_t high;
- t1 = 1.0L;
- GET_LDOUBLE_MSW64(high,t1);
- SET_LDOUBLE_MSW64(t1,high+(k<<48));
- return t1*w;
- } else return w;
- }
|