Something went wrong. Try again.
Cider Isn't Darwin Emulation, Really
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718
/* * powf.s * * by Ian Ollmann * * Copyright (c) 2007, Apple Inc. All Rights Reserved. * * Implementation of C99 powf function for MacOS X __i386__ and __x86_64__ architectures. * */
#define LOCAL_STACK_SIZE 3*FRAME_SIZE #include "machine/asm.h"#include "abi.h"
.constgMaskShift: .byte 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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,31,31,31,31,31,31,31,31, \ 31,31,31,31, 31,31,31,31,31,31,31,31,31,31,31,31, \ 31,31,31,31, 31,31,31,31,31,31,31,31,31,31,31,31, \ 31,31,31,31, 31,31,31,31,31,31,31,31,31,31,31,31, \ 31,31,31,31, 31,31,31,31,31,31,31,31,31,31,31,31, \ 31,31,31,31, 31,31,31,31,31,31,31,31,31,31,31,31, \ 31,31,31,31, 31,31,31,31,31,31,31,31,31,31,31,31
.align 4// 8th order minimax fit of exp2 on [-1.0,1.0]. |error| < 0.402865722354948566583852e-9:powf_exp_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
// The reduction for the log2 stage is done as://// for log2(x):// x = 2**i * 1.m 1.0 <= 1.m < 2.0// index = top 7 bits of m// reduced = 1.m * powf_log2_c[2*index] 1-2**-8 < reduced <= 1+2**-7// log2x = exp2(i) + powf_log2_c[2*index+1] + log1p(reduced - 1) / ln(2);//// exp2(i) is exact// powf_log2_c[2*index] has 53 bits of precision, and is exact for the first and last entries// powf_log2_c[2*index+1] has 53 bits of precision, and is exact for the first and last entries// log1p( -2**-8 < x <= 2**-7 ) is done as a 5 term Taylor series. Error should be proportional// to the missing 6th order term or < 2**(6*-7)/6 ~ 2**-44, in the worst case. For smaller x, // it will obviously be better. // // Error in powf is rougly porportional to error_in_exp2_stage + y * error_in_log2_stage.// y is bounded by the point that powf overflows (underflow loses precision a bit sooner) // For our 2**-7 worst case, y is about 10000 at overflow. 10000 * 2**-44 = ~6e-10, ~2**-31// which should give us a good margin of safety. For cases that are in 1.0-2**-8 < x < 1.0+2**-7, // we expect the precision of the Taylor series to improve faster than y's ability to disrupt // the precision in ylog2(x) up to the point that we run out of precision in the double. The // worst case in this regard is probably 0x1.ffffep-1f, the number closest to 1. This overflows // somewhere in the range y = [ -2**30, -2**31 ]. We predict we'll need 53-55 bits of precision// here, which may slightly exceed the precision of the double. While this number may come out // wrong by a few ulps, all the other ones should be within tolerance since the next closest // number is twice as far from 1.0. For these cases, the values in powf_log2_c are exact, so // the only source of error is the Taylor series for log1p, and post-scaling by 1/ln(2).// For cases not in 1.0-2**-8 < x < 1.0+2**-7 the point at which y overflows is much smaller, so // we don't need so much precision.//// Reduction table for log2 stage of power prepared as:// #include <stdio.h>// #include <stdint.h>// #include <math.h>//// int main( void )// {// int i;//// for( i = 0; i < 128; i++ )// {// long double a = 1.0L /( 1.0L + (long double) i / (long double) 127 );// union{ double d; uint64_t u;}u, v;// u.d = a;// v.d = -log2l(a);//// printf( ".quad 0x%llx,\t0x%llx\t// %Lg, -log2l(%Lg)\n", u.u, v.u, a, a );// }//// return 0;// }//// We use /127 here rather than /128 to allow the cases where the most precision is needed to be reduced by// exact powers of two. (These are 1.0 + 1 ulp and 1.0 - 1ulp.) The other values will land somewhere in the// range [ 1.0 - 2**-8, 1.0 + 2**-7 ]. (Experimentally verified for all floats 1.0 <= x < 2.0.)//.align 3powf_log2_c: .quad 0x3ff0000000000000, 0x8000000000000000 // 1, -log2l(1) .quad 0x3fefc00000000000, 0x3f872c7ba20f7327 // 0.992188, -log2l(0.992188) .quad 0x3fef80fe03f80fe0, 0x3f9715662c7f3dbc // 0.984496, -log2l(0.984496) .quad 0x3fef42f42f42f42f, 0x3fa13eea2b6545df // 0.976923, -log2l(0.976923) .quad 0x3fef05dcd30dadec, 0x3fa6e7f0bd9710dd // 0.969466, -log2l(0.969466) .quad 0x3feec9b26c9b26ca, 0x3fac85f25e12da51 // 0.962121, -log2l(0.962121) .quad 0x3fee8e6fa39be8e7, 0x3fb10c8cd0c74414 // 0.954887, -log2l(0.954887) .quad 0x3fee540f4898d5f8, 0x3fb3d0c813e48e00 // 0.947761, -log2l(0.947761) .quad 0x3fee1a8c536fe1a9, 0x3fb68fbf5169e028 // 0.940741, -log2l(0.940741) .quad 0x3fede1e1e1e1e1e2, 0x3fb949866f0b017b // 0.933824, -log2l(0.933824) .quad 0x3fedaa0b3630957d, 0x3fbbfe30e28821c0 // 0.927007, -log2l(0.927007) .quad 0x3fed7303b5cc0ed7, 0x3fbeadd1b4ef9a1f // 0.92029, -log2l(0.92029) .quad 0x3fed3cc6e80ebbdb, 0x3fc0ac3dc2e0ca0c // 0.913669, -log2l(0.913669) .quad 0x3fed075075075075, 0x3fc1ff2046fb7116 // 0.907143, -log2l(0.907143) .quad 0x3fecd29c244fe2f3, 0x3fc34f99517622ae // 0.900709, -log2l(0.900709) .quad 0x3fec9ea5dbf193d5, 0x3fc49db19c99a54d // 0.894366, -log2l(0.894366) .quad 0x3fec6b699f5423ce, 0x3fc5e971b3a4ee80 // 0.888112, -log2l(0.888112) .quad 0x3fec38e38e38e38e, 0x3fc732e1f41ccdba // 0.881944, -log2l(0.881944) .quad 0x3fec070fe3c070fe, 0x3fc87a0a8f0ff9b2 // 0.875862, -log2l(0.875862) .quad 0x3febd5eaf57abd5f, 0x3fc9bef38a4ffae5 // 0.869863, -log2l(0.869863) .quad 0x3feba5713280dee9, 0x3fcb01a4c19f6811 // 0.863946, -log2l(0.863946) .quad 0x3feb759f2298375a, 0x3fcc4225e7d5e3c6 // 0.858108, -log2l(0.858108) .quad 0x3feb4671655e7f24, 0x3fcd807e87fa4521 // 0.852349, -log2l(0.852349) .quad 0x3feb17e4b17e4b18, 0x3fcebcb6065350a2 // 0.846667, -log2l(0.846667) .quad 0x3feae9f5d3eba7d7, 0x3fcff6d3a16f617f // 0.84106, -log2l(0.84106) .quad 0x3feabca1af286bca, 0x3fd0976f3991af9e // 0.835526, -log2l(0.835526) .quad 0x3fea8fe53a8fe53b, 0x3fd1326eb8c0aba3 // 0.830065, -log2l(0.830065) .quad 0x3fea63bd81a98ef6, 0x3fd1cc6bb7e3870f // 0.824675, -log2l(0.824675) .quad 0x3fea3827a3827a38, 0x3fd265698fa26c0a // 0.819355, -log2l(0.819355) .quad 0x3fea0d20d20d20d2, 0x3fd2fd6b881e82d3 // 0.814103, -log2l(0.814103) .quad 0x3fe9e2a65187566c, 0x3fd39474d95e1649 // 0.808917, -log2l(0.808917) .quad 0x3fe9b8b577e61371, 0x3fd42a88abb54986 // 0.803797, -log2l(0.803797) .quad 0x3fe98f4bac46d7c0, 0x3fd4bfaa182b7fe3 // 0.798742, -log2l(0.798742) .quad 0x3fe9666666666666, 0x3fd553dc28dd9724 // 0.79375, -log2l(0.79375) .quad 0x3fe93e032e1c9f02, 0x3fd5e721d95d124d // 0.78882, -log2l(0.78882) .quad 0x3fe9161f9add3c0d, 0x3fd6797e170c5221 // 0.783951, -log2l(0.783951) .quad 0x3fe8eeb9533d4065, 0x3fd70af3c177f740 // 0.779141, -log2l(0.779141) .quad 0x3fe8c7ce0c7ce0c8, 0x3fd79b85aaad8878 // 0.77439, -log2l(0.77439) .quad 0x3fe8a15b8a15b8a1, 0x3fd82b36978f76d5 // 0.769697, -log2l(0.769697) .quad 0x3fe87b5f9d4d1bc2, 0x3fd8ba09402697ed // 0.76506, -log2l(0.76506) .quad 0x3fe855d824ca58e9, 0x3fd948004ff12dbf // 0.760479, -log2l(0.760479) .quad 0x3fe830c30c30c30c, 0x3fd9d51e662f92a2 // 0.755952, -log2l(0.755952) .quad 0x3fe80c1e4bbd595f, 0x3fda6166162e9ec8 // 0.751479, -log2l(0.751479) .quad 0x3fe7e7e7e7e7e7e8, 0x3fdaecd9e78fdbea // 0.747059, -log2l(0.747059) .quad 0x3fe7c41df1077c42, 0x3fdb777c568f9ae2 // 0.74269, -log2l(0.74269) .quad 0x3fe7a0be82fa0be8, 0x3fdc014fd448fe3a // 0.738372, -log2l(0.738372) .quad 0x3fe77dc7c4cf2aea, 0x3fdc8a56c6f80bca // 0.734104, -log2l(0.734104) .quad 0x3fe75b37e875b37f, 0x3fdd12938a39d6f0 // 0.729885, -log2l(0.729885) .quad 0x3fe7390d2a6c405e, 0x3fdd9a086f4ad416 // 0.725714, -log2l(0.725714) .quad 0x3fe71745d1745d17, 0x3fde20b7bd4365a8 // 0.721591, -log2l(0.721591) .quad 0x3fe6f5e02e4850ff, 0x3fdea6a3b152b1e6 // 0.717514, -log2l(0.717514) .quad 0x3fe6d4da9b536a6d, 0x3fdf2bce7ef7d06b // 0.713483, -log2l(0.713483) .quad 0x3fe6b4337c6cb157, 0x3fdfb03a50395dba // 0.709497, -log2l(0.709497) .quad 0x3fe693e93e93e93f, 0x3fe019f4a2edc134 // 0.705556, -log2l(0.705556) .quad 0x3fe673fa57b0cbab, 0x3fe05b6ebbca3d9a // 0.701657, -log2l(0.701657) .quad 0x3fe6546546546546, 0x3fe09c8c7a1fd74c // 0.697802, -log2l(0.697802) .quad 0x3fe63528917c80b3, 0x3fe0dd4ee107ae0a // 0.693989, -log2l(0.693989) .quad 0x3fe61642c8590b21, 0x3fe11db6ef5e7873 // 0.690217, -log2l(0.690217) .quad 0x3fe5f7b282135f7b, 0x3fe15dc59fdc06b7 // 0.686486, -log2l(0.686486) .quad 0x3fe5d9765d9765d9, 0x3fe19d7be92a2310 // 0.682796, -log2l(0.682796) .quad 0x3fe5bb8d015e75bc, 0x3fe1dcdabdfad537 // 0.679144, -log2l(0.679144) .quad 0x3fe59df51b3bea36, 0x3fe21be30d1e0ddb // 0.675532, -log2l(0.675532) .quad 0x3fe580ad602b580b, 0x3fe25a95c196bef3 // 0.671958, -log2l(0.671958) .quad 0x3fe563b48c20563b, 0x3fe298f3c2af6595 // 0.668421, -log2l(0.668421) .quad 0x3fe5470961d7ca63, 0x3fe2d6fdf40e09c5 // 0.664921, -log2l(0.664921) .quad 0x3fe52aaaaaaaaaab, 0x3fe314b535c7b89e // 0.661458, -log2l(0.661458) .quad 0x3fe50e97366227cb, 0x3fe3521a64737cf3 // 0.658031, -log2l(0.658031) .quad 0x3fe4f2cddb0d3225, 0x3fe38f2e593cda73 // 0.654639, -log2l(0.654639) .quad 0x3fe4d74d74d74d75, 0x3fe3cbf1e9f5cf2f // 0.651282, -log2l(0.651282) .quad 0x3fe4bc14e5e0a72f, 0x3fe40865e9285f33 // 0.647959, -log2l(0.647959) .quad 0x3fe4a12316176410, 0x3fe4448b2627ade3 // 0.64467, -log2l(0.64467) .quad 0x3fe48676f31219dc, 0x3fe480626d20a876 // 0.641414, -log2l(0.641414) .quad 0x3fe46c0f6feb6ac6, 0x3fe4bbec872a4505 // 0.638191, -log2l(0.638191) .quad 0x3fe451eb851eb852, 0x3fe4f72a3a555958 // 0.635, -log2l(0.635) .quad 0x3fe4380a3065e3fb, 0x3fe5321c49bc0c91 // 0.631841, -log2l(0.631841) .quad 0x3fe41e6a74981447, 0x3fe56cc37590e6c5 // 0.628713, -log2l(0.628713) .quad 0x3fe4050b59897548, 0x3fe5a7207b2d815a // 0.625616, -log2l(0.625616) .quad 0x3fe3ebebebebebec, 0x3fe5e1341520db00 // 0.622549, -log2l(0.622549) .quad 0x3fe3d30b3d30b3d3, 0x3fe61afefb3d5201 // 0.619512, -log2l(0.619512) .quad 0x3fe3ba68636adfb0, 0x3fe65481e2a6477b // 0.616505, -log2l(0.616505) .quad 0x3fe3a2027932b48f, 0x3fe68dbd7ddd6e15 // 0.613527, -log2l(0.613527) .quad 0x3fe389d89d89d89e, 0x3fe6c6b27ccfc698 // 0.610577, -log2l(0.610577) .quad 0x3fe371e9f3c04e64, 0x3fe6ff618ce24cd7 // 0.607656, -log2l(0.607656) .quad 0x3fe35a35a35a35a3, 0x3fe737cb58fe5716 // 0.604762, -log2l(0.604762) .quad 0x3fe342bad7f64b39, 0x3fe76ff0899daa49 // 0.601896, -log2l(0.601896) .quad 0x3fe32b78c13521d0, 0x3fe7a7d1c4d64520 // 0.599057, -log2l(0.599057) .quad 0x3fe3146e92a10d38, 0x3fe7df6fae65e424 // 0.596244, -log2l(0.596244) .quad 0x3fe2fd9b8396ba9e, 0x3fe816cae7bd40b1 // 0.593458, -log2l(0.593458) .quad 0x3fe2e6fecf2e6fed, 0x3fe84de4100b0ce2 // 0.590698, -log2l(0.590698) .quad 0x3fe2d097b425ed09, 0x3fe884bbc446ae3f // 0.587963, -log2l(0.587963) .quad 0x3fe2ba6574cae996, 0x3fe8bb529f3ab8f3 // 0.585253, -log2l(0.585253) .quad 0x3fe2a46756e62a46, 0x3fe8f1a9398f2d58 // 0.582569, -log2l(0.582569) .quad 0x3fe28e9ca3a728ea, 0x3fe927c029d3798a // 0.579909, -log2l(0.579909) .quad 0x3fe27904a7904a79, 0x3fe95d980488409a // 0.577273, -log2l(0.577273) .quad 0x3fe2639eb2639eb2, 0x3fe993315c28e8fb // 0.574661, -log2l(0.574661) .quad 0x3fe24e6a171024e7, 0x3fe9c88cc134f3c3 // 0.572072, -log2l(0.572072) .quad 0x3fe239662b9f91cb, 0x3fe9fdaac2391e1c // 0.569507, -log2l(0.569507) .quad 0x3fe2249249249249, 0x3fea328bebd84e80 // 0.566964, -log2l(0.566964) .quad 0x3fe20fedcba98765, 0x3fea6730c8d44efa // 0.564444, -log2l(0.564444) .quad 0x3fe1fb78121fb781, 0x3fea9b99e21655eb // 0.561947, -log2l(0.561947) .quad 0x3fe1e7307e4ef157, 0x3feacfc7beb75e94 // 0.559471, -log2l(0.559471) .quad 0x3fe1d31674c59d31, 0x3feb03bae40852a0 // 0.557018, -log2l(0.557018) .quad 0x3fe1bf295cc93903, 0x3feb3773d59a05ff // 0.554585, -log2l(0.554585) .quad 0x3fe1ab68a0473c1b, 0x3feb6af315450638 // 0.552174, -log2l(0.552174) .quad 0x3fe197d3abc65f4f, 0x3feb9e3923313e58 // 0.549784, -log2l(0.549784) .quad 0x3fe18469ee58469f, 0x3febd1467ddd70a7 // 0.547414, -log2l(0.547414) .quad 0x3fe1712ad98b8957, 0x3fec041ba2268731 // 0.545064, -log2l(0.545064) .quad 0x3fe15e15e15e15e1, 0x3fec36b90b4ebc3a // 0.542735, -log2l(0.542735) .quad 0x3fe14b2a7c2fee92, 0x3fec691f33049ba0 // 0.540426, -log2l(0.540426) .quad 0x3fe1386822b63cbf, 0x3fec9b4e9169de22 // 0.538136, -log2l(0.538136) .quad 0x3fe125ce4feeb7a1, 0x3feccd479d1a1f94 // 0.535865, -log2l(0.535865) .quad 0x3fe1135c81135c81, 0x3fecff0acb3170e3 // 0.533613, -log2l(0.533613) .quad 0x3fe10112358e75d3, 0x3fed30988f52c6d3 // 0.531381, -log2l(0.531381) .quad 0x3fe0eeeeeeeeeeef, 0x3fed61f15bae4663 // 0.529167, -log2l(0.529167) .quad 0x3fe0dcf230dcf231, 0x3fed9315a1076fa2 // 0.526971, -log2l(0.526971) .quad 0x3fe0cb1b810ecf57, 0x3fedc405cebb27dc // 0.524793, -log2l(0.524793) .quad 0x3fe0b96a673e2808, 0x3fedf4c252c5a3e1 // 0.522634, -log2l(0.522634) .quad 0x3fe0a7de6d1d6086, 0x3fee254b99c83339 // 0.520492, -log2l(0.520492) .quad 0x3fe096771e4d528c, 0x3fee55a20f0eecf9 // 0.518367, -log2l(0.518367) .quad 0x3fe0853408534085, 0x3fee85c61c963f0d // 0.51626, -log2l(0.51626) .quad 0x3fe07414ba8f0741, 0x3feeb5b82b10609b // 0.51417, -log2l(0.51417) .quad 0x3fe06318c6318c63, 0x3feee578a1eaa83f // 0.512097, -log2l(0.512097) .quad 0x3fe0523fbe3367d7, 0x3fef1507e752c6c8 // 0.51004, -log2l(0.51004) .quad 0x3fe04189374bc6a8, 0x3fef4466603be71d // 0.508, -log2l(0.508) .quad 0x3fe030f4c7e7859c, 0x3fef73947063b3fd // 0.505976, -log2l(0.505976) .quad 0x3fe0208208208208, 0x3fefa2927a574422 // 0.503968, -log2l(0.503968) .quad 0x3fe0103091b51f5e, 0x3fefd160df77ed7a // 0.501976, -log2l(0.501976) .quad 0x3fe0000000000000, 0x3ff0000000000000 // 0.5, -log2l(0.5)
// Taylor series coefficients for log2 stagepowf_logTaylor: .double -0.5, 0.33333333333333333333333333333333, -0.25, 0.2
.literal8oneD: .double 1.0d128: .double 128.0dm150: .double -150.0recip_ln2: .quad 0x3ff71547652b82fe // 1.0 / ln(2)
.literal4infF: .long 0x7f800000 // infminfF: .long 0xff800000 // -infoneF: .long 0x3f800000 // 1.0fmoneF: .long 0xbf800000 // -1.0fmaxy: .long 0x4effffff // 0x1.0p31f - 1 ulpminy: .long 0xCeffffff // -0x1.0p31f + 1 ulpmantissaMask: .long 0x007fffff
.text
#if defined( __x86_64__ ) #define SI_P %rsi #define DI_P %rdi #define RELATIVE_ADDR( _a) (_a)( %rip )#else #define SI_P %esi #define DI_P %edi #define RELATIVE_ADDR( _a) (_a)-0b( BX_P )#endif
ENTRY( powf )#if defined( __i386__ ) movl FRAME_SIZE( STACKP ), %eax movl 4+FRAME_SIZE( STACKP ), %edx movss FRAME_SIZE( STACKP ), %xmm0 movss 4+FRAME_SIZE( STACKP ), %xmm1#else movd %xmm0, %eax movd %xmm1, %edx#endif
//early out for x == 1.0 cmpl $0x3f800000, %eax //if( x == 1.0 ) je 6f // goto 6
//early out for y == 1.0 (costs 1 cycle for x86_64, free for i386) cmpl $0x3f800000, %edx //if( y == 1.0 ) je 6f // goto 6
andl $0x7fffffff, %edx // |y| // Find out if y is an integer without raising inexact // Note tested over entire range. Fails for Inf/NaN, but we don't care about that here. push BX_P push SI_P push DI_P
#if defined( __i386__ ) call 0f0: pop BX_P#else xorq %rdi, %rdi#endif
// check to see if we fell into an edge case subl $1, %eax subl $1, %edx cmpl $0x7f7fffff, %eax // if( x < 0 || x == inf || isnan(x) ) jae 7f // goto 7 cmpl $0x4affffff, %edx // if( |y| >= 0x1.0p23 || 0 == y || isnan(y) ) jae 7f // goto 7
cmpl $0x3effffff, %edx // if( |y| == 0.5f ) je 8f // goto 8 // The main part of pow:// 0 < x < inf, |y| < 0x1.0p31, x != 1, y != 0
addl $1, %eax andl $0x7fffffff, %eax // |x|
#if 0 // if y is integer, call ipowf instead addl $1, %edx movl %edx, %edi // |y| lea RELATIVE_ADDR(gMaskShift), CX_P // gMaskShift ptr shrl $23, %edi // |y| >> 23 movzbl (CX_P, DI_P, 1), %ecx // gMaskShift[ |y| >> 23 ] mov $0x3fffffff, DI_P // 0x3fffffff shrl %cl, %edi // 0x3fffffff >> gMaskShift[ |y| >> 23 ] andl %edx, %edi // fractional part of y cmpl $0, %edi je ___ipowf#endif
//separate |x| into 2**i * 1.m movss RELATIVE_ADDR( mantissaMask), %xmm3 movss RELATIVE_ADDR( oneF), %xmm2 andps %xmm3, %xmm0 // m orps %xmm2, %xmm0 // 1.m shrl $23, %eax // exponent + bias cmpl $0, %eax jne 1f // normalize denormal x subss %xmm2, %xmm0 // 1.m - 1.0 movd %xmm0, %eax shrl $23, %eax // exponent + bias andps %xmm3, %xmm0 // m orps %xmm2, %xmm0 // 1.m subl $126, %eax
1: subl $127, %eax // i = exponent - bias cvtsi2sd %eax, %xmm3 // log2x = (double) i //check for unit mantissa ucomiss %xmm2, %xmm0 // if( 1.m == 1.0 ) je 2f // skip to 2
//handle non-unit mantissa here movd %xmm0, %eax // set aside 1.m#if defined( __x86_64__ ) cdqe#endif cvtss2sd %xmm0, %xmm0 // r = (double) 1.m lea RELATIVE_ADDR( powf_log2_c ), CX_P // use the top 7 bits of the mantissa to index the powf_log2_c table shr $(23-7-4), AX_P and $0x7f0, AX_P // reduce r to 1-2**7 < r < 1+2**-7 mulsd (CX_P, AX_P, 1), %xmm0 // r *= powf_log2_c[ 2 * index ]
// compensate in log2x by adding powf_log2_c[ 2 * index + 1] // do this early so that we force -1.0 + 1.0 to avoid (-1.0 + tiny) + 1.0 later. // Precision loss from this is at most 7 bits, which is acceptable addsd 8(CX_P, AX_P, 1), %xmm3 // log2x + powf_log2_c[ 2 * index + 1] // we calculate log2(r) as log1p( r-1 ) / ln(2) subsd RELATIVE_ADDR(oneD), %xmm0 // r -= 1.0
// log(1+r) = r - rr/2 + rrr/3 - rrrr/4 + rrrrr/5 // with -2**-7 < r < 2**-7, should be good to (5+1)*7 +2 = 44 bits of accuracy or so // (5+1) because the error term is roughly equal to the missing r**6/6 term lea RELATIVE_ADDR( powf_logTaylor ), CX_P movsd 8(CX_P), %xmm4 movsd 24(CX_P), %xmm5 movapd %xmm0, %xmm2 // r mulsd %xmm0, %xmm0 // rr mulsd %xmm2, %xmm4 // 0.333333333333333333333r mulsd %xmm2, %xmm5 // 0.2r addsd (CX_P), %xmm4 // -0.5 + 0.333333333333333333333r addsd 16(CX_P), %xmm5 // -0.25 + 0.2r mulsd %xmm0, %xmm4 // -0.5rr + 0.333333333333333333333rrr mulsd %xmm0, %xmm0 // rrrr addsd %xmm2, %xmm4 // r - 0.5rr + 0.333333333333333333333rrr mulsd %xmm0, %xmm5 // -0.25rrrr + 0.2rrrrr addsd %xmm5, %xmm4 // r - 0.5rr + 0.333333333333333333333rrr - 0.25rrrr + 0.2rrrrr mulsd RELATIVE_ADDR( recip_ln2), %xmm4 // ( r - 0.5rr + 0.333333333333333333333rrr - 0.25rrrr + 0.2rrrrr ) * (1/ln(2)) addsd %xmm4, %xmm3 // log2x + powf_log2_c[ 2 * index + 1] + ( r - 0.5rr + 0.333333333333333333333rrr - 0.25rrrr + 0.2rrrrr ) * (1/ln(2))
// multiply by y2: cvtss2sd %xmm1, %xmm0 mulsd %xmm3, %xmm0 // y * log2(x)
ucomisd RELATIVE_ADDR( d128), %xmm0 // if( ylog2(x) >= 128 ) jae 4f // goto 4
ucomisd RELATIVE_ADDR( dm150), %xmm0 // if( ylog2(x) <= -150 jbe 4f // goto 4
// separate ylog2(x) into i + f cvttsd2si %xmm0, %eax // i = (int) ylog2(x) cvtsi2sd %eax, %xmm1 // trunc( ylog2(x) ) subsd %xmm1, %xmm0 // f
// calculate 2**i addl $1023, %eax // exponent + bias movd %eax, %xmm7 // move to vector register psllq $52, %xmm7 // shift exponent + bias into place
// early out for power of 2 xorpd %xmm6, %xmm6 ucomisd %xmm0, %xmm6 movsd RELATIVE_ADDR( oneD), %xmm1 je 3f //f = exp2(f)#if defined( __SSE3__ ) movddup %xmm0, %xmm1 // { f, f }#else movapd %xmm0, %xmm1 unpcklpd %xmm1, %xmm1 // { f, f }#endif mulsd %xmm0, %xmm0 // ff = f*f movapd %xmm1, %xmm3 // { f, f } lea RELATIVE_ADDR( powf_exp_c ), CX_P mulpd 48(CX_P), %xmm1 // { c3f, (c7/c8)f } mulpd 16(CX_P), %xmm3 // { c1f, (c5/c8)f }#if defined( __SSE3__ ) movddup %xmm0, %xmm4 // { ff, ff }#else movapd %xmm0, %xmm4 unpcklpd %xmm4, %xmm4 // { ff, ff }#endif mulsd %xmm0, %xmm0 // ffff = ff * ff addpd 32(CX_P), %xmm1 // { c2 + c3f, (c6/c8) + (c7/c8)f } addpd (CX_P), %xmm3 // { c0 + c1f, (c4/c8) + (c5/c8)f } mulpd %xmm4, %xmm1 // { c2ff + c3fff, (c6/c8)ff + (c7/c8)fff } addsd %xmm0, %xmm3 // { c0 + c1x, (c4/c8) + (c5/c8)f + ffff } mulsd 64(CX_P), %xmm0 // c8ffff addpd %xmm1, %xmm3 // { c0 + c1f + c2ff + c3fff, (c4/c8) + (c5/c8)f + (c6/c8)ff + (c7/c8)fff + ffff } movhlps %xmm3, %xmm1 // { ?, c0 + c1f + c2ff + c3fff } mulsd %xmm0, %xmm3 // { ..., c8ffff* ((c4/c8) + (c5/c8)f + (c6/c8)ff + (c7/c8)fff + ffff) } addsd %xmm3, %xmm1 // c0 + c1f + c2ff + c3fff + c4ffff + c5fffff + c6ffffff + c7fffffff + c8fffffffff
// scale by 2**i, and convert to float3: mulsd %xmm1, %xmm7 xorps %xmm0, %xmm0 cvtsd2ss %xmm7, %xmm0
pop DI_P pop SI_P pop BX_P#if defined( __i386__ ) movss %xmm0, FRAME_SIZE( STACKP ) flds FRAME_SIZE( STACKP )#endif ret
// overflow / underflow4: xorpd %xmm1, %xmm1 // 0 cmpltsd %xmm0, %xmm1 // 0 < ylog2(x) ? -1LL : 0 movd %xmm1, %eax // 0 < ylog2(x) ? -1U : 0 andl $0xfff, %eax // 0 < ylog2(x) ? 0xfff : 0 xorl $0x801, %eax // 0 < ylog2(x) ? 0x7fe : 0x801 movd %eax, %xmm2 // 0 < ylog2(x) ? 0x7fe : 0x801 psllq $52, %xmm2 // 0 < ylog2(x) ? 0x1.0p+1023 : -0x1.0p-1022 mulsd %xmm0, %xmm2 // result = ylog2(x) * (0 < ylog2(x) ? 0x1.0p+1023 : -0x1.0p-1022) xorps %xmm0, %xmm0 // 0 cvtsd2ss %xmm2, %xmm0 // convert result to float jmp 9f
// ( x < 0 && isfinite(x) && |y| is not in { 0, inf, NaN } ) or x is unknown, but |y| >= 0x1.0p235: cmpl $0, %edi // if( y is not an integer ) jne 8f // since we know y is an integer, we can just call ipowf jmp ___ipowf 6: // x == 1.0f return x#if defined( __i386__ ) flds FRAME_SIZE( STACKP )#endif ret // A whole basket of special cases lands here // (x <= 0 || x == Inf || isnan(x)) or ( |y| >= 0x1.0p23f || y == 0 || isnan(y) ) // all we have to do is figure out which one!7: addl $1, %eax // |y| addl $1, %edx // |y| andl $0x7fffffff, %eax // |x| cmpl $0, %edx // if( |y| == 0 ) je 4f // goto 4
// (x <= 0 || x == Inf || isnan(x)) or ( |y| >= 0x1.0p23 || isnan(y) )
//check for NaNs ucomiss %xmm0, %xmm1 jp 7f
// (x <= 0 || x == Inf ) or |y| >= 0x1.0p23f // calculate fractional part of y and ones bit of y movl %edx, %edi // |y| lea RELATIVE_ADDR(gMaskShift), CX_P // gMaskShift ptr shrl $23, %edi // |y| >> 23 movzbl (CX_P, DI_P, 1), %ecx // gMaskShift[ |y| >> 23 ] mov $0x3fffffff, DI_P // 0x3fffffff mov $0x40000000, SI_P // 0x40000000 shrl %cl, %edi // 0x3fffffff >> gMaskShift[ |y| >> 23 ] shrl %cl, %esi // 0x40000000 >> gMaskShift[ |y| >> 23 ] andl %edx, %edi // fractional part of y andl %edx, %esi // ones bit of y
// if( x == 0 ) goto 2 xorps %xmm2, %xmm2 ucomiss %xmm0, %xmm2 je 2f // (x < 0 || x == Inf ) or |y| >= 0x1.0p23f
// if( |y| == inf ) goto 3 cmpl $0x7f800000, %edx je 3f
// (x < 0 || x == Inf) or ( 0x1.0p23f <= |y| < inf )
// if( x == inf ) goto 5 ucomiss RELATIVE_ADDR( infF ), %xmm0 je 5f
// x < 0 or ( 0x1.0p23f <= |y| < inf )
// negative finite x or large y go off to be considered for ipowf ucomiss RELATIVE_ADDR( minfF ), %xmm0 // if( x != -inf ) ja 5b // goto the other 5
// x == -inf
// At this point, we know that x is -Inf and |y| is not in { 0, Inf, NaN }. // Deal with y is odd integer cases // if( 0 == fractionalBits && 0 != onesBit ) movl %edi, %ecx // fractional Bits subl $1, %ecx // fractionalBits == 0 ? -1 : some non-negative number sarl $31, %ecx // fractionalBits == 0 ? -1 : 0 andl %esi, %ecx // fractionalBits == 0 ? onesBit : 0 cmpl $0, %ecx // if( 0 == fractionalBits && 0 != onesBit ) jne 6f // goto 6 // x = -inf, |y| is not in { 0, Inf, NaN }, and y is not an odd integer
// if( 0.0f < y ) return -x; else return 0 cmpltss %xmm1, %xmm2 // 0.0f < y ? -1 : 0 andps %xmm2, %xmm0 // 0.0f < y ? x : 0 pslld $31, %xmm2 // 0.0f < y ? 0x80000000 : 0 xorps %xmm2, %xmm0 // 0.0f < y ? -x : 0 jmp 9f //return 0
// x < 0 && y is not an integer, or |y| == 0.5f 8: xorps %xmm2, %xmm2 cmpless %xmm0, %xmm2 // 0 <= x ? -1 : 0 andps %xmm1, %xmm2 // 0 <= x ? y : 0 xorps %xmm3, %xmm3 ucomiss %xmm2, %xmm3 // if( x >= 0 && 0 > y ) ja 1f // goto 1
sqrtss %xmm0, %xmm0 jmp 9f //return // y == -0.5f && x > 01: cvtss2sd %xmm0, %xmm0 movsd RELATIVE_ADDR( oneD ), %xmm1 divsd %xmm0, %xmm1 sqrtsd %xmm1, %xmm1 xorps %xmm0, %xmm0 cvtsd2ss %xmm1, %xmm0 jmp 9f // return
// x == 02: // if( y is an odd integer ) goto 8 movl %edi, %ecx // fractional Bits subl $1, %ecx // fractionalBits == 0 ? -1 : some non-negative number sarl $31, %ecx // fractionalBits == 0 ? -1 : 0 andl %esi, %ecx // fractionalBits == 0 ? onesBit : 0 cmpl $0, %ecx // if( fractionalBits == 0 && 0 != onesBit ) jne 8f // y is an odd integer, goto 8
xorps %xmm0, %xmm0 // x = fabsf(x) ucomiss %xmm1, %xmm0 // if( 0 < y ) jb 9f // return x
//return 1.0 / f movss RELATIVE_ADDR( oneF ), %xmm1 divss %xmm0, %xmm1 // return inf and set div/0 movaps %xmm1, %xmm0 jmp 9f
// |y| == inf3: ucomiss RELATIVE_ADDR( moneF ), %xmm0 // if( -1.0f == x ) je 4f // return 1.0f
cmpl $0x3f7fffff, %eax // if( |x| > 1.0f ) ja 1f // goto 1f xorps %xmm0, %xmm0 // 0.0f cmpnless %xmm1, %xmm0 // y == inf ? 0 : -1 psrld $1, %xmm0 // y == inf ? 0 : 0x7fffffff andps %xmm1, %xmm0 // y == inf ? 0 : inf jmp 9f // return // return 1.0f4: movl $1, %ecx xorps %xmm0, %xmm0 cvtsi2ss %ecx, %xmm0 jmp 9f
// x == inf5: xorps %xmm2, %xmm2 cmpltss %xmm1, %xmm2 // 0 < y ? -1 : 0 andps %xmm2, %xmm0 // 0 < y ? x : 0 jmp 9f
// 0 == fractionalBits && 0 != onesBit 6: xorps %xmm2, %xmm2 ucomiss %xmm1, %xmm2 // if( 0 < y ) jb 9f //return x
movl $0x80000000, %ecx // -0.0f movd %ecx, %xmm0 // copy to xmm, zero high part of register jmp 9f //return -0.0
7: // NaN addss %xmm1, %xmm0 jmp 9f
// x == 0, y is an odd integer8: ucomiss %xmm1, %xmm2 // if( 0 < y ) jb 9f // return x //return 1.0 / f movss RELATIVE_ADDR( oneF ), %xmm1 divss %xmm0, %xmm1 // return inf and set div/0 movaps %xmm1, %xmm0 jmp 9f
//|y| == inf, |x| > 1.0f1: xorps %xmm0, %xmm0 cmpltss %xmm1, %xmm0 andps %xmm1, %xmm0 jmp 9f
.align 4 // return value in %xmm09: pop DI_P pop SI_P pop BX_P#if defined( __i386__ ) movss %xmm0, FRAME_SIZE( STACKP ) flds FRAME_SIZE( STACKP )#endif ret // x and y passed in in xmm0 and xmm1// result returned in xmm0 // BX_P already points to label 0 above___ipowf: // clamp INT_MIN <= y < INT_MAX. Values outside this range can't be odd numbers. maxss RELATIVE_ADDR( miny ), %xmm1 minss RELATIVE_ADDR( maxy ), %xmm1 cvttss2si %xmm1, %edx // (int) y cvtss2sd %xmm0, %xmm0 // x movsd RELATIVE_ADDR( oneD ), %xmm2 // r = 1.0
cmpl $0, %edx // if( y >= 0 ) jge 1f // goto 4 // y < 0 movapd %xmm0, %xmm1 // x movapd %xmm2, %xmm0 // 1.0 divsd %xmm1, %xmm0 // 1.0 / x negl %edx 1: test $1, %edx jz 3f // if( |y| is odd ) movapd %xmm0, %xmm2 // r = x jmp 3f
.align 4// do{2: mulsd %xmm0, %xmm0 // x *= x test $1, %edx jz 3f // if( |y| is odd ) continue mulsd %xmm0, %xmm2 // r *= x
3: shrl $1, %edx // |y| >>= 1 test $-1, %edx jnz 2b // if( y ) continue// }while( |y| )
// round to float xorps %xmm0, %xmm0 // 0 cvtsd2ss %xmm2, %xmm0 // (float) r
//exit jmp 9b