newlib/libm/mathfp/e_hypot.c

   1
   2 /* @(#)e_hypot.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14 /*
  15 FUNCTION
  16         <<hypot>>, <<hypotf>>---distance from origin
  17 INDEX
  18         hypot
  19 INDEX
  20         hypotf
  21
  22 SYNOPSIS
  23         #include <math.h>
  24         double hypot(double <[x]>, double <[y]>);
  25         float hypotf(float <[x]>, float <[y]>);
  26
  27 DESCRIPTION
  28         <<hypot>> calculates the Euclidean distance
  29         @tex
  30         $\sqrt{x^2+y^2}$
  31         @end tex
  32         @ifnottex
  33         <<sqrt(<[x]>*<[x]> + <[y]>*<[y]>)>>
  34         @end ifnottex
  35         between the origin (0,0) and a point represented by the
  36         Cartesian coordinates (<[x]>,<[y]>).  <<hypotf>> differs only
  37         in the type of its arguments and result.
  38
  39 RETURNS
  40         Normally, the distance value is returned.  On overflow,
  41         <<hypot>> returns <<HUGE_VAL>> and sets <<errno>> to
  42         <<ERANGE>>.
  43
  44 PORTABILITY
  45         <<hypot>> and <<hypotf>> are not ANSI C.  */
  46
  47 /* hypot(x,y)
  48  *
  49  * Method :
  50  *      If (assume round-to-nearest) z=x*x+y*y
  51  *      has error less than sqrt(2)/2 ulp, than
  52  *      sqrt(z) has error less than 1 ulp (exercise).
  53  *
  54  *      So, compute sqrt(x*x+y*y) with some care as
  55  *      follows to get the error below 1 ulp:
  56  *
  57  *      Assume x>y>0;
  58  *      (if possible, set rounding to round-to-nearest)
  59  *      1. if x > 2y  use
  60  *              x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
  61  *      where x1 = x with lower 32 bits cleared, x2 = x-x1; else
  62  *      2. if x <= 2y use
  63  *              t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
  64  *      where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
  65  *      y1= y with lower 32 bits chopped, y2 = y-y1.
  66  *
  67  *      NOTE: scaling may be necessary if some argument is too
  68  *            large or too tiny
  69  *
  70  * Special cases:
  71  *      hypot(x,y) is INF if x or y is +INF or -INF; else
  72  *      hypot(x,y) is NAN if x or y is NAN.
  73  *
  74  * Accuracy:
  75  *      hypot(x,y) returns sqrt(x^2+y^2) with error less
  76  *      than 1 ulps (units in the last place)
  77  */
  78
  79 #include "fdlibm.h"
  80
  81 #ifndef _DOUBLE_IS_32BITS
  82
  83 #ifdef __STDC__
  84         double hypot(double x, double y)
  85 #else
  86         double hypot(x,y)
  87         double x, y;
  88 #endif
  89 {
  90         double a=x,b=y,t1,t2,y1,y2,w;
  91         __int32_t j,k,ha,hb;
  92
  93         GET_HIGH_WORD(ha,x);
  94         ha &= 0x7fffffff;
  95         GET_HIGH_WORD(hb,y);
  96         hb &= 0x7fffffff;
  97         if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
  98         SET_HIGH_WORD(a,ha);    /* a <- |a| */
  99         SET_HIGH_WORD(b,hb);    /* b <- |b| */
 100         if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
 101         k=0;
 102         if(ha > 0x5f300000) {   /* a>2**500 */
 103            if(ha >= 0x7ff00000) {       /* Inf or NaN */
 104                __uint32_t low;
 105                w = a+b;                 /* for sNaN */
 106                GET_LOW_WORD(low,a);
 107                if(((ha&0xfffff)|low)==0) w = a;
 108                GET_LOW_WORD(low,b);
 109                if(((hb^0x7ff00000)|low)==0) w = b;
 110                return w;
 111            }
 112            /* scale a and b by 2**-600 */
 113            ha -= 0x25800000; hb -= 0x25800000;  k += 600;
 114            SET_HIGH_WORD(a,ha);
 115            SET_HIGH_WORD(b,hb);
 116         }
 117         if(hb < 0x20b00000) {   /* b < 2**-500 */
 118             if(hb <= 0x000fffff) {      /* subnormal b or 0 */
 119                 __uint32_t low;
 120                 GET_LOW_WORD(low,b);
 121                 if((hb|low)==0) return a;
 122                 t1=0;
 123                 SET_HIGH_WORD(t1,0x7fd00000);   /* t1=2^1022 */
 124                 b *= t1;
 125                 a *= t1;
 126                 k -= 1022;
 127             } else {            /* scale a and b by 2^600 */
 128                 ha += 0x25800000;       /* a *= 2^600 */
 129                 hb += 0x25800000;       /* b *= 2^600 */
 130                 k -= 600;
 131                 SET_HIGH_WORD(a,ha);
 132                 SET_HIGH_WORD(b,hb);
 133             }
 134         }
 135     /* medium size a and b */
 136         w = a-b;
 137         if (w>b) {
 138             t1 = 0;
 139             SET_HIGH_WORD(t1,ha);
 140             t2 = a-t1;
 141             w  = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
 142         } else {
 143             a  = a+a;
 144             y1 = 0;
 145             SET_HIGH_WORD(y1,hb);
 146             y2 = b - y1;
 147             t1 = 0;
 148             SET_HIGH_WORD(t1,ha+0x00100000);
 149             t2 = a - t1;
 150             w  = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 151         }
 152         if(k!=0) {
 153             __uint32_t high;
 154             t1 = 1.0;
 155             GET_HIGH_WORD(high,t1);
 156             SET_HIGH_WORD(t1,high+(k<<20));
 157             return t1*w;
 158         } else return w;
 159 }
 160
 161 #endif /* defined(_DOUBLE_IS_32BITS) */