winsup/cygwin/math/hypotl.c

   1 /**
   2  * This file has no copyright assigned and is placed in the Public Domain.
   3  * This file is part of the mingw-w64 runtime package.
   4  * No warranty is given; refer to the file DISCLAIMER.PD within this package.
   5  */
   6 #include <math.h>
   7 #include <float.h>
   8 #include <errno.h>
   9
  10 /*
  11 This implementation is based largely on Cephes library
  12 function cabsl (cmplxl.c), which bears the following notice:
  13
  14 Cephes Math Library Release 2.1:  January, 1989
  15 Copyright 1984, 1987, 1989 by Stephen L. Moshier
  16 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  17 */
  18
  19 /*
  20    Modified for use in libmingwex.a
  21    02 Sept 2002  Danny Smith  <dannysmith@users.sourceforege.net>
  22    Calls to ldexpl replaced by logbl and calls to frexpl replaced
  23    by scalbnl to avoid duplicated range checks.
  24 */
  25
  26 #define PRECL 32
  27
  28 long double
  29 hypotl (long double x, long double y)
  30 {
  31   int exx;
  32   int eyy;
  33   int  scale;
  34   long double xx =fabsl(x);
  35   long double yy =fabsl(y);
  36   if (!isfinite(xx) || !isfinite(yy))
  37     {
  38       /* Annex F.9.4.3, hypot returns +infinity if
  39          either component is an infinity, even when the
  40          other component is NaN.  */
  41       return (isinf(xx) || isinf(yy)) ? INFINITY : NAN;
  42     }
  43
  44   if (xx == 0.0L)
  45      return yy;
  46   if (yy == 0.0L)
  47      return xx;
  48
  49   /* Get exponents */
  50   exx =  logbl (xx);
  51   eyy =  logbl (yy);
  52
  53   /* Check if large differences in scale */
  54   scale = exx - eyy;
  55   if ( scale > PRECL)
  56      return xx;
  57   if ( scale < -PRECL)
  58      return yy;
  59
  60   /* Exponent of approximate geometric mean (x 2) */
  61   scale = (exx + eyy) >> 1;
  62
  63   /*  Rescale: Geometric mean is now about 2 */
  64   x = scalbnl(xx, -scale);
  65   y = scalbnl(yy, -scale);
  66
  67   xx = sqrtl(x * x  + y * y);
  68
  69   /* Check for overflow and underflow */
  70   exx = logbl(xx);
  71   exx += scale;
  72     if (exx > LDBL_MAX_EXP)
  73     {
  74       errno = ERANGE;
  75       return INFINITY;
  76     }
  77   if (exx < LDBL_MIN_EXP)
  78     return 0.0L;
  79
  80   /* Undo scaling */
  81   return (scalbnl (xx, scale));
  82 }