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