overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22
23/****************************************************************************
24 double cabs(double complex z) returns the absolute value (magnitude) of its
25 complex argument z, avoiding spurious overflow, underflow, and invalid
26 exceptions. The algorithm is from Kahan's paper.
27
28 CONSTANTS: FPKSQT2 = sqrt(2.0) to double precision
29 FPKR2P1 = sqrt(2.0) + 1.0 to double precision
30 FPKT2P1 = sqrt(2.0) + 1.0 - FPKR2P1 to double precision, so
31 that FPKR2P1 + FPKT2P1 = sqrt(2.0) + 1.0 to double
32 double precision.
33
34 Calls: fpclassify, fabs, sqrt, feholdexcept, fesetround, feclearexcept,
35 and feupdateenv.
36****************************************************************************/
37
38#include "math.h"
39#include "fenv.h"
40#include "fp_private.h"
41#include "complex.h"
42
43#define complex _Complex
44
45#define Real(z) (__real__ z)
46#define Imag(z) (__imag__ z)
47
48static const /* sqrt(2.0) */
49 hexdouble FPKSQT2 = HEXDOUBLE(0x3ff6a09e,0x667f3bcd);
50
51static const /* sqrt(2.0) + 1.0 to double */
52 hexdouble FPKR2P1 = HEXDOUBLE(0x4003504f,0x333f9de6);
53
54static const /* sqrt(2.0) + 1.0 - FPKR2P1 to double */
55 hexdouble FPKT2P1 = HEXDOUBLE(0x3ca21165,0xf626cdd6);
56
57static const
58 hexdouble infinity = HEXDOUBLE(0x7ff00000, 0x00000000);
59
60/****************************************************************************
61 double cabs(double complex z) returns the absolute value (magnitude) of its
62 complex argument z, avoiding spurious overflow, underflow, and invalid
63 exceptions. The algorithm is from Kahan's paper.
64
65 CONSTANTS: FPKSQT2 = sqrt(2.0) to double precision
66 FPKR2P1 = sqrt(2.0) + 1.0 to double precision
67 FPKT2P1 = sqrt(2.0) + 1.0 - FPKR2P1 to double precision, so
68 that FPKR2P1 + FPKT2P1 = sqrt(2.0) + 1.0 to double
69 double precision.
70
71 Calls: fpclassify, fabs, sqrt, feholdexcept, fesetround, feclearexcept,
72 and feupdateenv.
73****************************************************************************/
74
75double cabs ( double complex z )
76{
77 double a,b,s,t;
78 fenv_t env;
79 double FPR_inf = infinity.d;
80
81 a = fabs(Real(z));
82 b = fabs(Imag(z));
83
84 if (unlikely( (a == FPR_inf) || (b == FPR_inf) ))
85 return FPR_inf;
86
87 if (unlikely( (a != a) || (b != b) ))
88 return __FABS ( a + b );
89
90 if (unlikely((a == 0.0) || (b == 0.0) ))
91 return __FABS ( a + b );
92
93 /* both components of z are finite, nonzero */
94 {
95 (void)feholdexcept(&env); /* save environment, clear flags */
96 (void)fesetround(FE_TONEAREST); /* set default rounding */
97
98 s = 0.0;
99 if (a < b) /* order a >= b */
100 {
101 t = a;
102 a = b;
103 b = t;
104 }
105 t = a - b; /* magnitude difference */
106
107 if (t != a) /* b not negligible relative to a */
108 {
109 if (t > b) /* a - b > b */
110 {
111 s = a/b;
112 s += sqrt(1.0 + s*s);
113 }
114 else /* a - b <= b */
115 {
116 s = t/b;
117 t = (2.0 + s)*s;
118 s = ((FPKT2P1.d+t/(FPKSQT2.d+sqrt(2.0+t)))+s)+FPKR2P1.d;
119 }
120
121 s = b/s; /* may spuriously underflow */
122 feclearexcept(FE_UNDERFLOW);
123 }
124
125 feupdateenv(&env); /* restore environment */
126 return (a + s); /* deserved overflow occurs here */
127 } /* finite, nonzero case */
128}
129