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;
- u_int32_t j,k,ea,eb;
- GET_LDOUBLE_EXP(ea,x);
- ea &= 0x7fff;
- GET_LDOUBLE_EXP(eb,y);
- eb &= 0x7fff;
- if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;}
- SET_LDOUBLE_EXP(a,ea);
- SET_LDOUBLE_EXP(b,eb);
- if((ea-eb)>0x46) {return a+b;}
- k=0;
- if(ea > 0x5f3f) {
- if(ea == 0x7fff) {
- u_int32_t es,high,low;
- w = a+b;
- GET_LDOUBLE_WORDS(es,high,low,a);
- if(((high&0x7fffffff)|low)==0) w = a;
- GET_LDOUBLE_WORDS(es,high,low,b);
- if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b;
- return w;
- }
-
- ea -= 0x2580; eb -= 0x2580; k += 9600;
- SET_LDOUBLE_EXP(a,ea);
- SET_LDOUBLE_EXP(b,eb);
- }
- if(eb < 0x20bf) {
- if(eb == 0) {
- u_int32_t es,high,low;
- GET_LDOUBLE_WORDS(es,high,low,b);
- if((high|low)==0) return a;
- SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0);
- b *= t1;
- a *= t1;
- k -= 16382;
- } else {
- ea += 0x2580;
- eb += 0x2580;
- k -= 9600;
- SET_LDOUBLE_EXP(a,ea);
- SET_LDOUBLE_EXP(b,eb);
- }
- }
-
- w = a-b;
- if (w>b) {
- u_int32_t high;
- GET_LDOUBLE_MSW(high,a);
- SET_LDOUBLE_WORDS(t1,ea,high,0);
- t2 = a-t1;
- w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
- } else {
- u_int32_t high;
- GET_LDOUBLE_MSW(high,b);
- a = a+a;
- SET_LDOUBLE_WORDS(yy1,eb,high,0);
- y2 = b - yy1;
- GET_LDOUBLE_MSW(high,a);
- SET_LDOUBLE_WORDS(t1,ea+1,high,0);
- t2 = a - t1;
- w = sqrtl(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
- }
- if(k!=0) {
- u_int32_t es;
- t1 = 1.0;
- GET_LDOUBLE_EXP(es,t1);
- SET_LDOUBLE_EXP(t1,es+k);
- return t1*w;
- } else return w;
- }
|