Something went wrong. Try again.
Cider Isn't Darwin Emulation, Really
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130/* * Copyright (c) 2002 Apple Computer, Inc. All rights reserved. * * @APPLE_LICENSE_HEADER_START@ * * The contents of this file constitute Original Code as defined in and * are subject to the Apple Public Source License Version 1.1 (the * "License"). You may not use this file except in compliance with the * License. Please obtain a copy of the License at * http://www.apple.com/publicsource and read it before using this file. * * This Original Code and all software distributed under the License are * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES, * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the * License for the specific language governing rights and limitations * under the License. * * @APPLE_LICENSE_HEADER_END@ */
/**************************************************************************** double cabs(double complex z) returns the absolute value (magnitude) of its complex argument z, avoiding spurious overflow, underflow, and invalid exceptions. The algorithm is from Kahan's paper. CONSTANTS: FPKSQT2 = sqrt(2.0) to double precision FPKR2P1 = sqrt(2.0) + 1.0 to double precision FPKT2P1 = sqrt(2.0) + 1.0 - FPKR2P1 to double precision, so that FPKR2P1 + FPKT2P1 = sqrt(2.0) + 1.0 to double double precision. Calls: fpclassify, fabs, sqrt, feholdexcept, fesetround, feclearexcept, and feupdateenv.****************************************************************************/
#include "math.h"#include "fenv.h"#include "fp_private.h"#include "complex.h"
#define complex _Complex
#define Real(z) (__real__ z)#define Imag(z) (__imag__ z)
static const /* sqrt(2.0) */ hexdouble FPKSQT2 = HEXDOUBLE(0x3ff6a09e,0x667f3bcd);
static const /* sqrt(2.0) + 1.0 to double */ hexdouble FPKR2P1 = HEXDOUBLE(0x4003504f,0x333f9de6);
static const /* sqrt(2.0) + 1.0 - FPKR2P1 to double */ hexdouble FPKT2P1 = HEXDOUBLE(0x3ca21165,0xf626cdd6);
static const hexdouble infinity = HEXDOUBLE(0x7ff00000, 0x00000000);
/**************************************************************************** double cabs(double complex z) returns the absolute value (magnitude) of its complex argument z, avoiding spurious overflow, underflow, and invalid exceptions. The algorithm is from Kahan's paper. CONSTANTS: FPKSQT2 = sqrt(2.0) to double precision FPKR2P1 = sqrt(2.0) + 1.0 to double precision FPKT2P1 = sqrt(2.0) + 1.0 - FPKR2P1 to double precision, so that FPKR2P1 + FPKT2P1 = sqrt(2.0) + 1.0 to double double precision. Calls: fpclassify, fabs, sqrt, feholdexcept, fesetround, feclearexcept, and feupdateenv.****************************************************************************/
double cabs ( double complex z ){ double a,b,s,t; fenv_t env; double FPR_inf = infinity.d; a = fabs(Real(z)); b = fabs(Imag(z)); if (unlikely( (a == FPR_inf) || (b == FPR_inf) )) return FPR_inf; if (unlikely( (a != a) || (b != b) )) return __FABS ( a + b ); if (unlikely((a == 0.0) || (b == 0.0) )) return __FABS ( a + b );
/* both components of z are finite, nonzero */ { (void)feholdexcept(&env); /* save environment, clear flags */ (void)fesetround(FE_TONEAREST); /* set default rounding */
s = 0.0; if (a < b) /* order a >= b */ { t = a; a = b; b = t; } t = a - b; /* magnitude difference */ if (t != a) /* b not negligible relative to a */ { if (t > b) /* a - b > b */ { s = a/b; s += sqrt(1.0 + s*s); } else /* a - b <= b */ { s = t/b; t = (2.0 + s)*s; s = ((FPKT2P1.d+t/(FPKSQT2.d+sqrt(2.0+t)))+s)+FPKR2P1.d; } s = b/s; /* may spuriously underflow */ feclearexcept(FE_UNDERFLOW); } feupdateenv(&env); /* restore environment */ return (a + s); /* deserved overflow occurs here */ } /* finite, nonzero case */}