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
hypotl(long double x,long double y)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 }
83