Something went wrong. Try again.
Cider Isn't Darwin Emulation, Really
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655/* * 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@ *//******************************************************************************** ** File tableExpD.c, ** Functions exp, exp2 and expm1. ** Implementation of exp(x) based on IBM/Taligent table method. ** ** Copyright � 1997-2001 Apple Computer, Inc. All rights reserved. ** ** Written by Ali Sazegari, started on April 1997. ** Modified and ported by Robert A. Murley (ram) for Mac OS X. ** ** A MathLib v4 file. ** ** November 07 2001: removed exp2 to prevent conflict with CarbonCore. ** November 06 2001: commented out warning about Intel architectures. ** changed i386 stub to call abort(). ** November 02 2001: added stub for i386 version of exp2. ** April 28 1997: port of the ibm/taligent exp routines. ** July 16 1997: changed the rounding direction sensitivity of ** delivered results to default rounding at all times. ** August 28 2001: replaced __setflm with FEGETENVD/FESETENVD. ** replaced DblInHex typedef with hexdouble. ** used standard exception symbols from fenv.h. ** added #ifdef __ppc__. ** September 09 2001: added more comments. ** September 10 2001: added macros to detect PowerPC and correct compiler. ** September 18 2001: added <CoreServices/CoreServices.h> to get <fenv.h>. ** October 08 2001: removed <CoreServices/CoreServices.h>. ** changed compiler errors to warnings. ** ** W A R N I N G: ** These routines require a 64-bit double precision IEEE-754 model. ** They are written for PowerPC only and are expecting the compiler ** to generate the correct sequence of multiply-add fused instructions. ** ** These routines are not intended for 32-bit Intel architectures. ** ** A version of gcc higher than 932 is required. ** ** GCC compiler options: ** optimization level 3 (-O3) ** -fschedule-insns -finline-functions -funroll-all-loops ** ********************************************************************************/
#include "math.h"#include "fenv_private.h"#include "fp_private.h"
/******************************************************************************** Floating-point constants. ********************************************************************************/
static const double rintFactor = 6.7553994410557440000e+15; /* 0x43380000, 0x00000000 */static const double oneOverLn2 = 1.4426950408889633000e+00; /* 0x3ff71547, 0x652b82fe */static const double ln2Head = 6.9314718055994530000e-01; /* 0x3fe62e42, 0xfefa39ef */static const double ln2Tail = 2.3190468138462996000e-17; /* 0x3c7abc9e, 0x3b39803f */static const double maxExp = 7.0978271289338397000e+02; /* 0x40862e42, 0xfefa39ef */static const double denormal = 2.9387358770557188000e-39; /* 0x37f00000, 0x00000000 */static const hexdouble infinity = HEXDOUBLE(0x7ff00000, 0x00000000);
/******************************************************************************** Approximation coefficients. ********************************************************************************/
static const double cm1 = 8.3333336309523691e-03; /* 0x3f811111,0x1b4af38e */static const double c0 = 4.1666668402777808000e-02; /* 0x3fa55555,0x643f1505 */static const double c1 = 1.6666666666666655000e-01; /* 0x3fc55555,0x55555551 */static const double c2 = 4.9999999999999955000e-01; /* 0x3fdfffff,0xfffffff8 */
extern const uint32_t expTable[];
struct expTableEntry { double x; double f; };
static const double kMinNormal = 0x1.0p-1022; // 2.2250738585072014e-308;static const double kMaxNormal = 1.7976931348623157e308;
#if !defined(BUILDING_FOR_CARBONCORE_LEGACY)/******************************************************************************** Floating-point constants. ********************************************************************************/static const double minExp = -7.4513321910194122000e+02; /* 0xc0874910, 0xd52d3052 */static const double oneOverDenorm = 3.402823669209384635e+38; /* 0x47f00000, 0x00000000 */
/******************************************************************************** Approximation coefficients. ********************************************************************************/
static const double cc4 = 0.5000000000000000000;static const double cc3 = 0.16666666666663809522820;static const double cc2 = 0.04166666666666111110939;static const double cc1 = 0.008333338095239329810170;static const double cc0 = 0.001388889583333492938381;
/******************************************************************************** ** The base e exponential function. ** ********************************************************************************** ** Raised exceptions are inexact, overflow & inexact and underflow & inexact.** ********************************************************************************/
double exp ( double x ) { hexdouble scale, xInHex, yInHex; register double d, y, yTail, z, zTail, z2, temp1, temp2, power, result; register int32_t i; register struct expTableEntry *tablePointer, *pT; register double FPR_oneOverLn2, FPR_rintFactor, FPR_ln2Head, FPR_ln2Tail, FPR_h, FPR_z, FPR_512, FPR_scale; register double FPR_env, FPR_x, FPR_f, FPR_cm1, FPR_c0, FPR_c1, FPR_c2; FPR_x = __FABS( x ); FPR_z = 0.0; FPR_512 = 512.0; FPR_oneOverLn2 = oneOverLn2; FPR_rintFactor = rintFactor; FPR_ln2Head = ln2Head; FPR_ln2Tail = ln2Tail; tablePointer = ( struct expTableEntry * ) expTable + 177; FEGETENVD ( FPR_env ); // save old environment, set default __ENSURE( FPR_z, FPR_512, FPR_oneOverLn2 ); __ENSURE( FPR_ln2Head, FPR_ln2Tail, FPR_rintFactor ); FESETENVD ( FPR_z ); // Compute FPR_h and store it to yInHex as early as possible, even for cases that don't need it FPR_h = __FMADD( x, FPR_oneOverLn2, FPR_rintFactor ); yInHex.d = FPR_h; FPR_h -= FPR_rintFactor;
if (likely( FPR_x < FPR_512 )) { if (likely( x != FPR_z )) { scale.i.hi = ( yInHex.i.lo + 1023 ) << 20; scale.i.lo = 0;
y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h );
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
__NOOP; __NOOP; __NOOP; i = xInHex.i.lo; pT = &(tablePointer[i]); FPR_x = pT->x; FPR_f = pT->f; d = y - FPR_x;
z = d - yTail; FPR_cm1 = cm1; FPR_c1 = c1; FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 ); FPR_scale = scale.d;
temp2 = __FMADD( temp1, z2, z );
temp1 = __FMADD( zTail, temp2, zTail ); result = __FMUL( FPR_scale, FPR_f ); temp1 = temp1 + temp2;
result = __FMADD( result, temp1, result );
FESETENVD ( FPR_env ); __PROG_INEXACT( FPR_cm1 ); return result; } else { FESETENVD ( FPR_env ); return 1.0; } } if (likely( ( x <= maxExp ) && ( x > minExp ) )) { if ( x >= FPR_512 ) { power = 2.0; scale.i.hi = ( yInHex.i.lo + 1022 ) << 20; scale.i.lo = 0; } else { power = denormal; scale.i.hi = ( yInHex.i.lo + 1023+128 ) << 20; scale.i.lo = 0; }
y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h ); FPR_scale = scale.d;
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
__NOOP; __NOOP; __NOOP; i = xInHex.i.lo; pT = &(tablePointer[i]); FPR_x = pT->x; FPR_f = pT->f; d = y - FPR_x; result = __FMUL( FPR_scale, FPR_f );
z = d - yTail; FPR_cm1 = cm1; FPR_c1 = c1; FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 );
temp2 = __FMADD( temp1, z2, z );
temp1 = __FMADD( zTail, temp2, zTail ); temp1 = temp1 + temp2;
result = __FMADD( result, temp1, result ); result = __FMUL( result, power ); FESETENVD ( FPR_env ); __PROG_INEXACT( FPR_oneOverLn2 ); return result; } FESETENVD ( FPR_env ); if ( x != x ) return x; else if ( x == infinity.d ) return infinity.d; else if ( x == -infinity.d ) return FPR_z; else if ( x > maxExp ) { __PROG_OF_INEXACT( kMaxNormal ); return infinity.d; } else { __PROG_UF_INEXACT( kMinNormal ); return FPR_z; }}
/******************************************************************************** ** The expm1 function. ** ********************************************************************************/
static const hexdouble k2M55 = HEXDOUBLE(0x3c800000, 0x00000000);
double expm1 ( double x ) { hexdouble scale, invScale, xInHex, yInHex; register double d, y, yTail, z, zTail, z2, temp1, temp2, power, result, f; register int32_t i; register struct expTableEntry *tablePointer, *pT;
register double FPR_oneOverLn2, FPR_rintFactor, FPR_ln2Head, FPR_ln2Tail, FPR_h, FPR_z, FPR_512, FPR_scale; register double FPR_env, FPR_x, FPR_f, FPR_cc0, FPR_cc1, FPR_cc2, FPR_cc3, FPR_cc4, FPR_iscale; FPR_x = __FABS( x ); FPR_z = 0.0; FPR_512 = 512.0; FPR_oneOverLn2 = oneOverLn2; FPR_rintFactor = rintFactor; FPR_ln2Head = ln2Head; FPR_ln2Tail = ln2Tail; tablePointer = ( struct expTableEntry * ) expTable + 177;
FEGETENVD ( FPR_env ); // save old environment, set default __ENSURE( FPR_z, FPR_512, FPR_oneOverLn2 ); __ENSURE( FPR_ln2Head, FPR_ln2Tail, FPR_rintFactor ); FESETENVD ( FPR_z ); // Compute FPR_h and store it to yInHex as early as possible, even for cases that don't need it FPR_h = __FMADD( x, FPR_oneOverLn2, FPR_rintFactor ); yInHex.d = FPR_h; FPR_h -= FPR_rintFactor;
if (likely( FPR_x < FPR_512 )) { if (unlikely( FPR_x < k2M55.d )) { FESETENVD ( FPR_env ); if ( x == FPR_z ) { /* NOTHING */ } else { if ( FPR_x < kMinNormal ) __PROG_UF_INEXACT( kMinNormal ); else __PROG_INEXACT( FPR_oneOverLn2 ); } return x; } scale.i.hi = ( yInHex.i.lo + 1023 ) << 20; scale.i.lo = 0;
y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h ); xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor ); invScale.i.hi = 0x7fe00000 - scale.i.hi; invScale.i.lo = 0;
i = xInHex.i.lo; pT = &(tablePointer[i]); FPR_x = pT->x; FPR_f = pT->f; d = y - FPR_x; FPR_cc4 = cc4; z = d - yTail; FPR_cc0 = cc0; FPR_cc2 = cc2; FPR_cc1 = cc1; FPR_cc3 = cc3;
z2 = __FMUL( z, z ); zTail = d - z - yTail; temp1 = __FMADD( FPR_cc0, z2, FPR_cc2 ); temp2 = __FMADD( FPR_cc1, z2, FPR_cc3 ); temp1 = __FMADD( temp1, z2, FPR_cc4 ); temp2 = __FMADD( temp2, z, temp1 ); FPR_iscale = invScale.d; temp1 = __FMADD( temp2, z2, z ); temp2 = __FMADD( zTail, temp1, zTail ); temp2 = temp1 + temp2; d = FPR_f - FPR_iscale; FPR_scale = scale.d; temp1 = __FMADD( FPR_f, temp2, d ); result = __FMUL( temp1, FPR_scale ); FESETENVD ( FPR_env ); __PROG_INEXACT( FPR_oneOverLn2 ); return result; } if (likely( ( x <= maxExp ) && ( x > minExp ) )) { if ( x >= FPR_512 ) { power = 2.0; f = 0.5; scale.i.hi = ( yInHex.i.lo + 1022 ) << 20; scale.i.lo = 0; } else { power = denormal; f = oneOverDenorm; scale.i.hi = ( yInHex.i.lo + 1023+128 ) << 20; scale.i.lo = 0; if ( scale.i.hi < ( 168<<20 ) ) { FESETENVD ( FPR_env ); __PROG_INEXACT( FPR_oneOverLn2 ); return -1.0; } } y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h ); FPR_scale = scale.d;
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor ); invScale.i.hi = 0x7fe00000 - scale.i.hi; invScale.i.lo = 0;
i = xInHex.i.lo; pT = &(tablePointer[i]); FPR_x = pT->x; FPR_f = pT->f; d = y - FPR_x; FPR_iscale = invScale.d; FPR_cc4 = cc4;
z = d - yTail; FPR_cc0 = cc0; FPR_cc2 = cc2; FPR_cc1 = cc1; FPR_cc3 = cc3;
z2 = __FMUL( z, z ); zTail = d - z - yTail; temp1 = __FMADD( FPR_cc0, z2, FPR_cc2 ); temp2 = __FMADD( FPR_cc1, z2, FPR_cc3 ); temp1 = __FMADD( temp1, z2, FPR_cc4 ); d = __FNMSUB( f, FPR_iscale, FPR_f ); temp2 = __FMADD( temp2, z, temp1 ); temp1 = __FMADD( temp2, z2, z ); temp2 = __FMADD( zTail, temp1, zTail ); temp2 = temp1 + temp2; temp1 = __FMADD( FPR_f, temp2, d ); result = __FMUL( temp1, FPR_scale );
result = __FMUL( result, power ); FESETENVD ( FPR_env ); __PROG_INEXACT( FPR_oneOverLn2 ); return result; } FESETENVD ( FPR_env ); if ( x != x ) return x; else if ( x == infinity.d ) return infinity.d; else if ( x == -infinity.d ) return -1.0; else if ( x > maxExp ) { __PROG_OF_INEXACT( kMaxNormal ); return infinity.d; } else {#if 0 /* XXX scp: test vectors don't want this to underflow */ __PROG_UF_INEXACT( kMinNormal );#else __PROG_INEXACT( FPR_oneOverLn2 );#endif return -1.0; }}
#else /* BUILDING_FOR_CARBONCORE_LEGACY *//******************************************************************************** Floating-point constants. ********************************************************************************/
static const double maxExp2 = 1024.0; static const double minNormExp2 = -1022.0;static const double minExp2 = -1075.0;
/******************************************************************************** ** The base 2 exponential function. ** ********************************************************************************/
double exp2 ( double x ) { hexdouble OldEnvironment, scale, xInHex, yInHex; register double d, y, yTail, z, zTail, z2, temp1, temp2, power, result; register int32_t i; register struct expTableEntry *tablePointer, *pT; register double FPR_oneOverLn2, FPR_rintFactor, FPR_ln2Head, FPR_ln2Tail, FPR_h, FPR_z, FPR_512, FPR_scale; register double FPR_env, FPR_diff, FPR_x, FPR_f, FPR_cm1, FPR_c0, FPR_c1, FPR_c2; FPR_x = __FABS( x ); FPR_z = 0.0; FPR_512 = 512.0; FPR_oneOverLn2 = oneOverLn2; FPR_rintFactor = rintFactor; FPR_ln2Head = ln2Head; FPR_ln2Tail = ln2Tail; tablePointer = ( struct expTableEntry * ) expTable + 177; FEGETENVD ( FPR_env ); // save old environment, set default __ENSURE( FPR_z, FPR_512, FPR_oneOverLn2 ); __ENSURE( FPR_ln2Head, FPR_ln2Tail, FPR_rintFactor ); FESETENVD ( FPR_z ); // Compute FPR_h and store it to yInHex as early as possible, even for cases that don't need it FPR_h = x + FPR_rintFactor; yInHex.d = FPR_h; FPR_h -= FPR_rintFactor; FPR_diff = x - FPR_h; if (likely( FPR_x < FPR_512 )) { if (likely( FPR_x != FPR_z )) { scale.i.hi = ( yInHex.i.lo + 1023 ) << 20; scale.i.lo = 0; y = __FMUL( FPR_ln2Head, FPR_diff ); yTail = __FMUL( FPR_ln2Tail, FPR_diff ); xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor ); __NOOP; __NOOP; __NOOP; i = xInHex.i.lo; pT = &(tablePointer[i]); FPR_x = pT->x; FPR_f = pT->f; d = y - FPR_x;
z = d - yTail; FPR_cm1 = cm1; FPR_c1 = c1; FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 );
temp2 = __FMADD( temp1, z2, z ); FPR_scale = scale.d;
temp1 = __FMADD( zTail, temp2, zTail ); temp1 = temp1 + temp2; result = __FMUL( FPR_scale, FPR_f );
result = __FMADD( result, temp1, result ); FESETENVD ( FPR_env ); if ( FPR_diff != FPR_z) __PROG_INEXACT( FPR_oneOverLn2 ); return result; } else { FESETENVD ( FPR_env ); return 1.0; } }
if (likely( ( x < maxExp2 ) && ( x > minExp2 ) )) { if ( x >= FPR_512 ) { power = 2.0; scale.i.hi = ( yInHex.i.lo + 1022 ) << 20; scale.i.lo = 0; } else { power = denormal; scale.i.hi = ( yInHex.i.lo + 1023+128 ) << 20; scale.i.lo = 0; }
FPR_diff = x - FPR_h; FPR_scale = scale.d; y = __FMUL( FPR_ln2Head, FPR_diff ); yTail = __FMUL( FPR_ln2Tail, FPR_diff ); xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
__NOOP; __NOOP; __NOOP; i = xInHex.i.lo; pT = &(tablePointer[i]); FPR_x = pT->x; FPR_f = pT->f; d = y - FPR_x; result = __FMUL( FPR_scale, FPR_f );
z = d - yTail; FPR_cm1 = cm1; FPR_c1 = c1; FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 );
temp2 = __FMADD( temp1, z2, z );
temp1 = __FMADD( zTail, temp2, zTail ); temp1 = temp1 + temp2;
result = __FMADD( result, temp1, result ); result = __FMUL( result, power ); FESETENVD ( FPR_env ); if ( x < minNormExp2 ) { OldEnvironment.d = FPR_env; __NOOP; __NOOP; __NOOP; OldEnvironment.i.lo |= FE_UNDERFLOW; // set underflow flag if ( FPR_h != x ) OldEnvironment.i.lo |= FE_INEXACT; // set inexact flag FESETENVD_GRP ( OldEnvironment.d ); } else if ( FPR_h != x ) __PROG_INEXACT( FPR_oneOverLn2 );
return result; } FESETENVD ( FPR_env ); if ( x != x ) return x; else if ( x == infinity.d ) return infinity.d; else if ( x == -infinity.d ) return FPR_z; else if ( x > maxExp ) { __PROG_OF_INEXACT( kMaxNormal ); return infinity.d; } else { __PROG_UF_INEXACT( kMinNormal ); return FPR_z; }}
#endif /* BUILDING_FOR_CARBONCORE_LEGACY */