| The Android Open Source Project | 1dc9e47 | 2009-03-03 19:28:35 -0800 | [diff] [blame] | 1 | /* e_hypotf.c -- float version of e_hypot.c. | 
|  | 2 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. | 
|  | 3 | */ | 
|  | 4 |  | 
|  | 5 | /* | 
|  | 6 | * ==================================================== | 
|  | 7 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 
|  | 8 | * | 
|  | 9 | * Developed at SunPro, a Sun Microsystems, Inc. business. | 
|  | 10 | * Permission to use, copy, modify, and distribute this | 
|  | 11 | * software is freely granted, provided that this notice | 
|  | 12 | * is preserved. | 
|  | 13 | * ==================================================== | 
|  | 14 | */ | 
|  | 15 |  | 
|  | 16 | #ifndef lint | 
|  | 17 | static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_hypotf.c,v 1.9 2002/05/28 18:15:03 alfred Exp $"; | 
|  | 18 | #endif | 
|  | 19 |  | 
|  | 20 | #include "math.h" | 
|  | 21 | #include "math_private.h" | 
|  | 22 |  | 
|  | 23 | float | 
|  | 24 | __ieee754_hypotf(float x, float y) | 
|  | 25 | { | 
|  | 26 | float a=x,b=y,t1,t2,y1,y2,w; | 
|  | 27 | int32_t j,k,ha,hb; | 
|  | 28 |  | 
|  | 29 | GET_FLOAT_WORD(ha,x); | 
|  | 30 | ha &= 0x7fffffff; | 
|  | 31 | GET_FLOAT_WORD(hb,y); | 
|  | 32 | hb &= 0x7fffffff; | 
|  | 33 | if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} | 
|  | 34 | SET_FLOAT_WORD(a,ha);	/* a <- |a| */ | 
|  | 35 | SET_FLOAT_WORD(b,hb);	/* b <- |b| */ | 
|  | 36 | if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */ | 
|  | 37 | k=0; | 
|  | 38 | if(ha > 0x58800000) {	/* a>2**50 */ | 
|  | 39 | if(ha >= 0x7f800000) {	/* Inf or NaN */ | 
|  | 40 | w = a+b;			/* for sNaN */ | 
|  | 41 | if(ha == 0x7f800000) w = a; | 
|  | 42 | if(hb == 0x7f800000) w = b; | 
|  | 43 | return w; | 
|  | 44 | } | 
|  | 45 | /* scale a and b by 2**-68 */ | 
|  | 46 | ha -= 0x22000000; hb -= 0x22000000;	k += 68; | 
|  | 47 | SET_FLOAT_WORD(a,ha); | 
|  | 48 | SET_FLOAT_WORD(b,hb); | 
|  | 49 | } | 
|  | 50 | if(hb < 0x26800000) {	/* b < 2**-50 */ | 
|  | 51 | if(hb <= 0x007fffff) {	/* subnormal b or 0 */ | 
|  | 52 | if(hb==0) return a; | 
|  | 53 | SET_FLOAT_WORD(t1,0x7e800000);	/* t1=2^126 */ | 
|  | 54 | b *= t1; | 
|  | 55 | a *= t1; | 
|  | 56 | k -= 126; | 
|  | 57 | } else {		/* scale a and b by 2^68 */ | 
|  | 58 | ha += 0x22000000; 	/* a *= 2^68 */ | 
|  | 59 | hb += 0x22000000;	/* b *= 2^68 */ | 
|  | 60 | k -= 68; | 
|  | 61 | SET_FLOAT_WORD(a,ha); | 
|  | 62 | SET_FLOAT_WORD(b,hb); | 
|  | 63 | } | 
|  | 64 | } | 
|  | 65 | /* medium size a and b */ | 
|  | 66 | w = a-b; | 
|  | 67 | if (w>b) { | 
|  | 68 | SET_FLOAT_WORD(t1,ha&0xfffff000); | 
|  | 69 | t2 = a-t1; | 
|  | 70 | w  = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1))); | 
|  | 71 | } else { | 
|  | 72 | a  = a+a; | 
|  | 73 | SET_FLOAT_WORD(y1,hb&0xfffff000); | 
|  | 74 | y2 = b - y1; | 
|  | 75 | SET_FLOAT_WORD(t1,ha+0x00800000); | 
|  | 76 | t2 = a - t1; | 
|  | 77 | w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b))); | 
|  | 78 | } | 
|  | 79 | if(k!=0) { | 
|  | 80 | SET_FLOAT_WORD(t1,0x3f800000+(k<<23)); | 
|  | 81 | return t1*w; | 
|  | 82 | } else return w; | 
|  | 83 | } |