sysdeps/libm-ieee754/e_hypotf.c

   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 #if defined(LIBM_SCCS) && !defined(lint)
  17 static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $";
  18 #endif
  19
  20 #include "math.h"
  21 #include "math_private.h"
  22
  23 #ifdef __STDC__
  24         float __ieee754_hypotf(float x, float y)
  25 #else
  26         float __ieee754_hypot(x,y)
  27         float x, y;
  28 #endif
  29 {
  30         float a=x,b=y,t1,t2,y1,y2,w;
  31         int32_t j,k,ha,hb;
  32
  33         GET_FLOAT_WORD(ha,x);
  34         ha &= 0x7fffffff;
  35         GET_FLOAT_WORD(hb,y);
  36         hb &= 0x7fffffff;
  37         if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
  38         SET_FLOAT_WORD(a,ha);   /* a <- |a| */
  39         SET_FLOAT_WORD(b,hb);   /* b <- |b| */
  40         if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */
  41         k=0;
  42         if(ha > 0x58800000) {   /* a>2**50 */
  43            if(ha >= 0x7f800000) {       /* Inf or NaN */
  44                w = a+b;                 /* for sNaN */
  45                if(ha == 0x7f800000) w = a;
  46                if(hb == 0x7f800000) w = b;
  47                return w;
  48            }
  49            /* scale a and b by 2**-60 */
  50            ha -= 0x5d800000; hb -= 0x5d800000;  k += 60;
  51            SET_FLOAT_WORD(a,ha);
  52            SET_FLOAT_WORD(b,hb);
  53         }
  54         if(hb < 0x26800000) {   /* b < 2**-50 */
  55             if(hb <= 0x007fffff) {      /* subnormal b or 0 */
  56                 if(hb==0) return a;
  57                 SET_FLOAT_WORD(t1,0x3f000000);  /* t1=2^126 */
  58                 b *= t1;
  59                 a *= t1;
  60                 k -= 126;
  61             } else {            /* scale a and b by 2^60 */
  62                 ha += 0x5d800000;       /* a *= 2^60 */
  63                 hb += 0x5d800000;       /* b *= 2^60 */
  64                 k -= 60;
  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  = __ieee754_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);
  80             t2 = a - t1;
  81             w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
  82         }
  83         if(k!=0) {
  84             SET_FLOAT_WORD(t1,0x3f800000+(k<<23));
  85             return t1*w;
  86         } else return w;
  87 }