Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

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