#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
static const double kBig = 6.755399441055744e+15;
static const double kPiBy2 = 1.570796326794896619231322; static const double kPiBy2Tail1 = 6.123233995736766036e-17; static const double kPiBy2Tail1m = 6.123233995736764803e-17; static const double kPiBy2Tail2 = -1.4973849048591698e-33; static const double kRecip = 6.869637070016756334588131e-12; static const double kLimit1 = 0.0859375; static const double kLimit2 = 1.015625;
static const double kLimit3 = 13.0;
static const double kLimit4 = 8.11296384146066817e+31; static const Float kInfinityu = {0x7f800000};
static const float factor[12] = {
32.0,16.0,8.0,4.0,4.0,2.0,2.0,1.0,1.0,1.0,1.0,1.0
};
static const float adjust[12] = {
32.0,64.0,88.0,104.0,104.0,116.0,116.0,124.0,124.0,124.0,124.0,124.0
};
static const double atancoeff[15] = {
145568097675.0, -48522699225.0,
29113619535.0, -20795442525.0, 16174233075.0,
-13233463425.0, 11197545975.0, -9704539845.0,
8562829275.0, -7661478825.0, 6931814175.0,
-6329047725.0, 5822723907.0, -5391411025.0,
5019589575.0
};
static const uint32_t atanTable[] = {
0x3FB80000, 0x96AC4A5A, 0x3FB7EE18, 0xBB5F2BDC, 0x3B1F7151, 0xA60EEDFD,
0x3FBC0000, 0xAB2E5990, 0x3FBBE39F, 0x679755D8, 0x3B5AF7C1, 0xBBB2A139,
0x3FC00000, 0x7F748304, 0x3FBFD5BB, 0x95A94110, 0x3B5B0A0C, 0x811C880E,
0x3FC20000, 0x3DF2157F, 0x3FC1E1FB, 0x37C2C7E5, 0x3B617207, 0xB776AF4F,
0x3FC40000, 0x30754446, 0x3FC3D6EF, 0x1814018D, 0xBB63A092, 0x88503B16,
0x3FC60000, 0x7996ECD4, 0x3FC5C981, 0x94588C25, 0xBB5B6979, 0x558D8B64,
0x3FC80000, 0x76D31B22, 0x3FC7B97B, 0xBE985C07, 0x3B6708AD, 0xCC1C4519,
0x3FCA0000, 0xEADCE47A, 0x3FC9A6A9, 0xCAFAF9A8, 0x3B62971C, 0xA562CA82,
0x3FCC0000, 0x7E8C3F27, 0x3FCB90D7, 0xCB57348B, 0xBB493C17, 0xA875159C,
0x3FCE0000, 0x5B063A6A, 0x3FCD77D6, 0x3569311A, 0xBB5A33AB, 0xA192F16A,
0x3FD00000, 0x64EF22FF, 0x3FCF5B76, 0xB72AE095, 0x3B52E368, 0xB50909D6,
0x3FD10000, 0x87DE9AAA, 0x3FD09DC6, 0x16C29778, 0xBB55B0DC, 0xF3418D68,
0x3FD20000, 0x0AAD7C84, 0x3FD18BF5, 0xACF10E71, 0x3B41EAC1, 0x346154B4,
0x3FD30000, 0x1D92088F, 0x3FD27837, 0x3B84D2FF, 0x3B48BBE9, 0xDC9998E5,
0x3FD40000, 0xCC4E6285, 0x3FD36277, 0xF12910F6, 0xBB6585F3, 0xFA8B319B,
0x3FD50000, 0x766A7592, 0x3FD44AA4, 0xA1AA8D78, 0x3B6C6289, 0x13F656CF,
0x3FD60000, 0xCFA33544, 0x3FD530AE, 0x5303BB5D, 0x3B63F533, 0xC1117E20,
0x3FD70000, 0xC925EED6, 0x3FD61484, 0xB52DF321, 0xBB70059A, 0x66E93FA9,
0x3FD80000, 0x54AA10D9, 0x3FD6F619, 0x8C1ECA84, 0x3B60C60C, 0x08BED883,
0x3FD90000, 0x113AC5B7, 0x3FD7D560, 0x5A568BA4, 0xBB541817, 0xB5654BFC,
0x3FDA0000, 0x6D101F2E, 0x3FD8B24D, 0x96E71249, 0xBB55A721, 0x1379C827,
0x3FDB0000, 0xE246C017, 0x3FD98CD6, 0x05641F92, 0x3B6918FF, 0xCA4E1DF9,
0x3FDC0000, 0x39B58346, 0x3FDA64EE, 0xF43C341A, 0x3B6633D0, 0xEECC747A,
0x3FDD0000, 0x4C36F71C, 0x3FDB3A91, 0x5CE1B431, 0xBB61E754, 0x0B1A4FB4,
0x3FDE0000, 0x66CA921B, 0x3FDC0DB5, 0x1D94F195, 0x3B571DE1, 0xE6981964,
0x3FDF0000, 0x3AD3639E, 0x3FDCDE53, 0x72D1AC94, 0xBB736718, 0x370B3962,
0x3FE00000, 0x7815CCCD, 0x3FDDAC67, 0xC5849B77, 0x3B52FF7C, 0xC7D84489,
0x3FE08000, 0xF40DDDFC, 0x3FDE77ED, 0x00AEBFCD, 0x3B62E958, 0xE33FF707,
0x3FE10000, 0xA6C93DAD, 0x3FDF40DE, 0x0F7AA915, 0x3B77641D, 0x15CA0475,
0x3FE18000, 0x85EBF4D7, 0x3FE0039C, 0xDAD8BEF9, 0x3B722996, 0x8391BF09,
0x3FE20000, 0xC129DA76, 0x3FE0657F, 0x27977717, 0xBB761869, 0xC045B9BD,
0x3FE28000, 0x93F751BC, 0x3FE0C614, 0xCA41B0AE, 0xBB7A8A92, 0x3C961154,
0x3FE30000, 0x7C1EF6CE, 0x3FE1255D, 0xF7C0A55C, 0x3B591CBF, 0x478B1C38,
0x3FE38000, 0x169FA349, 0x3FE1835A, 0x993DD546, 0x3B6FBFAB, 0x698175A0,
0x3FE40000, 0x27D6B511, 0x3FE1E00B, 0xC884E585, 0xBB248C8D, 0x9BB39A1C,
0x3FE48000, 0x6BE89876, 0x3FE23B72, 0x2F4EF788, 0xBB815690, 0x7DF099CF,
0x3FE50000, 0x280F4ABE, 0x3FE2958E, 0x7530C25F, 0xBB72B7E0, 0xEE07EAE7,
0x3FE58000, 0x2FEDFD47, 0x3FE2EE62, 0xA50C9A07, 0xBB2B522D, 0xBFFA3A8C,
0x3FE60000, 0x875C63BA, 0x3FE345F0, 0x78B8C001, 0xBB408ED6, 0xD056E4C4,
0x3FE68000, 0x8A368DBA, 0x3FE39C39, 0x79511750, 0x3B77CB64, 0x0C777681,
0x3FE70000, 0xEE7CA285, 0x3FE3F140, 0x55DECB32, 0xBB5AA9DC, 0xBC68CF79,
0x3FE78000, 0xA4904A32, 0x3FE44506, 0xC661B510, 0x3B7EBCE3, 0x3EB6C0CD,
0x3FE80000, 0x4CA9F6E2, 0x3FE4978F, 0xD4373CAA, 0x3AFF9675, 0xAB73F345,
0x3FE88000, 0xD7B95553, 0x3FE4E8DE, 0xE3B770C7, 0xBB715BD3, 0xBFA6C0AD,
0x3FE90000, 0x849C7E04, 0x3FE538F5, 0xCDE27023, 0xBB6B45C2, 0x3E5A5CC8,
0x3FE98000, 0xC171C662, 0x3FE587D8, 0x95C38C2E, 0x3B7887A1, 0x345B0752,
0x3FEA0000, 0x4321104C, 0x3FE5D589, 0xAF86142F, 0xBB8438EB, 0x1AE94AC1,
0x3FEA8000, 0x838D8D26, 0x3FE6220D, 0x5F66C783, 0x3B7CA422, 0x8F96E72F,
0x3FEB0000, 0xE4DA897E, 0x3FE66D66, 0xBED2B224, 0xBB35FB4D, 0xFD76AD75,
0x3FEB8000, 0xD775765C, 0x3FE6B799, 0x0DF9D017, 0x3B7B5383, 0xE0D80102,
0x3FEC0000, 0xAA053CA7, 0x3FE700A8, 0x25C3BB63, 0xBB7D5234, 0x29157DFC,
0x3FEC8000, 0x773AAC9B, 0x3FE74897, 0xD237C658, 0x3B6ACFB2, 0xF4E46A15,
0x3FED0000, 0x3FA5ABCF, 0x3FE78F6B, 0xE04F6C27, 0xBB8680F4, 0x55D1C7DF,
0x3FED8000, 0xD31E46AA, 0x3FE7D528, 0x9AC0235F, 0xBB61A1B0, 0x1EABC199,
0x3FEE0000, 0xBB676A41, 0x3FE819D1, 0x1AD33A87, 0xBB7726EA, 0x735248A8,
0x3FEE8000, 0x76CB5B54, 0x3FE85D69, 0x95ABE3CC, 0x3B772D9C, 0x044A3CEB,
0x3FEF0000, 0x26506206, 0x3FE89FF6, 0x131BC92E, 0x3B860889, 0x51E5F051,
0x3FEF8000, 0xED69BF50, 0x3FE8E17B, 0x2230274E, 0xBB7FB7C5, 0x5CA2A11A,
0x3FF00000, 0x38A5805B, 0x3FE921FB, 0x8CE9AD0F, 0xBB728518, 0x56659FDA,
0x3FF08000, 0xE4CDB94C, 0x3FE9A001, 0x86F9991F, 0xBB7F4F21, 0xC242FD58,
0x3FF10000, 0x29DA60EA, 0x3FEA1A26, 0x1A19C30A, 0xBB87F5DE, 0xE73B7B3A,
0x3FF18000, 0x9BBB03C5, 0x3FEA908B, 0x882B14EE, 0x3B8A364F, 0xC319894C,
0x3FF20000, 0x42A338B5, 0x3FEB034F, 0x7337C8D1, 0xBB768B38, 0xD9BDAC43,
0x3FF28000, 0xED0428FC, 0x3FEB7292, 0x7FBAC8B6, 0x3B7B400A, 0x53FD9AAF,
0x3FF30000, 0x36946784, 0x3FEBDE71, 0x1A8E3A61, 0x3B6C2B37, 0x13D123E1,
0x3FF38000, 0x21F16647, 0x3FEC470A, 0xDA7DBAA7, 0xBB523353, 0xD7A5886A,
0x3FF40000, 0xF06227B4, 0x3FECAC7D, 0x132231D2, 0xBB8AEFE3, 0x0558B34B,
0x3FF48000, 0xF9AB1CF7, 0x3FED0EE2, 0xE23FB012, 0xBB4465DC, 0xD4A6C05D,
0x3FF50000, 0x1DAE2F90, 0x3FED6E57, 0xE51C7E16, 0xBB85A9B9, 0x83503D52,
0x3FF58000, 0xD679A51F, 0x3FEDCAF8, 0xC6A4C801, 0x3B8576BE, 0xCEAEFC33,
0x3FF60000, 0xCE367C5E, 0x3FEE24DD, 0xD375A188, 0xBB78F5B1, 0x1E4802A0,
0x3FF68000, 0xCD05D052, 0x3FEE7C20, 0xCBEB8A43, 0xBB8AD216, 0x43D5F156,
0x3FF70000, 0x98AA0697, 0x3FEED0D9, 0xE022B144, 0x3B267908, 0xB6C976B5,
0x3FF78000, 0xA3FB1493, 0x3FEF2320, 0xDB8F1029, 0xBB85B081, 0x83128CFA,
0x3FF80000, 0x681E0063, 0x3FEF730C, 0x12946C3F, 0xBB8AA69F, 0xBF359098,
0x3FF88000, 0x75E901B6, 0x3FEFC0B1, 0xB86DE15F, 0x3B8EAA0C, 0x4663839A,
0x3FF90000, 0x736A4590, 0x3FF00613, 0x4FBE5B7E, 0xBB4E3551, 0x42306F18,
0x3FF98000, 0xC5B189FC, 0x3FF02ABF, 0xA107BE50, 0x3B63FD1C, 0xCE1D6AEB,
0x3FFA0000, 0x41411D33, 0x3FF04E67, 0x3966892F, 0xBB86ED08, 0x00B1EF7B,
0x3FFA8000, 0x4DDFE3DA, 0x3FF07113, 0xDB772C28, 0x3B746CEB, 0x807EA687,
0x3FFB0000, 0xA7AB4BE2, 0x3FF092CE, 0x72AC060E, 0xBB72F0AF, 0x054AD94E,
0x3FFB8000, 0x35EE2E87, 0x3FF0B39F, 0x5C6DBFD7, 0xBB8A699D, 0x14CF8D89,
0x3FFC0000, 0x79493CCE, 0x3FF0D38F, 0x4A3683E2, 0xBB7C178E, 0x4E0FB167,
0x3FFC8000, 0x7C2C16CD, 0x3FF0F2A5, 0xF7C1C52E, 0xBB77BBEC, 0x0F5FF2D9,
0x3FFD0000, 0xF7925B8C, 0x3FF110EB, 0x3A456E0E, 0x3B86B850, 0xB94665FD,
0x3FFD8000, 0xD4EDAFAA, 0x3FF12E66, 0x2A98E0EA, 0x3B8FC04D, 0xF88FF47B,
0x3FFE0000, 0x1822F5F9, 0x3FF14B1D, 0xDB5166C5, 0x3B400FEE, 0xA21D9DDA,
0x3FFE8000, 0x5F375D4F, 0x3FF16719, 0x6EAC8C79, 0xBB45D7BA, 0x99C8508E,
0x3FFF0000, 0x983AF36D, 0x3FF1825F, 0x2745DBAE, 0xBB8849C2, 0x328CE199,
0x3FFF8000, 0xDF7044EE, 0x3FF19CF5, 0x48D90CBC, 0xBB876AC4, 0x23C76E06,
0x40000000, 0xC89A22AD, 0x3FF1B6E1, 0xE3296298, 0xBB83E57A, 0x4E33C68C,
0x40008000, 0xAE85BCDD, 0x3FF1E8D4, 0xB635433E, 0x3B8FDE86, 0xD38237A1,
0x40010000, 0x1E66B4AD, 0x3FF21862, 0xFF00ECE6, 0x3B74E47A, 0x940D38B0,
0x40018000, 0x257005B2, 0x3FF245B5, 0x07E38198, 0xBB90341D, 0xF7AEB5A2,
0x40020000, 0x55656B58, 0x3FF270EF, 0x71D13D73, 0x3B88C680, 0x13906BD2,
0x40028000, 0x466225D9, 0x3FF29A34, 0x0FA8F5E6, 0x3B91F854, 0xF983041D,
0x40030000, 0x717EC3D3, 0x3FF2C1A2, 0x640509A5, 0x3B8E1C75, 0xAF6F1BFF,
0x40038000, 0x75A7B542, 0x3FF2E757, 0x4A697F72, 0x3B92AB3B, 0x412F53FB,
0x40040000, 0x00512EF1, 0x3FF30B6D, 0x7980B2E2, 0x3B8AD702, 0x40329E0C,
0x40048000, 0x09ED9B8A, 0x3FF32DFE, 0x0460EC6D, 0x3B6B92BF, 0x54B284E9,
0x40050000, 0x1B9F04CA, 0x3FF34F1F, 0xC21A2D1B, 0xBB82C1A6, 0xD9CE49BA,
0x40058000, 0x697C1A2B, 0x3FF36EE8, 0x0C4A7DD1, 0x3B73DD4E, 0x5C021394,
0x40060000, 0xABB8A398, 0x3FF38D6A, 0x94FD5FEE, 0x3B373349, 0x3B4FC3A6,
0x40068000, 0x18374FED, 0x3FF3AAB9, 0x8BB1767A, 0x3B72834E, 0xCBC4CA13,
0x40070000, 0xBCD94FDA, 0x3FF3C6E6, 0x7976E8D9, 0xBB931322, 0x7228DCD0,
0x40078000, 0x7BA476B6, 0x3FF3E200, 0xC84E84B4, 0xBB87D21A, 0xC44FB391,
0x40080000, 0xD7180754, 0x3FF3FC17, 0x967F5249, 0xBB64B771, 0x5FA0957B,
0x40090000, 0x6DFDCD42, 0x3FF42D70, 0x558EAD29, 0x3B84B923, 0xDD410F01,
0x400A0000, 0xE7C427AE, 0x3FF45B54, 0xAB8A2C51, 0x3B91D694, 0xE38B7BFF,
0x400B0000, 0x10E3223E, 0x3FF4861B, 0x4FB5B5B4, 0xBB7CB027, 0x50BA194A,
0x400C0000, 0x8EA13F5E, 0x3FF4AE11, 0x11ECF83F, 0x3B8B6616, 0xD89A7B16,
0x400D0000, 0x06D105F7, 0x3FF4D378, 0xC2906964, 0x3B93659C, 0x0737481A,
0x400E0000, 0x727DD5BC, 0x3FF4F68D, 0xF99AE908, 0x3B8514AE, 0x4636F80B,
0x400F0000, 0x3FDFD1C5, 0x3FF51785, 0x020F4064, 0xBB65DB78, 0x8577ADC5,
0x40100000, 0xCF49E240, 0x3FF5368C, 0xC5E4B1C8, 0xBB5E5499, 0xA9C4BF51,
0x40110000, 0x94E90733, 0x3FF56F6F, 0x52E31047, 0x3B8F3B21, 0x724AC4D7,
0x40120000, 0x46578102, 0x3FF5A250, 0x5F4EF405, 0x3B8657A1, 0x1EBE66F0,
0x40130000, 0xCA3D6D0F, 0x3FF5D013, 0xE66FF8CA, 0x3B60B423, 0x26CAA3CA,
0x40140000, 0x32964FAE, 0x3FF5F973, 0x1CEDA34B, 0x3B635D66, 0xC3ED1C18,
0x40150000, 0xA6AA4915, 0x3FF61F06, 0xDE005085, 0x3B6E8C7F, 0x000B0CDD,
0x40160000, 0x8906EAE3, 0x3FF6414D, 0x5593660A, 0x3B95FADD, 0xB8BF7A98,
0x40170000, 0x831FA60C, 0x3FF660B0, 0x3BD94D47, 0x3B96491D, 0x5B9C83F6,
0x40180000, 0x8DCF9463, 0x3FF67D88, 0x73114F7D, 0x3B928119, 0xC1230051,
0x401A0000, 0xCBB9EEAD, 0x3FF6B0BA, 0xFB083C0D, 0x3B93280F, 0xAE0D69E3,
0x401C0000, 0xDEF559E1, 0x3FF6DCC5, 0x8D8B9598, 0x3B86722D, 0x1022C04A,
0x401E0000, 0x744BF1EC, 0x3FF7030D, 0x01605442, 0xBB736E8B, 0xF92A87CA,
0x40200000, 0x456D4528, 0x3FF7249F, 0xB324E4B7, 0x3B675690, 0xDD79F342,
0x40220000, 0x43E2466D, 0x3FF75CBA, 0xD9437CC1, 0x3B46FCCD, 0xEB45ACCD,
0x40240000, 0x18A681EA, 0x3FF789BD, 0x2E09D7EA, 0x3B6C324D, 0xB28DD9D9,
0x40260000, 0xF3FE273A, 0x3FF7AEA3, 0x9C1AAC21, 0x3B8D6F66, 0xFC2D1461,
0x40280000, 0x455A6241, 0x3FF7CD6F, 0x71992B07, 0x3B8B2FE8, 0xBAEBACBC,
0x40280000, 0x455A6241, 0x3FF7CD6F, 0x71992B07, 0x3B8B2FE8, 0xBAEBACBC
};
struct atantblEntry {
double One;
double Two;
double Three;
} atantblEntry;
long double atanl( long double x )
{
double p1, p2, p3, p2a, p4, p5, p6, p7, p8, p29, p30, p33, p35;
double quot, quot2;
double rarg, result, restemp;
double top, top2, topextra, top3;
double bot, bot2;
double r, rarg2, rarg2a;
double temps;
double remainder, residual;
double partial;
double arg, argl, argsq, argsq2;
double temp,temp2;
double sum, sumt, suml;
double absx;
double xHead, xTail;
double res, reslow, resmid, restiny, reshi, resbot;
double p, q;
double prod, prodlow;
Double DeN ;
double fenv;
DblDbl ddx;
int i, rflag;
uint32_t index;
struct atantblEntry *atanPtr = (struct atantblEntry *)atanTable - 6;
FEGETENVD(fenv);
FESETENVD(0.0);
ddx.f = x;
xHead = ddx.d[0];
xTail = ddx.d[1];
if ((xHead != xHead)||(xHead == 0.0)) { FESETENVD(fenv);
return x;
}
absx = fabs(xHead);
if (absx >= kLimit4) { ddx.d[0] = kPiBy2;
ddx.d[1] = kPiBy2Tail1m; if (xHead < 0.0) {
ddx.d[0] = -ddx.d[0];
ddx.d[1] = -ddx.d[1];
}
FESETENVD(fenv);
return (ddx.f);
}
rflag = 0;
if (absx >= kLimit3) { rarg = 1.0/xHead;
rflag = 1; p1 = __FNMSUB( xHead, rarg, 1.0 ); p2 = xTail*rarg; p3 = p1 - p2; rarg2a = p3*rarg;
rarg2 = rarg2a + __FNMSUB( xHead, rarg2a, p3 )*rarg; xHead = rarg;
xTail = rarg2;
}
if( (absx < kLimit1) || (absx >= kLimit3) ) { temp = 2.0*xHead*xTail; argsq = __FMADD( xHead, xHead, temp ); argsq2 = __FMADD( xTail, xTail, __FMSUB( xHead, xHead, argsq ) + temp ); sum = atancoeff[14];
sum = __FMADD( argsq, sum, atancoeff[13] ); sum = __FMADD( argsq, sum, atancoeff[12] ); sum = __FMADD( argsq, sum, atancoeff[11] ); sumt = __FMADD( argsq, sum, atancoeff[10] ); sum = __FMADD( argsq, sum, atancoeff[9] ); suml = __FMADD( argsq, sumt, atancoeff[9] - sum );
for (i = 8; i >= 1; i--) {
temps = __FMADD( argsq, sum, atancoeff[i] );
suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
sum=temps;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, xHead, prod*xTail ); temp = __FMADD( prod, xHead, temp2 ); temp2 = __FMSUB( prod, xHead, temp ) + temp2; sum = temp*kRecip; remainder = __FNMSUB( sum, atancoeff[0], temp ); partial = __FADD( remainder, temp2 );
residual = remainder - partial + temp2;
suml = __FMADD( partial, kRecip, (residual*kRecip) ); res = __FADD( xHead, sum ); reslow = (xHead - res + sum); resmid = __FADD( xTail, suml );
restiny = xTail - resmid + suml;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = res + p;
resbot = (res - reshi + p) + (q + restiny);
if (rflag == 1) { if (xHead > 0) { p1 = kPiBy2;
p2 = kPiBy2Tail1;
p2a = kPiBy2Tail2;
}
else {
p1 = -kPiBy2;
p2 = -kPiBy2Tail1;
p2a = -kPiBy2Tail2;
}
p3 = __FSUB( p1, reshi ); p4 = p1 - p3 - reshi;
p5 = __FSUB( p2, resbot );
p6 = p2 - p5 - resbot + p2a;
p7 = __FADD( p4, p5 );
reshi = __FADD( p3, p7 );
if (fabs(p4) > fabs(p5))
p8 = p4 - p7 + p5;
else
p8 = p5 - p7 + p4;
resbot = p3 - reshi + p7 + (p6 + p8 - restiny);
}
ddx.d[0] = __FADD( reshi, resbot );
ddx.d[1] = (reshi - ddx.d[0]) + resbot;
FESETENVD(fenv);
return (ddx.f);
}
else {
DeN.f = __FMADD( 64.0, absx, kBig );
if (absx > kLimit2) { FESETENVD(kBig+1.0); DeN.f = absx + kBig; FESETENVD(kBig+1.0); DeN.f = absx + kBig; FESETENVD(kBig); index = DeN.i[1];
i = index - 1;
DeN.f = __FMADD( factor[i], absx, kBig + adjust[i] ); }
index = DeN.i[1];
if (xHead < 0) { arg = -xHead;
argl = -xTail;
}
else {
arg = xHead;
argl = xTail;
}
top = __FSUB( arg, atanPtr[(int32_t)index].One );
top2 = __FADD( top, argl );
topextra = arg - top - atanPtr[(int32_t)index].One;
top3 = top - top2 + argl + topextra; bot =__FMADD( arg, atanPtr[(int32_t)index].One, 1.0 ); bot2 = __FMADD( argl, atanPtr[(int32_t)index].One, __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 - bot ) );
r = 1.0/bot;
quot = r*top2;
quot = __FMADD( __FNMSUB( bot, quot, top2), r, quot );
p29 = __FNMSUB( bot, quot, top2); p30 = quot*bot2; p33 = p29 - p30;
p35 = p33 + top3;
quot2 = p35*r;
quot2 = __FMADD( __FNMSUB( quot2, bot, p35 ), r, quot2 ); temp = 2.0*quot*quot2; argsq = __FMADD( quot, quot, temp ); argsq2 = __FMADD( quot2, quot2, __FMSUB( quot, quot, argsq) + temp ); sumt = __FMADD( argsq, atancoeff[7], atancoeff[6] ); sum = __FMADD( argsq, sumt, atancoeff[5] ); suml = __FMADD( argsq, sumt, atancoeff[5] - sum );
for (i = 4; i >= 1; i--) {
temps = __FMADD( sum, argsq, atancoeff[i] );
suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
sum = temps;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, quot, prod*quot2 ); temp = __FMADD( prod, quot, temp2 ); temp2 = __FMSUB( prod, quot, temp ) + temp2; sum = temp*kRecip; remainder = __FNMSUB( sum, atancoeff[0], temp ); partial = __FADD( remainder, temp2 );
residual = remainder - partial + temp2;
suml = __FMADD( partial, kRecip, (residual*kRecip) ); res = __FADD( quot, sum ); reslow = (quot - res + sum); resmid = __FADD( quot2, suml );
restiny = quot2 - resmid + suml;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
result = __FADD( atanPtr[(int32_t)index].Two, reshi );
reslow = __FADD( atanPtr[(int32_t)index].Three, resbot );
resmid = atanPtr[(int32_t)index].Two - result + reshi;
restiny = resbot - reslow + atanPtr[(int32_t)index].Three;
restemp = __FADD( resmid, reslow );
reshi = __FADD( result, restemp );
if (fabs(resmid) > fabs(reslow))
restiny = resmid - restemp + reslow + restiny;
else
restiny = reslow - restemp + resmid + restiny;
resbot = result - reshi + restemp + restiny;
if (xHead < 0.0) {
reshi = -reshi;
resbot = -resbot;
result = -result;
restemp = -restemp;
restiny = -restiny;
}
ddx.d[0] = __FADD( reshi, resbot );
ddx.d[1] = (reshi - ddx.d[0]) + resbot;
FESETENVD(fenv);
return (ddx.f); }
}
static const double kTiny = 1.805194375864829576e-276;
static const double kRescale = 2.923003274661805836e+48;
static const double kPi = 3.141592653589793116; static const double kPiTail1 = 1.224646799147353207e-16; static const double kPiTail1m = 1.224646799147352961e-16; static const double kPiTail2 = -2.994769809718339666e-33;
static const double k3PiBy4 = 2.356194490192344837; static const double k3PiBy4Tail1 = 9.184850993605148438e-17;
long double atan2l( long double y, long double x )
{
double p1, p2, p3, p2a, p4, p5, p6, p7, p8, p29, p30, p33, p35;
double p31, p32, p34, p36, p37, p38;
double quot, quot2;
double rarg, result, restemp;
double top, top2, topextra, top3;
double extra;
double bot, bot2;
double r, rarg2, rarg2a;
double temps;
double remainder, residual;
double partial;
double arg, argl, argsq, argsq2;
double temp,temp2;
double sum, sumt, suml;
double absx;
double xHead, xTail;
double res, reslow, resmid, restiny, reshi, resbot;
double p, q;
double prod, prodlow;
double xDen, xDenTail, yNum, yNumTail;
double frac3, frac2, frac, fract;
double pin, pin2, pin3;
double ressmall, restop, resint;
Double DeN ;
double fpenv;
DblDbl ddx, ddy;
int i, rflag;
uint32_t index;
uint32_t signx, signy, expx, expy, sigx, sigy;
struct atantblEntry *atanPtr = (struct atantblEntry *)atanTable - 6;
FEGETENVD(fpenv);
FESETENVD(0.0);
ddx.f = x; ddy.f = y;
xDen = ddx.d[0]; xDenTail = ddx.d[1]; yNum = ddy.d[0]; yNumTail = ddy.d[1];
expy = (ddy.i[0] >> 20) & 0x7ffu;
expx = (ddx.i[0] >> 20) & 0x7ffu; signy = (ddy.i[0] >> 31) & 1u;
signx = (ddx.i[0] >> 31) & 1u; sigy = (ddy.i[0] & 0xfffffu) | ddy.i[1];
sigx = (ddx.i[0] & 0xfffffu) | ddx.i[1];
if (expy == 0x7ffu) { if (sigy != 0) {
FESETENVD(fpenv);
return (y); }
else if (expx == 0x7ffu) { if (sigx != 0) {
FESETENVD(fpenv);
return (x); }
else if (signx == 0) { if (signy) { ddx.d[0] = -0.5*kPiBy2 ;
ddx.d[1] = -0.5*kPiBy2Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
else { ddx.d[0] = 0.5*kPiBy2 ;
ddx.d[1] = 0.5*kPiBy2Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
}
else { if (signy) { ddx.d[0] = -k3PiBy4 ;
ddx.d[1] = -k3PiBy4Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
else { ddx.d[0] = k3PiBy4 ;
ddx.d[1] = k3PiBy4Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
}
}
else { if (signy == 0.0) {
if (signx == 0.0) { ddx.d[0] = kPiBy2 ;
ddx.d[1] = kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
else { ddx.d[0] = kPiBy2 ;
ddx.d[1] = kPiBy2Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
}
else { if (signx == 0.0) { ddx.d[0] = -kPiBy2 ;
ddx.d[1] = -kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
else { ddx.d[0] = -kPiBy2 ;
ddx.d[1] = -kPiBy2Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
}
}
} if (expx == 0x7ffu) { if (sigx != 0) { FESETENVD(fpenv);
return (x);
}
else if (signx) { if (signy == 0) {
ddx.d[0] = 2.0*kPiBy2 ;
ddx.d[1] = 2.0*kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
else {
ddx.d[0] = -2.0*kPiBy2 ;
ddx.d[1] = -2.0*kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
}
else if (signy == 0) {
ddx.d[0] = 0.0 ;
ddx.d[1] = 0.0;
FESETENVD(fpenv);
return (ddx.f);
}
else {
ddx.d[0] = -0.0 ;
ddx.d[1] = -0.0;
FESETENVD(fpenv);
return (ddx.f);
}
} if ((sigx | expx) == 0u) { if ((sigy | expy) == 0u) { if (signx) {
if (signy) {
ddx.d[0] = -2.0*kPiBy2 ;
ddx.d[1] = -2.0*kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
else {
ddx.d[0] = 2.0*kPiBy2 ;
ddx.d[1] = 2.0*kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
}
else {
FESETENVD(fpenv);
return (y);
}
}
else if (signy) { if (signx) {
ddx.d[0] = -kPiBy2 ;
ddx.d[1] = -kPiBy2Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
else {
ddx.d[0] = -kPiBy2 ;
ddx.d[1] = -kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
}
else { if (signx) {
ddx.d[0] = kPiBy2 ;
ddx.d[1] = kPiBy2Tail1;
FESETENVD(fpenv);
return (ddx.f);
}
else {
ddx.d[0] = kPiBy2 ;
ddx.d[1] = kPiBy2Tail1m;
FESETENVD(fpenv);
return (ddx.f);
}
}
} if ((sigy | expy) == 0u) { if (signx == 0.0) { ddx.d[0] = yNum; ddx.d[1] = 0.0;
FESETENVD(fpenv);
return (ddx.f);
}
else { if (signy == 0) { ddx.d[0] = kPi ;
ddx.d[1] = kPiTail1m;
FESETENVD(fpenv);
return (ddx.f);
}
else { ddx.d[0] = -kPi ;
ddx.d[1] = -kPiTail1m;
FESETENVD(fpenv);
return (ddx.f);
}
}
}
if((((ddy.i[0]) & 0x7fffffff) < 0x06B00000)|| (((ddx.i[0]) & 0x7fffffff) < 0x06B00000)) {
xDen = xDen*kRescale; xDenTail = xDenTail*kRescale; yNum = yNum*kRescale; yNumTail = yNumTail*kRescale; r = 1.0/xDen;
if (fabs(r) < kTiny) { if (xDen > 0) { xHead = 0.0; xTail = 0.0;
frac3 = 0.0;
}
else { xHead = yNum; xTail = 0.0; frac3 = 0.0; }
}
}
r = 1.0/xDen; fract = __FMUL( yNum, r );
frac = __FMADD( r, __FNMSUB( xDen, fract, yNum ), fract ); p29 = __FNMSUB( xDen, frac, yNum ); p30 = __FMUL( frac, xDenTail ); p31 = __FNMSUB( frac, xDenTail, p30 ); p32 = p31;
p33 = __FSUB( p29, p30 );
p35 = __FADD( p33, yNumTail );
if (fabs(fract) == kInfinityu.f) {
xHead = fract;
xTail = 0.0;
frac3 = 0.0;
}
else if (fabs(fract)< kInfinityu.f) {
frac2 = __FMUL( p35, r );
frac2 = __FMADD( r, __FNMSUB( frac2, xDen, p35 ), frac2 ); if (fabs(p29) > fabs(p30))
p34 = p29 - p33 - p30 + p32;
else
p34 = p29 - (p33+p30)+p32;
p36 = __FNMSUB( frac2, xDen, p35 ) + __FNMSUB( frac2, xDenTail, p34 ); if (fabs(p33) > fabs(yNumTail))
p37 = p33 - p35 + yNumTail;
else
p37 = yNumTail - p35 + p33;
p38 = p37 + p36;
frac3 = p38*r/__FMADD( frac, frac, 1.0 ); xHead = frac; xTail = frac2; }
absx = fabs(xHead);
if (absx >= kLimit4) { reshi = kPiBy2; if (signx == 0.0) { resbot = kPiBy2Tail1m; ressmall = 0; if (xHead < 0.0) { reshi = -reshi;
resbot = -resbot;
}
}
else { reshi = kPiBy2; resbot = kPiBy2Tail1; ressmall = kPiBy2Tail2;
if (xHead < 0.0) {
reshi = -reshi;
resbot = -resbot;
ressmall = -ressmall;
}
}
}
rflag = 0;
if ((absx >= kLimit3) && (absx < kLimit4) ) { rarg = 1.0/xHead;
rflag = 1; p1 = __FNMSUB( xHead, rarg, 1.0 ); p2 = xTail*rarg; p3 = p1 - p2; rarg2a = p3*rarg;
rarg2 = rarg2a + __FNMSUB( xHead, rarg2a, p3 )*rarg; xHead = rarg;
xTail = rarg2;
}
if( (absx < kLimit1) || ((absx >= kLimit3) && (absx < kLimit4) )) { temp = 2.0*xHead*xTail; argsq = __FMADD( xHead, xHead, temp ); argsq2 = __FMADD( xTail, xTail, __FMSUB( xHead, xHead, argsq ) + temp ); sum = atancoeff[14];
sum = __FMADD( argsq, sum, atancoeff[13] ); sum = __FMADD( argsq, sum, atancoeff[12] ); sum = __FMADD( argsq, sum, atancoeff[11] ); sumt = __FMADD( argsq, sum, atancoeff[10] ); sum = __FMADD( argsq, sum, atancoeff[9] ); suml = __FMADD( argsq, sumt, atancoeff[9] - sum );
for (i = 8; i >= 1; i--) {
temps = __FMADD( argsq, sum, atancoeff[i] );
suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
sum=temps;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, xHead, prod*xTail ); temp = __FMADD( prod, xHead, temp2 ); temp2 = __FMSUB( prod, xHead, temp ) + temp2; sum = temp*kRecip; remainder = __FNMSUB( sum, atancoeff[0], temp ); partial = __FADD( remainder, temp2 );
residual = remainder - partial + temp2;
suml = __FMADD( partial, kRecip, (residual*kRecip) ); res = __FADD( xHead, sum ); reslow = (xHead - res + sum); resmid = __FADD( xTail, suml );
restiny = xTail - resmid + suml;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = res + p;
resbot = (res - reshi + p) + (q + restiny);
if (rflag == 1) { if (xHead > 0) { p1 = kPiBy2;
p2 = kPiBy2Tail1;
p2a = kPiBy2Tail2;
}
else {
p1 = -kPiBy2;
p2 = -kPiBy2Tail1;
p2a = -kPiBy2Tail2;
}
p3 = __FSUB( p1, reshi ); p4 = p1 - p3 - reshi;
p5 = __FSUB( p2, resbot );
p6 = p2 - p5 - resbot + p2a;
p7 = __FADD( p4, p5 );
reshi = p3 + p7;
if (fabs(p4) > fabs(p5))
p8 = p4 - p7 + p5;
else
p8 = p5 - p7 + p4;
resbot = __FADD( p3 - reshi + p7, (p6+p8-restiny) );
ressmall = (p3 - reshi + p7) - resbot + (p6 + p8 - restiny);
}
else ressmall = (res - reshi + p) - resbot + (q + restiny);
}
if( (absx >= kLimit1) && (absx < kLimit3) ) { DeN.f = __FMADD( 64.0, absx, kBig ); if (absx > kLimit2) { FESETENVD(kBig+1.0); DeN.f = absx + kBig; FESETENVD(kBig+1.0); DeN.f = absx + kBig; FESETENVD(kBig); index = DeN.i[1];
i = index - 1;
DeN.f = __FMADD( factor[i], absx, kBig+adjust[i] ); }
index = DeN.i[1];
if (xHead < 0) { arg = -xHead;
argl = -xTail;
}
else {
arg = xHead;
argl = xTail;
}
top = __FSUB( arg, atanPtr[(int32_t)index].One );
top2 = __FADD( top, argl );
topextra = arg - top - atanPtr[(int32_t)index].One;
top3 = top - top2 + argl + topextra; bot = __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 ); bot2 = __FMADD( argl, atanPtr[(int32_t)index].One, __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 - bot ) );
r = 1.0/bot;
quot = r*top2;
quot = __FMADD( __FNMSUB( bot, quot, top2 ) , r, quot );
p29 = __FNMSUB( bot, quot, top2 ); p30 = quot*bot2; p33 = p29-p30;
p35 = p33 + top3;
quot2 = p35*r;
quot2 = __FMADD( __FNMSUB( quot2, bot, p35 ), r, quot2 ); temp = 2.0*quot*quot2; argsq = __FMADD( quot, quot, temp ); argsq2 = __FMADD( quot2, quot2, __FMSUB( quot, quot, argsq) + temp ); sumt = __FMADD( argsq, atancoeff[7], atancoeff[6] ); sum = __FMADD( argsq, sumt, atancoeff[5] ); suml = __FMADD( argsq, sumt, atancoeff[5] - sum );
for (i = 4; i >= 1; i--) {
temps = __FMADD( sum, argsq, atancoeff[i] );
suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
sum = temps;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, quot, prod*quot2 ); temp = __FMADD( prod, quot, temp2 ); temp2 = __FMSUB( prod, quot, temp ) + temp2; sum = temp*kRecip; remainder = __FNMSUB( sum, atancoeff[0], temp ); partial = __FADD( remainder, temp2 );
residual = remainder - partial + temp2;
suml = __FMADD( partial, kRecip, (residual*kRecip) ); res = __FADD( quot, sum ); reslow = (quot - res + sum); resmid = __FADD( quot2, suml );
restiny = quot2 - resmid + suml;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
result = __FADD( atanPtr[(int32_t)index].Two, reshi );
reslow = __FADD( atanPtr[(int32_t)index].Three, resbot );
resmid = atanPtr[(int32_t)index].Two - result + reshi;
restiny = resbot - reslow + atanPtr[(int32_t)index].Three;
restemp = __FADD( resmid, reslow );
reshi = __FADD( result, restemp );
if (fabs(resmid) > fabs(reslow))
restiny = resmid - restemp + reslow + restiny;
else
restiny = reslow - restemp + resmid + restiny;
resbot = result - reshi + restemp + restiny;
if (xHead < 0.0) {
reshi = -reshi;
resbot = -resbot;
result = -result;
restemp = -restemp;
restiny = -restiny;
}
ressmall = (result - reshi + restemp) - resbot + restiny;
}
extra = ressmall + frac3;
if (signx == 0) { ddx.d[0] = __FADD( reshi, (resbot + extra) );
ddx.d[1] = (reshi - ddx.d[0]) + (resbot + extra);
}
else {
if (signy == 0) { pin = kPi; pin2 = kPiTail1; pin3 = kPiTail2;
}
else {
pin = -kPi; pin2 = -kPiTail1; pin3 = -kPiTail2;
}
restop = __FADD( pin, reshi );
resmid = pin - restop + reshi;
reslow = __FADD( pin2, resbot );
restiny = pin2 - reslow + resbot + extra;
resint = __FADD( resmid, reslow );
reshi = __FADD( restop, resint );
resbot = restop - reshi+resint;
if (fabs(resmid) > fabs(reslow))
resbot = resbot + (resmid-resint+reslow+restiny);
else
resbot = resbot + (reslow-resint+resmid+restiny);
ddx.d[0] = __FADD( reshi, resbot );
ddx.d[1] = (reshi - ddx.d[0]) + resbot;
}
FESETENVD(fpenv);
return (ddx.f);
}
#endif