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/*
11This implementation is based largely on Cephes library
12function cabsl (cmplxl.c), which bears the following notice:
13
14Cephes Math Library Release 2.1:  January, 1989
15Copyright 1984, 1987, 1989 by Stephen L. Moshier
16Direct 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
28long double
29hypotl (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}