Something went wrong. Try again.
Cider Isn't Darwin Emulation, Really
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393/* * expf.s * * by Ian Ollmann * * Copyright (c) 2007, Apple Inc. All Rights Reserved. * * This file implements the C99 expf function for the MacOS X __i386__ and __x86_64__ abis. */
#include <machine/asm.h>#include "abi.h"
//// Overall approach//// We calculate e**(y) as follows://// e**y = 2**x x = y / ln(2)//// ...and proceed largely as we do in exp2f()://// i = trunc(x)// f = x - i// // 2**(x) = 2**(f+i) = 2**f * 2**i//// We choose trunc here, so that 2**f always has the property of pushing 2**i that much closer to the // extrema of finite, non-zero floating point computation. That is://// x < 1 --> 2**f <= 1// x == 1 --> 2**f = 1// x > 1 --> 2**f >= 1//// This means that our 2**i calculation will not underflow/overflow prematurely before 2**f has had its say.// Round to nearest even would produce somewhat more accurate results, but accomplishing something like that // rounding without knowledge of the current rounding mode is expensive.// // In principle, the 2**f part can further be broken up as://// 2**f = 2**(c+r) = 2**c * 2**r//// Where c is chosen so as to make r as close to 0 as possible and c can be determined fairly exactly by // consulting a lookup table. Our choice of trunc above means we cant sure of the sign of f,// so our table ends up being twice as large as it would have if we had chosen floor or ceil. On the bright // side, we never have to worry about f = x - trunc(f) being inexact in single precision.//// Unfortunately, in practice this reduction seemed to be bad for throughput, so this implementation went with// a large polynomial instead that covers the range -1,1 or -0.5,0.5 depending on whether we know the rounding// mode is nearest or not.//// Most of the arithmetic will be done in double precision. I investigated calculating the fractional term as 2**r-1, so as // to preserve precision -- we expect a lot of zeros between the 1 and the first non-zero bit. Normally, we dont want// throw that away until the final operation so that non-default rounding modes arrive at something close to the // right answer. Doing float in double precision means that this is unnecessary however. //// This code is probably worth revisiting once more to see if we can come up with a cheap reduction. These huge// polynomials just arent good for latency.////
.align 4// 8th order minimax fit of exp2 on [-1.0,1.0]. |error| < 0.402865722354948566583852e-9:expf_c: .quad 0x40bc03f30399c376, 0x3ff000000001ea2a // c4/c8 = 0.961813690023115610862381719985771e-2 / 0.134107709538786543922336536865157e-5, c0 = 1.0 + 0.278626872016317130037181614004e-10 .quad 0x408f10e7f73e6d8f, 0x3fe62e42fd0933ee // c5/c8 = 0.133318252930790403741964203236548e-2 / 0.134107709538786543922336536865157e-5, c1 = .693147176943623740308984004029708 .quad 0x405cb616a9384e69, 0x3fcebfbdfd0f0afa // c6/c8 = 0.154016177542147239746127455226575e-3 / 0.134107709538786543922336536865157e-5, c2 = .240226505817268621584559118975830 .quad 0x4027173ebd288ba1, 0x3fac6b0a74f15403 // c7/c8 = 0.154832722143258821052933667742417e-4 / 0.134107709538786543922336536865157e-5, c3 = 0.555041568519883074165425891257052e-1 .quad 0x3eb67fe1dc3105ba // c8 = 0.134107709538786543922336536865157e-5
.align 3expf_nofenv_c: .quad 0x3ff0000000058fca // 1.0 + -8.09329727503262700660348520172e-11 .double 0.693147206709644041218074094717934 .double 0.240226515050550309232521176082490 .double 0.0555032721577134480296332498187050 .double 0.00961799451418197772157647729048837 .double 0.00134004316655903280893023589328151 .double 0.000154780222739739895295319810010147
.literal8reciprocalLn2: .quad 0x3ff71547652b82fe // 1.0 / ln(2)
.text
#if defined( __x86_64__ ) #define RELATIVE_ADDR( _a) (_a)( %rip ) #define RELATIVE_ADDR_B( _a) (_a)( %rip )#elif defined( __i386__ ) #define RELATIVE_ADDR( _a) (_a)-expf_body( CX_P ) #define RELATIVE_ADDR_B( _a) (_a)-expf_no_fenv_body( CX_P )
//a short routine to get the local address.align 4expf_pic: movl (%esp), %ecx //copy address of local_addr to %ecx ret#else #error arch not supported#endif
//Accurate to within 0.5068 ulp in round to nearest. ENTRY( expf )#if defined( __i386__ ) movl FRAME_SIZE(STACKP), %eax movss FRAME_SIZE(STACKP), %xmm0#else movd %xmm0, %eax#endif cvtss2sd %xmm0, %xmm2 // (double) x movl %eax, %ecx // x andl $0x7fffffff, %eax // |x|vl movl $1023, %edx // double precision bias subl $0x33b8aa3b, %eax // subtract 0x1.0p-24f / ln(2) cmpl $0x0f7d1d57, %eax // if( |x| > 126/ln(2) || isnan(x) || |x| < 0x1.0p-26f /ln(2) ) ja 3f // goto 3
//The test above was a little too aggressive. //values 126 < x < 128, and -150 < x < -126 will come back here: //For x < -126, we change %edx to hold the exponent for 2**-126, // and subtract -126 from x
1://PIC#if defined( __i386__ ) calll expf_pic // set %ecx to point to local_addrexpf_body: #endif
mulsd RELATIVE_ADDR( reciprocalLn2 ), %xmm2 lea RELATIVE_ADDR( expf_c), CX_P
//extract fractional part and exit early if it is zero cvttsd2si %xmm2, %eax // trunc(x) addl %eax, %edx // calculate biased exponent of result cvtsi2sd %eax, %xmm1 // trunc(x) movd %edx, %xmm7 // exponent part of result psllq $52, %xmm7 // shift result exponent into place subsd %xmm1, %xmm2 // get the fractional part
// c0 + c1x1 + c2x2 + c3x3 + c4x4 + c5x5 + c6x6 + c7x7 + c8x8 #if defined( __SSE3__ ) movddup %xmm2, %xmm1 // { x, x }#else movapd %xmm2, %xmm1 // x unpcklpd %xmm1, %xmm1 // { x, x }#endif mulsd %xmm2, %xmm2 // x*x movapd %xmm1, %xmm3 mulpd 48(CX_P), %xmm1 // { c3x, (c7/c8)x } mulpd 16(CX_P), %xmm3 // { c1x, (c5/c8)x }#if defined( __SSE3__ ) movddup %xmm2, %xmm4 // { xx, xx }#else movapd %xmm2, %xmm4 // xx unpcklpd %xmm4, %xmm4 // { xx, xx }#endif mulsd %xmm2, %xmm2 // xx*xx addpd 32(CX_P), %xmm1 // { c2 + c3x, (c6/c8) + (c7/c8)x } addpd (CX_P), %xmm3 // { c0 + c1x, (c4/c8) + (c5/c8)x } mulpd %xmm4, %xmm1 // { c2xx + c3xxx, (c6/c8)xx + (c7/c8)xxx } addsd %xmm2, %xmm3 // { c0 + c1x, (c4/c8) + (c5/c8)x + xxxx } mulsd 64(CX_P), %xmm2 // c8 * xxxx addpd %xmm1, %xmm3 // { c0 + c1x + c2xx + c3xxx, (c4/c8) + (c5/c8)x + (c6/c8)xx + (c7/c8)xxx + xxxx } movhlps %xmm3, %xmm6 // { ?, c0 + c1x + c2xx + c3xxx } mulsd %xmm2, %xmm3 // { ..., c8xxxx* ((c4/c8) + (c5/c8)x + (c6/c8)xx + (c7/c8)xxx + xxxx) } addsd %xmm6, %xmm3 // c0 + c1x + c2xx + c3xxx + c4xxxx + c5xxxxx + c6xxxxxx + c7xxxxxxx + c8xxxxxxxxx mulsd %xmm7, %xmm3 // 2**i * {c0 + c1x + c2xx + c3xxx + c4xxxx + c5xxxxx + c6xxxxxx + c7xxxxxxx + c8xxxxxxxxx} // convert to single precision and return cvtsd2ss %xmm3, %xmm0#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret
3: // |x| > 126 || isnan(x) || |x| < 0x1.0p-26f jge 4f // if( |x| > 126 || isnan(x) ) goto 4
//fall through for common case of |x| < 0x1.0p-26f movl $0x3f800000, %edx // 1.0f movd %edx, %xmm1 // 1.0f addss %xmm1, %xmm0 // round away from 1.0 as appropriate for current rounding mode#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret
// |x| > 126 || isnan(x)4: addl $0x33b8aa3b, %eax // restore |x| cmpl $0x7f800000, %eax // ja 5f // if isnan(x) goto 5 je 6f // if isinf(x) goto 6 cmpl $0xc358677d, %ecx // if( x <= -150 /ln(2) ) jae 7f // goto 7
cmpl $0x4338aa3b, %ecx // if( x >= 128 /ln(2) ) jge 8f // goto 8
jmp 1b
5: //x is nan addss %xmm0, %xmm0 //Quiet the NaN#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret 6: // |x| is inf movss %xmm0, %xmm1 // x psrad $31, %xmm0 // x == -Inf ? -1U : 0 andnps %xmm1, %xmm0 // x == -Inf ? 0 : Inf#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret 7: // x <= -150 movl $0x00800001, %eax movd %eax, %xmm0 mulss %xmm0, %xmm0 //Create correctly rounded result, set inexact/underflow#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret 8: // x >= 128.0 movl $0x7f7fffff, %eax // FLT_MAX movd %eax, %xmm0 // FLT_MAX addss %xmm0, %xmm0 // Inf, set overflow#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret
//Uses round to nearest to deliver a reduced value between [ -0.5, 0.5 ]//This allows us to use a smaller polynomial//In addition, we dont bother to check for x is integer. This saves a few cycles in the general case.//Accurate to within 0.534 ulpENTRY( expf$fenv_access_off )#if defined( __i386__ ) movl FRAME_SIZE(STACKP), %eax movss FRAME_SIZE(STACKP), %xmm0#else movd %xmm0, %eax#endif cvtss2sd %xmm0, %xmm2 // (double) x movl %eax, %ecx // x andl $0x7fffffff, %eax // |x|vl movl $1023, %edx // double precision bias subl $0x33b8aa3b, %eax // subtract 0x1.0p-24f / ln(2) cmpl $0x0f7d1d57, %eax // if( |x| > 126/ln(2) || isnan(x) || |x| < 0x1.0p-26f /ln(2) ) ja 3f // goto 3
//The test above was a little too aggressive. //values 126 < x < 128, and -150 < x < -126 will come back here: //For x < -126, we change %edx to hold the exponent for 2**-126, // and subtract -126 from x
1://PIC#if defined( __i386__ ) calll expf_pic // set %ecx to point to local_addrexpf_no_fenv_body: #endif
mulsd RELATIVE_ADDR( reciprocalLn2 ), %xmm2 lea RELATIVE_ADDR_B( expf_nofenv_c), CX_P
//extract fractional part and exit early if it is zero andl $0x80000000, %ecx // signof(x) orl $0x43300000, %ecx // copysign( 0x1.0p52, x ) >> 32 movd %ecx, %xmm6 // copysign( 0x1.0p52, x ) >> 32 movaps %xmm2, %xmm1 // x psllq $32, %xmm6 // copysign( 0x1.0p52, x ) addsd %xmm6, %xmm2 // x + copysign( 0x1.0p23f, x ) subsd %xmm6, %xmm2 // rint(x) cvttsd2si %xmm2, %eax // rint(x) addl %eax, %edx // rint(x) + bias movd %edx, %xmm7 // result exponent >> 32 subsd %xmm2, %xmm1 // fractional part psllq $52, %xmm7 // shift result exponent into place
movsd (3*8)(CX_P), %xmm3 // c3 movsd (5*8)(CX_P), %xmm5 // c5 movsd (6*8)(CX_P), %xmm6 // c6 movapd %xmm1, %xmm2 // x mulsd %xmm1, %xmm1 // x*x mulsd %xmm2, %xmm3 // c3*x mulsd %xmm2, %xmm5 // c5*x mulsd (1*8)(CX_P), %xmm2 // c1*x mulsd %xmm1, %xmm6 // c6*xx addsd (2*8)(CX_P), %xmm3 // c2+c3x movapd %xmm1, %xmm4 // xx mulsd %xmm1, %xmm1 // xx*xx addsd (4*8)(CX_P), %xmm5 // c4+c5x addsd (CX_P), %xmm2 // c0+c1x mulsd %xmm4, %xmm3 // c2xx+c3xxx addsd %xmm6, %xmm5 // c4+c5x+c6xx mulsd %xmm1, %xmm5 // c4xxxx+c5xxxxx+c6xxxxxx addsd %xmm2, %xmm3 // c0+c1x+c2xx+c3xxx addsd %xmm5, %xmm3 // c0+c1x+c2xx+c3xxx+c4xxxx+c5xxxxx+c6xxxxxx
mulsd %xmm7, %xmm3 // scale by 2**i // convert to single precision and return cvtsd2ss %xmm3, %xmm0#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret
3: // |x| > 126 || isnan(x) || |x| < 0x1.0p-26f jge 4f // if( |x| > 126 || isnan(x) ) goto 4
//fall through for common case of |x| < 0x1.0p-26f movl $0x3f800000, %edx // 1.0f movd %edx, %xmm1 // 1.0f addss %xmm1, %xmm0 // round away from 1.0 as appropriate for current rounding mode#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret
// |x| > 126 || isnan(x)4: addl $0x33b8aa3b, %eax // restore |x| cmpl $0x7f800000, %eax // ja 5f // if isnan(x) goto 5 je 6f // if isinf(x) goto 6 cmpl $0xc358677d, %ecx // if( x <= -150 /ln(2) ) jae 7f // goto 7
cmpl $0x4338aa3b, %ecx // if( x >= 128 /ln(2) ) jge 8f // goto 8
jmp 1b
5: //x is nan addss %xmm0, %xmm0 //Quiet the NaN#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret 6: // |x| is inf movss %xmm0, %xmm1 // x psrad $31, %xmm0 // x == -Inf ? -1U : 0 andnps %xmm1, %xmm0 // x == -Inf ? 0 : Inf#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret 7: // x <= -150 movl $0x00800001, %eax movd %eax, %xmm0 mulss %xmm0, %xmm0 //Create correctly rounded result, set inexact/underflow#if defined( __i386__ ) movss %xmm0, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#endif ret 8: // x >= 128.0 movl $0x7f800000, %eax // Inf#if defined( __i386__ ) movl %eax, FRAME_SIZE(STACKP) flds FRAME_SIZE(STACKP)#else movd %eax, %xmm0#endif ret