src/math/hypotf.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_hypotf.c */
   2 /*
   3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   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 #include "libm.h"
  17
  18 float hypotf(float x, float y)
  19 {
  20         float a,b,t1,t2,y1,y2,w;
  21         int32_t j,k,ha,hb;
  22
  23         GET_FLOAT_WORD(ha,x);
  24         ha &= 0x7fffffff;
  25         GET_FLOAT_WORD(hb,y);
  26         hb &= 0x7fffffff;
  27         if (hb > ha) {
  28                 a = y;
  29                 b = x;
  30                 j=ha; ha=hb; hb=j;
  31         } else {
  32                 a = x;
  33                 b = y;
  34         }
  35         a = fabsf(a);
  36         b = fabsf(b);
  37         if (ha - hb > 0xf000000)  /* x/y > 2**30 */
  38                 return a+b;
  39         k = 0;
  40         if (ha > 0x58800000) {    /* a > 2**50 */
  41                 if(ha >= 0x7f800000) {  /* Inf or NaN */
  42                         /* Use original arg order iff result is NaN; quieten sNaNs. */
  43                         w = fabsf(x+0.0f) - fabsf(y+0.0f);
  44                         if (ha == 0x7f800000) w = a;
  45                         if (hb == 0x7f800000) w = b;
  46                         return w;
  47                 }
  48                 /* scale a and b by 2**-68 */
  49                 ha -= 0x22000000; hb -= 0x22000000; k += 68;
  50                 SET_FLOAT_WORD(a, ha);
  51                 SET_FLOAT_WORD(b, hb);
  52         }
  53         if (hb < 0x26800000) {    /* b < 2**-50 */
  54                 if (hb <= 0x007fffff) {  /* subnormal b or 0 */
  55                         if (hb == 0)
  56                                 return a;
  57                         SET_FLOAT_WORD(t1, 0x7e800000);  /* t1 = 2^126 */
  58                         b *= t1;
  59                         a *= t1;
  60                         k -= 126;
  61                 } else {   /* scale a and b by 2^68 */
  62                         ha += 0x22000000;  /* a *= 2^68 */
  63                         hb += 0x22000000;  /* b *= 2^68 */
  64                         k -= 68;
  65                         SET_FLOAT_WORD(a, ha);
  66                         SET_FLOAT_WORD(b, hb);
  67                 }
  68         }
  69         /* medium size a and b */
  70         w = a - b;
  71         if (w > b) {
  72                 SET_FLOAT_WORD(t1, ha&0xfffff000);
  73                 t2 = a - t1;
  74                 w  = sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
  75         } else {
  76                 a  = a + a;
  77                 SET_FLOAT_WORD(y1, hb&0xfffff000);
  78                 y2 = b - y1;
  79                 SET_FLOAT_WORD(t1,(ha+0x00800000)&0xfffff000);
  80                 t2 = a - t1;
  81                 w  = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
  82         }
  83         if (k)
  84                 w = scalbnf(w, k);
  85         return w;
  86 }