#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
void _d3div(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low);
void _d3mul(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low);
void _d3add(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low);
long double _ExpInnerLD(double alpha, double beta, double gamma,
double *extra, int exptype);
static const uint32_t erfxtrData[] __attribute__ ((aligned(8))) = {
0x00000000, 0x00000000,
0xB9155555, 0x55555555,
0x38F99999, 0x9999999A
};
static const uint32_t erfpowData[] __attribute__ ((aligned(8))) = {
0x3FF00000, 0x00000000, 0x00000000, 0x00000000,
0xBFD55555, 0x55555555, 0xBC755555, 0x55555555,
0x3FB99999, 0x9999999A, 0xBC599999, 0x9999999A,
0xBF986186, 0x18618618, 0xBC386186, 0x18618613,
0x3F72F684, 0xBDA12F68, 0x3C12F684, 0xBDA15772,
0xBF48D301, 0x8D3018D3, 0xBB88D301, 0xBA2FF990,
0x3F1C01C0, 0x1C01C01C, 0x3B5C01FF, 0xEF4E45A0,
0xBEEBBD77, 0x9334EF0B, 0x3B84E5AE, 0x2FF08B10,
0x3EB87A00, 0x187A0018, 0x3B5EC7BC, 0xA752F683,
0xBE83777C, 0x55568CCD, 0xBB29453D, 0xDF9ADD3E,
0x3E4C2E30, 0x54870B52, 0xBADB98FC, 0x71C71C77,
0xBE12B673, 0x10AAA27C, 0x3AACE36B, 0x74F03296,
0x3DD6F448, 0xE13F19E7, 0xBA7D15ED, 0x097B425F,
0xBD9A289E, 0xE7F6D02A, 0xBA21F9AD, 0xD3C0CA47,
0x3D5BD577, 0xE7FBCC09, 0xB9E48B0F, 0xCD6E9E08,
0xBD1BC625, 0x25C514CA, 0xB9B161F9, 0xADD3C0CA,
0x3CDA173A, 0x4FC330FD, 0xB9648B0F, 0xCD6E9E03,
0xBC97270A, 0xC344BDA1, 0xB927B425, 0xED097B45,
0x3C537610, 0xC735BA78, 0x38D948B0, 0xFCD6E9DE,
0xBC0EF15A, 0xDD3C0CA4, 0xB8A61F9A, 0xDD3C0CA5,
0x3BC67194, 0x8B0FCD6F, 0xB8687E6B, 0x74F03292,
0xBB778E38, 0xE38E38E4, 0x381C71C7, 0x1C71C71C
};
static const uint32_t coeffData[] __attribute__ ((aligned(8))) = {
0x3FEFFFFF, 0xF87B4F81, 0xBC74B6DC, 0x82E2B33B,
0xBFDFFFF9, 0x651B9CCD, 0xBC77993A, 0x7E74C502,
0x3FE7FF4C, 0x63E2B209, 0x3C607D84, 0x888EE934,
0xBFFDF3AD, 0x307BCFA0, 0x3C80189F, 0x21239DB8,
0x4019F11D, 0x8747B2EA, 0x3CB282E6, 0x775914EA,
0xC03BFD72, 0x2186F691, 0x3CD1B8A8, 0xA64D28F3,
0x406123E2, 0xF4070A2E, 0xBD0AA69F, 0x832442C8,
0xC085E672, 0xCAE6452C, 0xBD264E3F, 0x5B4D5E70,
0x40AB48F2, 0xAF9DF505, 0x3D3CFDE8, 0x35619194,
0xC0CF99E0, 0xAC6AFA01, 0xBD525B2E, 0x7AEBA4B8,
0x40F07ABD, 0x450474EE, 0x3D9F4814, 0x21020E9F,
0xC10E4BB4, 0xEAA37C6D, 0x3DADDD80, 0x6CA8C8F2,
0x41282796, 0x01A20544, 0x3DC2C51E, 0x7D2AD0DC,
0xC1407A5E, 0xA1893145, 0xBDEB8D59, 0x4611C6D5,
0x4152F9A4, 0x9A0ED425, 0xBDF53C28, 0x002DF666,
0xC16227D6, 0x896B5A9A, 0x3DD50BAE, 0xCA782EDC,
0x416C4D5B, 0x1A0268E7, 0xBDE19C5A, 0x4D09E56C,
0xC1717E27, 0xC2AB2A36, 0x3E122CE5, 0xA75631BC,
0x41707C32, 0x89EC6733, 0xBE1F4B89, 0xEF2106D9,
0xC1663DCE, 0x61AF76F4, 0x3E028503, 0xD6CDBB9E,
0x41531F97, 0x0C38EEF3, 0xBDE20467, 0x5FB2BCF0,
0xC12F7BE6, 0x257D1D65, 0x3D8DFD02, 0xE8C06000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x3FEFFFFF, 0xFFEBC039, 0x3C529C69, 0x6CF540A8,
0xBFDFFFFF, 0xE38F77CF, 0x3BFCDCBB, 0x0CD77D00,
0x3FE7FFFB, 0x2DCCB103, 0x3C8175D4, 0xB96461CA,
0xBFFDFF79, 0x49317263, 0xBC906876, 0xAC321F5A,
0x401A3AAB, 0x9BF9D0BD, 0x3CABAD03, 0xD57E6058,
0xC03D5E32, 0x130E26B3, 0xBC9E17BB, 0x1F622DC0,
0x4063C8DE, 0x39F025C5, 0x3CEE695F, 0x35A91AF0,
0xC08E3E6E, 0x953DCEB7, 0x3D24E053, 0xE29EAEBE,
0x40B8AA52, 0x3AD56173, 0x3D50EDEF, 0x7FD5BCF1,
0xC0E440E6, 0xCEADBFBB, 0x3D71A20A, 0x51AF181A,
0x410FEE4B, 0xDCA4E1E0, 0x3D8CCED6, 0xFF021BF0,
0xC13752BC, 0xAAA09717, 0xBDDB430A, 0x583D9DFE,
0x415EC821, 0x1EEE2921, 0x3DDA8E50, 0xAC97C8D8,
0xC182031E, 0xEDA7A8C5, 0x3E23FDA5, 0x816A17A9,
0x41A26A35, 0x43351D14, 0x3E4A1C99, 0x4FE13BE4,
0xC1C0397D, 0x3FF8253E, 0xBE55B81D, 0xBAD43C68,
0x41D84B8E, 0x7B8414AB, 0x3E797BE5, 0x69A73302,
0xC1EE6B6B, 0x349E71AE, 0xBE7EB248, 0x45D29614,
0x41FF34DA, 0xF64919BF, 0xBE87460E, 0x4926C028,
0xC20984E7, 0x864FAEBA, 0x3EA1DD11, 0x07234030,
0x420FF82D, 0xCF64D788, 0x3E9B32FC, 0x90159470,
0xC20CCAA6, 0x72BB40E4, 0x3E619A97, 0x18E25640,
0x420096A8, 0x68135449, 0xBEAFFD28, 0xED6BBBFA,
0xC1E25D70, 0xCAB7C34C, 0xBE5B42DB, 0xEC71D440,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x3FEFFFFF, 0xFFFFE639, 0xBC76B4F8, 0xF22A1D6A,
0xBFDFFFFF, 0xFFC4926E, 0x3C74AA2B, 0x05602AA4,
0x3FE7FFFF, 0xEF8F13A8, 0xBC802572, 0x96640478,
0xBFFDFFFD, 0x175A90A3, 0x3C8448AC, 0x5447AC1C,
0x401A3FD0, 0xA3F30955, 0x3CBDC702, 0x0AF0E74C,
0xC03D85B0, 0xDE15683B, 0x3CD775BC, 0x8D3FFDA1,
0x406441FD, 0xA9731CF2, 0xBCEC6B02, 0xE82A19B0,
0xC0904FCE, 0xB4270577, 0x3CF73CA2, 0xBC7037A0,
0x40BDA90A, 0xA9B95F59, 0xBD5B5C9C, 0x5ABA6422,
0xC0ED1C94, 0x9666E63C, 0xBD5C789C, 0x5A7928DC,
0x411D56F4, 0x20DB5218, 0x3DB320CC, 0xDE074D1E,
0xC14CEC6A, 0x9FD76B11, 0x3DE9DE29, 0x03D42F2C,
0x417AC1F5, 0xF5446DF4, 0xBE031819, 0x3468BE58,
0xC1A67958, 0xAC7171B2, 0xBE3581A4, 0xB51D7A74,
0x41D0B077, 0x6E96B0B2, 0x3E574081, 0x998F5AA4,
0xC1F56C00, 0xFD9F8D53, 0xBE94D297, 0xFD7D3E08,
0x42173F0C, 0x0EB8D9D7, 0x3EBE2714, 0x71045A4A,
0xC234D11D, 0x4DA55FFA, 0x3EC4EA18, 0x03DADE62,
0x424DDAFE, 0x5B2DFE46, 0x3EE6D6D1, 0x58677012,
0xC2607209, 0x5B35C7C9, 0x3F04E6A4, 0x0B2429BA,
0x426A1766, 0x45D71E26, 0xBF0B66BC, 0xFE7F12C0,
0xC26A80F6, 0x5FAAB021, 0x3F042804, 0xC0D4DFB5,
0x4259DF2B, 0x9341D4CC, 0xBEE2816A, 0x58CBD636,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x3FF00000, 0x00000000, 0xBB75CA29, 0xC48EA000,
0xBFE00000, 0x00000000, 0x3C1A999F, 0xD44FCE00,
0x3FE7FFFF, 0xFFFFFFFE, 0x3C642A4A, 0xB64D8068,
0xBFFDFFFF, 0xFFFFFE99, 0xBC9ABDD0, 0x98308E2C,
0x401A3FFF, 0xFFFFA343, 0x3CB85B18, 0x87B93172,
0xC03D87FF, 0xFFEDDBFB, 0x3CD309B5, 0xEEE35E41,
0x40644D7F, 0xFE9A5199, 0x3CFC5B15, 0x3519E834,
0xC0907EF7, 0xE9B322F4, 0xBD35FCB1, 0x58754A9E,
0x40BEEE0E, 0xB41A926D, 0xBD57FB38, 0x347B51CF,
0xC0F06E6C, 0x572C5E36, 0x3D720A1F, 0x8B0D1254,
0x412382B8, 0xB55E07E1, 0xBDCABBAE, 0x93440AE8,
0xC1599872, 0x5E40B785, 0x3DF51615, 0xACBF5F19,
0x41925B4B, 0x43808517, 0xBE3C17EF, 0x5CEF40F8,
0xC1CC74A4, 0x3D911FB7, 0x3E172705, 0xD877BCC0,
0x42077587, 0x54FEECDC, 0xBEA98313, 0x07804B47,
0xC2441B83, 0xD7BACF31, 0xBEEA3073, 0x0D0DBA05,
0x428165FB, 0x78789F6C, 0xBED45B8E, 0xD3B8907C,
0xC2BD68C3, 0xB0580084, 0x3F1BFE68, 0xE7B64DE0,
0x42F78115, 0xFD9904B1, 0xBF8FED1E, 0x29D74809,
0xC3313FA3, 0xDBF81756, 0x3FDC570B, 0xC88AFF19,
0x4366A5A6, 0xB262F5E4, 0x3FFA9257, 0x0A3AD3F8,
0xC399F8EF, 0xB5A4E015, 0x401B7ECE, 0x1B301370,
0x43C9687F, 0xB88D7E83, 0xC06CE8BA, 0x85792A6C,
0xC3F4A8E4, 0x451D5574, 0x409B5C8E, 0x70046666,
0x441B0DCA, 0xB2FAB836, 0x40996671, 0xF8BB2152,
0xC43B52DB, 0xC1E7B2EE, 0x40DB5AF6, 0x615EB7EE,
0x4453ECD3, 0xF252D951, 0x40FC199B, 0x4FAD14EA,
0xC462A264, 0xDCCB32FF, 0xC0E17D1D, 0x52615E34,
0x4460C115, 0xAC203924, 0x40E780D0, 0xA1F70E00,
0x3FF00000, 0x00000000, 0x3A9A41C2, 0x18000000,
0xBFE00000, 0x00000000, 0xBB5DCD07, 0x8C200000,
0x3FE80000, 0x00000000, 0x3C0EBA46, 0xBCD71000,
0xBFFE0000, 0x00000001, 0xBC88EE2A, 0x6486B1E4,
0x401A4000, 0x00000101, 0x3CB25362, 0xC332A010,
0xC03D8800, 0x00009B65, 0x3CD119CB, 0x15915886,
0x40644D80, 0x0022EC9B, 0xBCE9E9C4, 0x0897A1AC,
0xC0907EF8, 0x05F94723, 0x3D2AB506, 0xB989F798,
0x40BEEE12, 0x939C8382, 0x3D5AD230, 0x743AEFD8,
0xC0F06E8D, 0xB86AB92B, 0x3D98E5AE, 0xBE0C9174,
0x412384D6, 0x161AD32F, 0x3DBB5BF4, 0xEF6A7DBF,
0xC159B642, 0x7281154F, 0x3DFECD4C, 0xAA7F47EC,
0x419305F9, 0xC7BF95B4, 0xBE349B0E, 0x9D9581DE,
0xC1D12C24, 0x5D3084BD, 0xBE5E4995, 0xB96B303C,
0x4215451C, 0xDF186F15, 0x3EA8548E, 0x669919A8,
0xC25F3153, 0x01C60944, 0x3ED326BD, 0x842F8568,
0x42A25C17, 0xD2287EF6, 0x3F493E23, 0x133CA4B2,
0xC2D613F5, 0x225EEB4A, 0x3F4A86AC, 0x4B33D440,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0xBFA4AFAB, 0x9B6E1105, 0x3C463A99, 0xD0AF98EE,
0xBFD4F526, 0x85059209, 0xBC73ECA6, 0xE2948B70,
0xBFEA28C7, 0xA68495F3, 0xBC81EAF6, 0x808C4753,
0xBFA980E7, 0x79855BE9, 0x3C30F43B, 0xBDDFAE98,
0x3F8895D7, 0xF2B826D4, 0x3C250FDF, 0xAF15221E,
0xBF6222D1, 0xDE8EFFF2, 0x3BFC3D53, 0x8B04AE74,
0x3F26522D, 0x3A826B8E, 0xBBBFB440, 0x1723702C,
0x3F15289F, 0x67803BEA, 0x3BBA9594, 0x12B030A4,
0xBF08627F, 0xFCF37770, 0x3BAD7A54, 0xCCB615F2,
0x3EEBEB7F, 0xBC45B14A, 0xBB715091, 0x8DB6D7B4,
0xBEBE0FDE, 0x30CC0705, 0xBB3ED352, 0x5BE0D640,
0xBE9B7052, 0x9A2F9B59, 0xBB29B1CC, 0xF66FE0B8,
0x3E97B4BD, 0x41755DF3, 0xBB10EE4A, 0xCBA6B5AC,
0xBE807536, 0x0AA20FDE, 0xBAE9E0B1, 0xF22797C0,
0x3E57D20D, 0xE4C07B40, 0x3AFADF83, 0x4F8A74E3,
0x3E1E7ADC, 0x1A5BC732, 0xBAACBD1B, 0xDB944806,
0xBE29CB52, 0xFF49FDEF, 0xBACE7432, 0x1EF280D9,
0x3E14D367, 0x091F84DD, 0x3ABD9BF9, 0x8C93A4B0,
0xBDF32CB0, 0xA2D996E6, 0x3A947D3D, 0xF3852A64,
0x3DC1332B, 0x63A393F5, 0xBA666775, 0x70308F4A,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0xBFD09B84, 0xB58C74B1, 0x3C649EAC, 0xBE7125A0,
0xBFE04818, 0x9E7A4AF5, 0xBC7247BB, 0x9C14B6B2,
0xBFEC9082, 0x7462A177, 0xBC65A4D5, 0x2B14D05A,
0xBF9A71A4, 0xBDE00D7D, 0xBC3D297D, 0x6BADB471,
0x3F789A87, 0xB48FE138, 0x3BF9AE10, 0x1027D54D,
0xBF549926, 0x2376DBFA, 0x3BDC0F4D, 0x2A750CAB,
0x3F2C4E86, 0xAA32FCDB, 0xBBCBD7A9, 0xE43684BF,
0xBEF790BC, 0x075477C9, 0x3B89898A, 0xADCD6E9C,
0xBEC63409, 0x8491E3A5, 0xBB572A20, 0x12F684BF,
0x3EC3E354, 0x0DCFCA75, 0x3B485BE3, 0x87E6B74C,
0xBEAB790E, 0xE236BFDB, 0x3B40DEF2, 0x9161F9AE,
0x3E89DA01, 0x98D7CE6C, 0xBB2327A2, 0xC3F35BA6,
0xBE5DEF8D, 0x200AB94C, 0x3AF95BA7, 0x81948B10,
0xBE129E35, 0x62080668, 0x3AB6F684, 0xBDA12F68,
0x3E24ACF9, 0x5B57BC54, 0x3AC55555, 0x55555556,
0xBE10CB5E, 0x9DE40846, 0x3AAC71C7, 0x1C71C71F,
0x3DF1C294, 0xD79A2DA8, 0xBA9F9ADD, 0x3C0CA45A,
0xBDC806EC, 0x95F17FF3, 0xBA66E9E0, 0x6522C3F3,
0x3D63D790, 0x6BE97B42, 0x3A07B425, 0xED097B40,
0x3D8A0F87, 0x44E09E06, 0x3A248B0F, 0xCD6E9E04,
0xBD780C79, 0x96291620, 0x3A1948B0, 0xFCD6E9E0,
0x3D5B9C07, 0xD9F03291, 0x39F87E6B, 0x74F03291,
0xBD34D565, 0xE06522C4, 0x39A948B0, 0xFCD6E9DF,
0x3CF1D800, 0xD6E9E065, 0x398161F9, 0xADD3C0C9,
0x3CF13B9A, 0xDD3C0CA4, 0x39961F9A, 0xDD3C0CA4,
0xBCE41FE1, 0xF9ADD3C1, 0x397ADD3C, 0x0CA45880,
0x3CC7D948, 0xB0FCD6EA, 0xB94F9ADD, 0x3C0CA458,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0x00000000, 0x00000000, 0x00000000, 0x00000000,
0xBF7E4EAA, 0x825F8588, 0x3C0E9883, 0x7D311FB0,
0xBF9B20DB, 0x38454E52, 0x3C282460, 0xA5A22510,
0xBFE92D49, 0x5DE062C7, 0x3C842242, 0x550B45F8,
0xBFAE856A, 0x106722F3, 0xBC3E5242, 0x6EC1567E,
0x3F8CF900, 0xD71B5137, 0xBC29B9D7, 0x4D921F9C,
0xBF630F14, 0x2DC98724, 0x3BEFD901, 0xC6097B43,
0x3EE6A533, 0x044E582C, 0x3B653307, 0x1C71C6A0,
0x3F261616, 0xBB9808F2, 0x3BC00838, 0xA4587E6A,
0xBF12A42E, 0x193904D9, 0x3B7FAA5E, 0xD097B432,
0x3EEFFC81, 0x1FC57602, 0xBB657CD6, 0xE9E0651E,
0x3E8108C1, 0x35437559, 0xBAFF9ADD, 0x3C0CA484,
0xBEBA9B0B, 0xE9B78F13, 0x3B555555, 0x55555554,
0x3EA7C3AD, 0xEB179410, 0xBB3948B0, 0xFCD6E9DF,
0xBE853AAC, 0x329694F0, 0xBB1948B0, 0xFCD6E9E2,
0xBE27591E, 0x10C71C72, 0x3ABC71C7, 0x1C71C700,
0x3E53D507, 0xBCC71C72, 0xBAEC71C7, 0x1C71C720,
0xBE4220D2, 0xDA12F685, 0x3AE097B4, 0x25ED097C,
0x3E2224EC, 0xBDA12F68, 0x3AC2F684, 0xBDA12F68,
0xBDF19BC0, 0xCA4587E7, 0x3A922C3F, 0x35BA781A
};
typedef double DD[2];
typedef DD DDarray[29];
static const Double ca = {{0x3FF20DD7, 0x50429B6D}};
static const Double cb = {{0x3C71AE3A, 0x914FED80}};
static const Double cc = {{0xB903CBBE, 0xBF65F145}};
static long double ErfCommon( double dhead, double dtail, int ierf )
{
double temp, temp2, temp3, argsq, argsq2, sum, suml, suminf, reslow, resmid, restiny;
double q, p, prodlow, addend, res, prod1,prod2,px, center, adjust, times;
double oarg, oarg2, a1, a2, a3, ab, ac, asq1,asq2,asq, arg2,argx;
double rx, arg, as, asx, resx, qx, exp;
DblDbl resdd;
int i, ir, is, ie;
double *erfxtr = (double *)erfxtrData;
DD *erfpow = (DD *)erfpowData;
DDarray *coeff = (DDarray *)(coeffData - 116);
if (fabs(dhead) <= 0.75) {
temp = (2.0*dhead)*dtail; argsq = __FMADD( dhead, dhead, temp ); argsq2 = __FMSUB( dhead, dhead, argsq + temp ) + __FMSUB( (2.0*dhead), dtail, temp ); sum = erfpow[21][0]; suml = erfpow[21][1];
for (i = 20; i >= 3; --i) { temp = __FMADD( sum, argsq, erfpow[i][0] );
suml = __FMADD( sum, argsq, erfpow[i][0] - temp ) +
__FMADD( suml, argsq, __FMADD( sum, argsq2, erfpow[i][1] ) );
sum = temp;
}
suminf = 0.0e0;
for (i = 2; i >= 0; --i) {
prodlow = __FMADD( suml, argsq, (sum*argsq2) );
addend = __FMSUB( suml, argsq, (suml*argsq) ) + __FMSUB( sum, argsq2, (sum*argsq2) ) +
__FMADD( suminf, argsq, __FMADD( sum, argsq2, __FMSUB( suml, argsq, prodlow ) ) );
temp = __FMADD( sum, argsq, prodlow ); temp2 = __FADD( __FMSUB( sum, argsq, temp), prodlow ); temp3 = __FMSUB( sum, argsq, temp) - temp2 + prodlow + addend; res = __FADD( erfpow[i][0], temp );
reslow = (erfpow[i][0] - res + temp);
resmid = __FADD( erfpow[i][1], temp2 );
restiny = erfpow[i][1] - resmid + temp2 + erfxtr[i] + temp3;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; sum = __FADD( res, p );
suml = __FADD( (res - sum + p), (q + restiny) );
suminf = (res - sum + p) - suml + (q + restiny);
}
_d3mul( sum, suml, suminf, dhead, dtail, 0.0, &prod1, &prod2, &px );
_d3mul( prod1, prod2, px, ca.f, cb.f, cc.f, &(resdd.d[0]), &(resdd.d[1]), &qx );
if (ierf == 0) {
_d3add( 1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -qx,
&(resdd.d[0]), &(resdd.d[1]), &qx );
}
return resdd.f;
}
else if (fabs(dhead) <= 2.0) {
if (fabs(dhead) < 1.00) { center = 0.625; adjust = 0.28125; times = 2.0; ir = 8;
is = 18;
ie = 15;
}
else if (fabs(dhead) < 1.25) {
center = 0.8125;
adjust = 0.28125;
times = 2.0;
ir = 6;
is = 19;
ie = 16;
}
else {
center = 1.5;
adjust = 1.375;
times = 3.0;
ir = 7;
is = 26;
ie = 14;
}
if (dhead < 0.0) {
oarg = dhead;
oarg2 = dtail;
temp = -center - dhead; arg2 = -dtail;
}
else {
temp = dhead - center;
arg2 = dtail;
oarg = -dhead;
oarg2 = -dtail;
}
arg = __FADD( temp, arg2 );
arg2 = temp - arg + arg2;
sum = coeff[ir][is][0];
suml = coeff[ir][is][1];
for (i = is-1; i >= ie; --i) { temp = __FMADD( sum, arg, coeff[ir][i][0] );
suml = __FMADD( sum, arg, coeff[ir][i][0] - temp ) +
__FMADD( suml, arg, __FMADD( sum, arg2, coeff[ir][i][1] ) );
sum = temp;
}
suminf = 0.0e0; for (i = ie-1; i >= 0; --i) { _d3mul( sum, suml, suminf, arg, arg2, 0.0, &prod1, &prod2, &px );
_d3add( prod1, prod2, px, coeff[ir][i][0], coeff[ir][i][1], 0.0,
&sum, &suml, &suminf );
}
a2 = times*oarg2;
a3 = __FMSUB( times, oarg2, a2 ); a1 = oarg*times; ab = __FADD( __FMSUB( times, oarg2, a1 ), a2 ); ac = __FMSUB( times, oarg2, a1 ) - ab + a2 + a3; _d3add( sum, suml, suminf, adjust, 0.0, 0.0, &sum, &suml, &as );
_d3add( sum, suml, as, a1, ab, ac, &sum, &suml, &asx );
resdd.f = _ExpInnerLD(sum, suml, asx, &resx, 4);
if (ierf == 0) {
if (dhead <= 0.0) _d3add( 2.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -resx,
&(resdd.d[0]), &(resdd.d[1]), &qx );
}
else { if (dhead > 0.0)
_d3add( 1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -resx,
&(resdd.d[0]), &(resdd.d[1]), &qx );
else
_d3add(-1.0, 0.0, 0.0, resdd.d[0], resdd.d[1], resx,
&(resdd.d[0]), &(resdd.d[1]), &qx );
}
return resdd.f;
}
else {
_d3mul( dhead, dtail, 0.0, dhead, dtail, 0.0,
&asq1, &asq2, &asq ); _d3div( 1.0, 0.0, 0.0, asq1, asq2, asq,
&arg, &arg2, &argx ); if (arg >= 0.18e0) { ir = 1; is = 21; ie = 10; }
else if (arg >= 0.11e0) {
ir = 2; is = 23;
ie = 12;
}
else if (arg >= 0.06e0) {
ir = 3; is = 22;
ie = 6;
if (ierf != 1) ie = 11;
}
else if (arg >= 0.01e0) {
ir = 4; is = 28;
ie = 10;
}
else if (arg >= 0.00138e0) {
ir = 5; is = 17;
ie = 4;
if (ierf != 1) ie = 8;
}
else {
resdd.d[0] = 0.0;
resdd.d[1] = 0.0;
rx = 0.0;
if (dhead == dhead) goto commontail;
resdd.d[0] = dhead;
return resdd.f;
}
if ((ierf == 1) && (arg < .01)) { resdd.d[0] = 0.0;
resdd.d[1] = 0.0; rx = 0.0; goto commontail;
}
sum = coeff[ir][is][0];
suml = coeff[ir][is][1]; for (i = is-1; i >= ie; --i) { temp = __FMADD( sum, arg, coeff[ir][i][0] );
suml = __FMADD( sum, arg, coeff[ir][i][0] - temp ) +
__FMADD( suml, arg, __FMADD( sum, arg2, coeff[ir][i][1] ) );
sum = temp;
}
suminf = 0.0e0; for (i = ie-1; i >= 0; --i) { _d3mul( sum, suml, suminf, arg, arg2, argx, &prod1, &prod2, &px);
_d3add( prod1, prod2, px, coeff[ir][i][0], coeff[ir][i][1], 0.0,
&sum, &suml, &suminf );
}
_d3mul( sum, suml, suminf, ca.f, cb.f, cc.f,
&prod1, &prod2, &px ); resdd.f = _ExpInnerLD(-asq1, -asq2, -asq, &exp, 4);
_d3mul( prod1, prod2, px, resdd.d[0], resdd.d[1], exp,
&prod1, &prod2, &px );
_d3div( prod1, prod2, px, dhead, dtail, 0.0,
&(resdd.d[0]), &(resdd.d[1]), &rx );
resdd.f = 0.5 * resdd.f;
rx = 0.5*rx;
commontail:
if (dhead > 0.0) { if (ierf == 1) _d3add( 1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -rx,
&(resdd.d[0]), &(resdd.d[1]), &qx );
}
else {
if (ierf == 1)
_d3add(-1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -rx,
&(resdd.d[0]), &(resdd.d[1]), &qx );
else
_d3add( 2.0, 0.0, 0.0, resdd.d[0], resdd.d[1], rx,
&(resdd.d[0]), &(resdd.d[1]), &qx );
}
return resdd.f;
}
}
long double erfl( long double x )
{
double head, fenv;
DblDbl t;
FEGETENVD(fenv);
FESETENVD(0.0);
t.f = x;
t.f = ErfCommon(t.d[0], t.d[1], 1);
head = __FADD( t.d[0], t.d[1] );
t.d[1] = (t.d[0] - head) + t.d[1]; t.d[0] = head;
FESETENVD(fenv);
return t.f;
}
long double erfcl( long double x )
{
double head, fenv;
DblDbl t;
FEGETENVD(fenv);
FESETENVD(0.0);
t.f = x;
t.f = ErfCommon(t.d[0], t.d[1], 0);
head = __FADD( t.d[0], t.d[1] );
t.d[1] = (t.d[0] - head) + t.d[1]; t.d[0] = head;
FESETENVD(fenv);
return t.f;
}
#endif