125c28e8Piotr Jasiukajtis/*
225c28e8Piotr Jasiukajtis * CDDL HEADER START
325c28e8Piotr Jasiukajtis *
425c28e8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e8Piotr Jasiukajtis *
825c28e8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e8Piotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e8Piotr Jasiukajtis * and limitations under the License.
1225c28e8Piotr Jasiukajtis *
1325c28e8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e8Piotr Jasiukajtis *
1925c28e8Piotr Jasiukajtis * CDDL HEADER END
2025c28e8Piotr Jasiukajtis */
2125c28e8Piotr Jasiukajtis
2225c28e8Piotr Jasiukajtis/*
2325c28e8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e8Piotr Jasiukajtis */
2525c28e8Piotr Jasiukajtis/*
2625c28e8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e8Piotr Jasiukajtis * Use is subject to license terms.
2825c28e8Piotr Jasiukajtis */
2925c28e8Piotr Jasiukajtis
30ddc0e0bRichard Lowe#pragma weak __tgammal = tgammal
3125c28e8Piotr Jasiukajtis
3225c28e8Piotr Jasiukajtis#include "libm.h"
3325c28e8Piotr Jasiukajtis#include <sys/isa_defs.h>
3425c28e8Piotr Jasiukajtis
3525c28e8Piotr Jasiukajtis#if defined(_BIG_ENDIAN)
3625c28e8Piotr Jasiukajtis#define	H0_WORD(x)	((unsigned *) &x)[0]
3725c28e8Piotr Jasiukajtis#define	H3_WORD(x)	((unsigned *) &x)[3]
3825c28e8Piotr Jasiukajtis#define	CHOPPED(x)	(long double) ((double) (x))
3925c28e8Piotr Jasiukajtis#else
4025c28e8Piotr Jasiukajtis#define	H0_WORD(x)	((((int *) &x)[2] << 16) | \
4125c28e8Piotr Jasiukajtis			(0x0000ffff & (((unsigned *) &x)[1] >> 15)))
4225c28e8Piotr Jasiukajtis#define	H3_WORD(x)	((unsigned *) &x)[0]
4325c28e8Piotr Jasiukajtis#define	CHOPPED(x)	(long double) ((float) (x))
4425c28e8Piotr Jasiukajtis#endif
4525c28e8Piotr Jasiukajtis
4625c28e8Piotr Jasiukajtisstruct LDouble {
4725c28e8Piotr Jasiukajtis	long double h, l;
4825c28e8Piotr Jasiukajtis};
4925c28e8Piotr Jasiukajtis
5025c28e8Piotr Jasiukajtis/* INDENT OFF */
5125c28e8Piotr Jasiukajtis/* Primary interval GTi() */
5225c28e8Piotr Jasiukajtisstatic const long double P1[] = {
5325c28e8Piotr Jasiukajtis	+0.709086836199777919037185741507610124611513720557L,
5425c28e8Piotr Jasiukajtis	+4.45754781206489035827915969367354835667391606951e-0001L,
5525c28e8Piotr Jasiukajtis	+3.21049298735832382311662273882632210062918153852e-0002L,
5625c28e8Piotr Jasiukajtis	-5.71296796342106617651765245858289197369688864350e-0003L,
5725c28e8Piotr Jasiukajtis	+6.04666892891998977081619174969855831606965352773e-0003L,
5825c28e8Piotr Jasiukajtis	+8.99106186996888711939627812174765258822658645168e-0004L,
5925c28e8Piotr Jasiukajtis	-6.96496846144407741431207008527018441810175568949e-0005L,
6025c28e8Piotr Jasiukajtis	+1.52597046118984020814225409300131445070213882429e-0005L,
6125c28e8Piotr Jasiukajtis	+5.68521076168495673844711465407432189190681541547e-0007L,
6225c28e8Piotr Jasiukajtis	+3.30749673519634895220582062520286565610418952979e-0008L,
6325c28e8Piotr Jasiukajtis};
6425c28e8Piotr Jasiukajtisstatic const long double Q1[] = {
6525c28e8Piotr Jasiukajtis	+1.0+0000L,
6625c28e8Piotr Jasiukajtis	+1.35806511721671070408570853537257079579490650668e+0000L,
6725c28e8Piotr Jasiukajtis	+2.97567810153429553405327140096063086994072952961e-0001L,
6825c28e8Piotr Jasiukajtis	-1.52956835982588571502954372821681851681118097870e-0001L,
6925c28e8Piotr Jasiukajtis	-2.88248519561420109768781615289082053597954521218e-0002L,
7025c28e8Piotr Jasiukajtis	+1.03475311719937405219789948456313936302378395955e-0002L,
7125c28e8Piotr Jasiukajtis	+4.12310203243891222368965360124391297374822742313e-0004L,
7225c28e8Piotr Jasiukajtis	-3.12653708152290867248931925120380729518332507388e-0004L,
7325c28e8Piotr Jasiukajtis	+2.36672170850409745237358105667757760527014332458e-0005L,
7425c28e8Piotr Jasiukajtis};
7525c28e8Piotr Jasiukajtisstatic const long double P2[] = {
7625c28e8Piotr Jasiukajtis	+0.428486815855585429730209907810650135255270600668084114L,
7725c28e8Piotr Jasiukajtis	+2.62768479103809762805691743305424077975230551176e-0001L,
7825c28e8Piotr Jasiukajtis	+3.81187532685392297608310837995193946591425896150e-0002L,
7925c28e8Piotr Jasiukajtis	+3.00063075891811043820666846129131255948527925381e-0003L,
8025c28e8Piotr Jasiukajtis	+2.47315407812279164228398470797498649142513408654e-0003L,
8125c28e8Piotr Jasiukajtis	+3.62838199917848372586173483147214880464782938664e-0004L,
8225c28e8Piotr Jasiukajtis	+3.43991105975492623982725644046473030098172692423e-0006L,
8325c28e8Piotr Jasiukajtis	+4.56902151569603272237014240794257659159045432895e-0006L,
8425c28e8Piotr Jasiukajtis	+2.13734755837595695602045100675540011352948958453e-0007L,
8525c28e8Piotr Jasiukajtis	+9.74123440547918230781670266967882492234877125358e-0009L,
8625c28e8Piotr Jasiukajtis};
8725c28e8Piotr Jasiukajtisstatic const long double Q2[] = {
8825c28e8Piotr Jasiukajtis	+1.0L,
8925c28e8Piotr Jasiukajtis	+9.18284118632506842664645516830761489700556179701e-0001L,
9025c28e8Piotr Jasiukajtis	-6.41430858837830766045202076965923776189154874947e-0003L,
9125c28e8Piotr Jasiukajtis	-1.24400885809771073213345747437964149775410921376e-0001L,
9225c28e8Piotr Jasiukajtis	+4.69803798146251757538856567522481979624746875964e-0003L,
9325c28e8Piotr Jasiukajtis	+7.18309447069495315914284705109868696262662082731e-0003L,
9425c28e8Piotr Jasiukajtis	-8.75812626987894695112722600697653425786166399105e-0004L,
9525c28e8Piotr Jasiukajtis	-1.23539972377769277995959339188431498626674835169e-0004L,
9625c28e8Piotr Jasiukajtis	+3.10019017590151598732360097849672925448587547746e-0005L,
9725c28e8Piotr Jasiukajtis	-1.77260223349332617658921874288026777465782364070e-0006L,
9825c28e8Piotr Jasiukajtis};
9925c28e8Piotr Jasiukajtisstatic const long double P3[] = {
10025c28e8Piotr Jasiukajtis	+0.3824094797345675048502747661075355640070439388902L,
10125c28e8Piotr Jasiukajtis	+3.42198093076618495415854906335908427159833377774e-0001L,
10225c28e8Piotr Jasiukajtis	+9.63828189500585568303961406863153237440702754858e-0002L,
10325c28e8Piotr Jasiukajtis	+8.76069421042696384852462044188520252156846768667e-0003L,
10425c28e8Piotr Jasiukajtis	+1.86477890389161491224872014149309015261897537488e-0003L,
10525c28e8Piotr Jasiukajtis	+8.16871354540309895879974742853701311541286944191e-0004L,
10625c28e8Piotr Jasiukajtis	+6.83783483674600322518695090864659381650125625216e-0005L,
10725c28e8Piotr Jasiukajtis	-1.10168269719261574708565935172719209272190828456e-0006L,
10825c28e8Piotr Jasiukajtis	+9.66243228508380420159234853278906717065629721016e-0007L,
10925c28e8Piotr Jasiukajtis	+2.31858885579177250541163820671121664974334728142e-0008L,
11025c28e8Piotr Jasiukajtis};
11125c28e8Piotr Jasiukajtisstatic const long double Q3[] = {
11225c28e8Piotr Jasiukajtis	+1.0L,
11325c28e8Piotr Jasiukajtis	+8.25479821168813634632437430090376252512793067339e-0001L,
11425c28e8Piotr Jasiukajtis	-1.62251363073937769739639623669295110346015576320e-0002L,
11525c28e8Piotr Jasiukajtis	-1.10621286905916732758745130629426559691187579852e-0001L,
11625c28e8Piotr Jasiukajtis	+3.48309693970985612644446415789230015515365291459e-0003L,
11725c28e8Piotr Jasiukajtis	+6.73553737487488333032431261131289672347043401328e-0003L,
11825c28e8Piotr Jasiukajtis	-7.63222008393372630162743587811004613050245128051e-0004L,
11925c28e8Piotr Jasiukajtis	-1.35792670669190631476784768961953711773073251336e-0004L,
12025c28e8Piotr Jasiukajtis	+3.19610150954223587006220730065608156460205690618e-0005L,
12125c28e8Piotr Jasiukajtis	-1.82096553862822346610109522015129585693354348322e-0006L,
12225c28e8Piotr Jasiukajtis};
12325c28e8Piotr Jasiukajtis
12425c28e8Piotr Jasiukajtisstatic const long double
12525c28e8Piotr Jasiukajtis#if defined(__x86)
12625c28e8Piotr JasiukajtisGZ1_h 	=  0.938204627909682449364570100414084663498215377L,
12725c28e8Piotr JasiukajtisGZ1_l   =  4.518346116624229420055327632718530617227944106e-20L,
12825c28e8Piotr JasiukajtisGZ2_h 	=  0.885603194410888700264725126309883762587560340L,
12925c28e8Piotr JasiukajtisGZ2_l   =  1.409077427270497062039119290776508217077297169e-20L,
13025c28e8Piotr JasiukajtisGZ3_h 	=  0.936781411463652321613537060640553022494714241L,
13125c28e8Piotr JasiukajtisGZ3_l   =  5.309836440284827247897772963887219035221996813e-21L,
13225c28e8Piotr Jasiukajtis#else
13325c28e8Piotr JasiukajtisGZ1_h 	=  0.938204627909682449409753561580326910854647031L,
13425c28e8Piotr JasiukajtisGZ1_l   =  4.684412162199460089642452580902345976446297037e-35L,
13525c28e8Piotr JasiukajtisGZ2_h 	=  0.885603194410888700278815900582588658192658794L,
13625c28e8Piotr JasiukajtisGZ2_l   =  7.501529273890253789219935569758713534641074860e-35L,
13725c28e8Piotr JasiukajtisGZ3_h 	=  0.936781411463652321618846897080837818855399840L,
13825c28e8Piotr JasiukajtisGZ3_l   =  3.088721217404784363585591914529361687403776917e-35L,
13925c28e8Piotr Jasiukajtis#endif
14025c28e8Piotr JasiukajtisTZ1	= -0.3517214357852935791015625L,
14125c28e8Piotr JasiukajtisTZ3	=  0.280530631542205810546875L;
14225c28e8Piotr Jasiukajtis/* INDENT ON */
14325c28e8Piotr Jasiukajtis
14425c28e8Piotr Jasiukajtis/* INDENT OFF */
14525c28e8Piotr Jasiukajtis/*
14625c28e8Piotr Jasiukajtis * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845]
14725c28e8Piotr Jasiukajtis * ...assume yh got 53 or 24(i386) significant bits
14825c28e8Piotr Jasiukajtis */
14925c28e8Piotr Jasiukajtis/* INDENT ON */
15025c28e8Piotr Jasiukajtisstatic struct LDouble
15125c28e8Piotr JasiukajtisGT1(long double yh, long double yl) {
15225c28e8Piotr Jasiukajtis	long double t3, t4, y;
15325c28e8Piotr Jasiukajtis	int i;
15425c28e8Piotr Jasiukajtis	struct LDouble r;
15525c28e8Piotr Jasiukajtis
15625c28e8Piotr Jasiukajtis	y = yh + yl;
15725c28e8Piotr Jasiukajtis	for (t4 = Q1[8], t3 = P1[8] + y * P1[9], i = 7; i >= 0; i--) {
15825c28e8Piotr Jasiukajtis		t4 = t4 * y + Q1[i];
15925c28e8Piotr Jasiukajtis		t3 = t3 * y + P1[i];
16025c28e8Piotr Jasiukajtis	}
16125c28e8Piotr Jasiukajtis	t3 = (y * y) * t3 / t4;
16225c28e8Piotr Jasiukajtis	t3 += (TZ1 * yl + GZ1_l);
16325c28e8Piotr Jasiukajtis	t4 = TZ1 * yh;
16425c28e8Piotr Jasiukajtis	r.h = CHOPPED((t4 + GZ1_h + t3));
16525c28e8Piotr Jasiukajtis	t3 += (t4 - (r.h - GZ1_h));
16625c28e8Piotr Jasiukajtis	r.l = t3;
16725c28e8Piotr Jasiukajtis	return (r);
16825c28e8Piotr Jasiukajtis}
16925c28e8Piotr Jasiukajtis
17025c28e8Piotr Jasiukajtis/* INDENT OFF */
17125c28e8Piotr Jasiukajtis/*
17225c28e8Piotr Jasiukajtis * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374]
17325c28e8Piotr Jasiukajtis * ...assume yh got 53 significant bits
17425c28e8Piotr Jasiukajtis */
17525c28e8Piotr Jasiukajtis/* INDENT ON */
17625c28e8Piotr Jasiukajtisstatic struct LDouble
17725c28e8Piotr JasiukajtisGT2(long double yh, long double yl) {
17825c28e8Piotr Jasiukajtis	long double t3, t4, y;
17925c28e8Piotr Jasiukajtis	int i;
18025c28e8Piotr Jasiukajtis	struct LDouble r;
18125c28e8Piotr Jasiukajtis
18225c28e8Piotr Jasiukajtis	y = yh + yl;
18325c28e8Piotr Jasiukajtis	for (t4 = Q2[9], t3 = P2[9], i = 8; i >= 0; i--) {
18425c28e8Piotr Jasiukajtis		t4 = t4 * y + Q2[i];
18525c28e8Piotr Jasiukajtis		t3 = t3 * y + P2[i];
18625c28e8Piotr Jasiukajtis	}
18725c28e8Piotr Jasiukajtis	t3 = GZ2_l + (y * y) * t3 / t4;
18825c28e8Piotr Jasiukajtis	r.h = CHOPPED((GZ2_h + t3));
18925c28e8Piotr Jasiukajtis	r.l = t3 - (r.h - GZ2_h);
19025c28e8Piotr Jasiukajtis	return (r);
19125c28e8Piotr Jasiukajtis}
19225c28e8Piotr Jasiukajtis
19325c28e8Piotr Jasiukajtis/* INDENT OFF */
19425c28e8Piotr Jasiukajtis/*
19525c28e8Piotr Jasiukajtis * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000]
19625c28e8Piotr Jasiukajtis * ...assume yh got 53 significant bits
19725c28e8Piotr Jasiukajtis */
19825c28e8Piotr Jasiukajtis/* INDENT ON */
19925c28e8Piotr Jasiukajtisstatic struct LDouble
20025c28e8Piotr JasiukajtisGT3(long double yh, long double yl) {
20125c28e8Piotr Jasiukajtis	long double t3, t4, y;
20225c28e8Piotr Jasiukajtis	int i;
20325c28e8Piotr Jasiukajtis	struct LDouble r;
20425c28e8Piotr Jasiukajtis
20525c28e8Piotr Jasiukajtis	y = yh + yl;
20625c28e8Piotr Jasiukajtis	for (t4 = Q3[9], t3 = P3[9], i = 8; i >= 0; i--) {
20725c28e8Piotr Jasiukajtis		t4 = t4 * y + Q3[i];
20825c28e8Piotr Jasiukajtis		t3 = t3 * y + P3[i];
20925c28e8Piotr Jasiukajtis	}
21025c28e8Piotr Jasiukajtis	t3 = (y * y) * t3 / t4;
21125c28e8Piotr Jasiukajtis	t3 += (TZ3 * yl + GZ3_l);
21225c28e8Piotr Jasiukajtis	t4 = TZ3 * yh;
21325c28e8Piotr Jasiukajtis	r.h = CHOPPED((t4 + GZ3_h + t3));
21425c28e8Piotr Jasiukajtis	t3 += (t4 - (r.h - GZ3_h));
21525c28e8Piotr Jasiukajtis	r.l = t3;
21625c28e8Piotr Jasiukajtis	return (r);
21725c28e8Piotr Jasiukajtis}
21825c28e8Piotr Jasiukajtis
21925c28e8Piotr Jasiukajtis/* INDENT OFF */
22025c28e8Piotr Jasiukajtis/* Hex value of GP[0] shoule be 3FB55555 55555555 */
22125c28e8Piotr Jasiukajtisstatic const long double GP[] = {
22225c28e8Piotr Jasiukajtis	+0.083333333333333333333333333333333172839171301L,
22325c28e8Piotr Jasiukajtis	-2.77777777777777777777777777492501211999399424104e-0003L,
22425c28e8Piotr Jasiukajtis	+7.93650793650793650793635650541638236350020883243e-0004L,
22525c28e8Piotr Jasiukajtis	-5.95238095238095238057299772679324503339241961704e-0004L,
22625c28e8Piotr Jasiukajtis	+8.41750841750841696138422987977683524926142600321e-0004L,
22725c28e8Piotr Jasiukajtis	-1.91752691752686682825032547823699662178842123308e-0003L,
22825c28e8Piotr Jasiukajtis	+6.41025641022403480921891559356473451161279359322e-0003L,
22925c28e8Piotr Jasiukajtis	-2.95506535798414019189819587455577003732808185071e-0002L,
23025c28e8Piotr Jasiukajtis	+1.79644367229970031486079180060923073476568732136e-0001L,
23125c28e8Piotr Jasiukajtis	-1.39243086487274662174562872567057200255649290646e+0000L,
23225c28e8Piotr Jasiukajtis	+1.34025874044417962188677816477842265259608269775e+0001L,
23325c28e8Piotr Jasiukajtis	-1.56803713480127469414495545399982508700748274318e+0002L,
23425c28e8Piotr Jasiukajtis	+2.18739841656201561694927630335099313968924493891e+0003L,
23525c28e8Piotr Jasiukajtis	-3.55249848644100338419187038090925410976237921269e+0004L,
23625c28e8Piotr Jasiukajtis	+6.43464880437835286216768959439484376449179576452e+0005L,
23725c28e8Piotr Jasiukajtis	-1.20459154385577014992600342782821389605893904624e+0007L,
23825c28e8Piotr Jasiukajtis	+2.09263249637351298563934942349749718491071093210e+0008L,
23925c28e8Piotr Jasiukajtis	-2.96247483183169219343745316433899599834685703457e+0009L,
24025c28e8Piotr Jasiukajtis	+2.88984933605896033154727626086506756972327292981e+0010L,
24125c28e8Piotr Jasiukajtis	-1.40960434146030007732838382416230610302678063984e+0011L,	/* 19 */
24225c28e8Piotr Jasiukajtis};
24325c28e8Piotr Jasiukajtis
24425c28e8Piotr Jasiukajtisstatic const long double T3[] = {
24525c28e8Piotr Jasiukajtis	+0.666666666666666666666666666666666634567834260213L,	/* T3[0] */
24625c28e8Piotr Jasiukajtis	+0.400000000000000000000000000040853636176634934140L,	/* T3[1] */
24725c28e8Piotr Jasiukajtis	+0.285714285714285714285696975252753987869020263448L,	/* T3[2] */
24825c28e8Piotr Jasiukajtis	+0.222222222222222225593221101192317258554772129875L,	/* T3[3] */
24925c28e8Piotr Jasiukajtis	+0.181818181817850192105847183461778186703779262916L,	/* T3[4] */
25025c28e8Piotr Jasiukajtis	+0.153846169861348633757101285952333369222567014596L,	/* T3[5] */
25125c28e8Piotr Jasiukajtis	+0.133033462889260193922261296772841229985047571265L,	/* T3[6] */
25225c28e8Piotr Jasiukajtis};
25325c28e8Piotr Jasiukajtis
25425c28e8Piotr Jasiukajtisstatic const long double c[] = {
25525c28e8Piotr Jasiukajtis0.0L,
25625c28e8Piotr Jasiukajtis1.0L,
25725c28e8Piotr Jasiukajtis2.0L,
25825c28e8Piotr Jasiukajtis0.5L,
25925c28e8Piotr Jasiukajtis1.0e-4930L,							/* tiny */
26025c28e8Piotr Jasiukajtis4.18937683105468750000e-01L,					/* hln2pim1_h */
26125c28e8Piotr Jasiukajtis8.50099203991780329736405617639861397473637783412817152e-07L,	/* hln2pim1_l */
26225c28e8Piotr Jasiukajtis0.418938533204672741780329736405617639861397473637783412817152L, /* hln2pim1 */
26325c28e8Piotr Jasiukajtis2.16608493865351192653179168701171875e-02L,			/* ln2_32hi */
26425c28e8Piotr Jasiukajtis5.96317165397058692545083025235937919875797669127130e-12L,	/* ln2_32lo */
26525c28e8Piotr Jasiukajtis46.16624130844682903551758979206054839765267053289554989233L,	/* invln2_32 */
26625c28e8Piotr Jasiukajtis#if defined(__x86)
26725c28e8Piotr Jasiukajtis1.7555483429044629170023839037639845628291e+03L,		/* overflow */
26825c28e8Piotr Jasiukajtis#else
26925c28e8Piotr Jasiukajtis1.7555483429044629170038892160702032034177e+03L,		/* overflow */
27025c28e8Piotr Jasiukajtis#endif
27125c28e8Piotr Jasiukajtis};
27225c28e8Piotr Jasiukajtis
27325c28e8Piotr Jasiukajtis#define	zero		c[0]
27425c28e8Piotr Jasiukajtis#define	one		c[1]
27525c28e8Piotr Jasiukajtis#define	two		c[2]
27625c28e8Piotr Jasiukajtis#define	half		c[3]
27725c28e8Piotr Jasiukajtis#define	tiny		c[4]
27825c28e8Piotr Jasiukajtis#define	hln2pim1_h	c[5]
27925c28e8Piotr Jasiukajtis#define	hln2pim1_l	c[6]
28025c28e8Piotr Jasiukajtis#define	hln2pim1	c[7]
28125c28e8Piotr Jasiukajtis#define	ln2_32hi	c[8]
28225c28e8Piotr Jasiukajtis#define	ln2_32lo	c[9]
28325c28e8Piotr Jasiukajtis#define	invln2_32	c[10]
28425c28e8Piotr Jasiukajtis#define	overflow	c[11]
28525c28e8Piotr Jasiukajtis
28625c28e8Piotr Jasiukajtis/*
28725c28e8Piotr Jasiukajtis * |exp(r) - (1+r+Et0*r^2+...+Et10*r^12)| <= 2^(-128.88) for |r|<=ln2/64
28825c28e8Piotr Jasiukajtis */
28925c28e8Piotr Jasiukajtisstatic const long double Et[] = {
29025c28e8Piotr Jasiukajtis	+5.0000000000000000000e-1L,
29125c28e8Piotr Jasiukajtis	+1.66666666666666666666666666666828835166292152466e-0001L,
29225c28e8Piotr Jasiukajtis	+4.16666666666666666666666666666693398646592712189e-0002L,
29325c28e8Piotr Jasiukajtis	+8.33333333333333333333331748774512601775591115951e-0003L,
29425c28e8Piotr Jasiukajtis	+1.38888888888888888888888845356011511394764753997e-0003L,
29525c28e8Piotr Jasiukajtis	+1.98412698412698413237140350092993252684198882102e-0004L,
29625c28e8Piotr Jasiukajtis	+2.48015873015873016080222025357442659895814371694e-0005L,
29725c28e8Piotr Jasiukajtis	+2.75573192239028921114572986441972140933432317798e-0006L,
29825c28e8Piotr Jasiukajtis	+2.75573192239448470555548102895526369739856219317e-0007L,
29925c28e8Piotr Jasiukajtis	+2.50521677867683935940853997995937600214167232477e-0008L,
30025c28e8Piotr Jasiukajtis	+2.08767928899010367374984448513685566514152147362e-0009L,
30125c28e8Piotr Jasiukajtis};
30225c28e8Piotr Jasiukajtis
30325c28e8Piotr Jasiukajtis/*
30425c28e8Piotr Jasiukajtis * long double precision coefficients for computing log(x)-1 in tgamma.
30525c28e8Piotr Jasiukajtis *  See "algorithm" for details
30625c28e8Piotr Jasiukajtis *
30725c28e8Piotr Jasiukajtis *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
30825c28e8Piotr Jasiukajtis *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
30925c28e8Piotr Jasiukajtis *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
31025c28e8Piotr Jasiukajtis *       T2(j) = T2[2j,2j+1] = log(z[j]),
31125c28e8Piotr Jasiukajtis *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + T3[2]s^7 + ... + T3[6]s^15
31225c28e8Piotr Jasiukajtis *  Note
31325c28e8Piotr Jasiukajtis *  (1) the leading entries are truncated to 24 binary point.
31425c28e8Piotr Jasiukajtis *  (2) Remez error for T3(s) is bounded by 2**(-136.54)
31525c28e8Piotr Jasiukajtis */
31625c28e8Piotr Jasiukajtisstatic const long double T1[] = {
31725c28e8Piotr Jasiukajtis-1.000000000000000000000000000000000000000000e+00L,
31825c28e8Piotr Jasiukajtis	+0.000000000000000000000000000000000000000000e+00L,
31925c28e8Piotr Jasiukajtis-3.068528175354003906250000000000000000000000e-01L,
32025c28e8Piotr Jasiukajtis-1.904654299957767878541823431924500011926579e-09L,
32125c28e8Piotr Jasiukajtis	+3.862943053245544433593750000000000000000000e-01L,
32225c28e8Piotr Jasiukajtis	+5.579533617547508924291635313615100141107647e-08L,
32325c28e8Piotr Jasiukajtis	+1.079441487789154052734375000000000000000000e+00L,
32425c28e8Piotr Jasiukajtis	+5.389068187551732136437452970422650211661470e-08L,
32525c28e8Piotr Jasiukajtis	+1.772588670253753662109375000000000000000000e+00L,
32625c28e8Piotr Jasiukajtis	+5.198602757555955348583270627230200282215294e-08L,
32725c28e8Piotr Jasiukajtis	+2.465735852718353271484375000000000000000000e+00L,
32825c28e8Piotr Jasiukajtis	+5.008137327560178560729088284037750352769117e-08L,
32925c28e8Piotr Jasiukajtis	+3.158883035182952880859375000000000000000000e+00L,
33025c28e8Piotr Jasiukajtis	+4.817671897564401772874905940845299849351090e-08L,
33125c28e8Piotr Jasiukajtis	+3.852030217647552490234375000000000000000000e+00L,
33225c28e8Piotr Jasiukajtis	+4.627206467568624985020723597652849919904913e-08L,
33325c28e8Piotr Jasiukajtis	+4.545177400112152099609375000000000000000000e+00L,
33425c28e8Piotr Jasiukajtis	+4.436741037572848197166541254460399990458737e-08L,
33525c28e8Piotr Jasiukajtis	+5.238324582576751708984375000000000000000000e+00L,
33625c28e8Piotr Jasiukajtis	+4.246275607577071409312358911267950061012560e-08L,
33725c28e8Piotr Jasiukajtis	+5.931471765041351318359375000000000000000000e+00L,
33825c28e8Piotr Jasiukajtis	+4.055810177581294621458176568075500131566384e-08L,
33925c28e8Piotr Jasiukajtis};
34025c28e8Piotr Jasiukajtis
34125c28e8Piotr Jasiukajtis/*
34225c28e8Piotr Jasiukajtis * T2[2i,2i+1] = log(1+i/64+1/128)
34325c28e8Piotr Jasiukajtis */
34425c28e8Piotr Jasiukajtisstatic const long double T2[] = {
34525c28e8Piotr Jasiukajtis	+7.7821016311645507812500000000000000000000e-03L,
34625c28e8Piotr Jasiukajtis	+3.8810890398166212900061136763678127453570e-08L,
34725c28e8Piotr Jasiukajtis	+2.3167014122009277343750000000000000000000e-02L,
34825c28e8Piotr Jasiukajtis	+4.5159525100885049160962289916579411752759e-08L,
34925c28e8Piotr Jasiukajtis	+3.8318812847137451171875000000000000000000e-02L,
35025c28e8Piotr Jasiukajtis	+5.1454999148021880325123797290345960518164e-08L,
35125c28e8Piotr Jasiukajtis	+5.3244471549987792968750000000000000000000e-02L,
35225c28e8Piotr Jasiukajtis	+4.2968824489897120193786528776939573415076e-08L,
35325c28e8Piotr Jasiukajtis	+6.7950606346130371093750000000000000000000e-02L,
35425c28e8Piotr Jasiukajtis	+5.5562377378300815277772629414034632394030e-08L,
35525c28e8Piotr Jasiukajtis	+8.2443654537200927734375000000000000000000e-02L,
35625c28e8Piotr Jasiukajtis	+1.4673873663533785068668307805914095366600e-08L,
35725c28e8Piotr Jasiukajtis	+9.6729576587677001953125000000000000000000e-02L,
35825c28e8Piotr Jasiukajtis	+4.9870874110342446056487463437015041543346e-08L,
35925c28e8Piotr Jasiukajtis	+1.1081433296203613281250000000000000000000e-01L,
36025c28e8Piotr Jasiukajtis	+3.3378253981382306169323211928098474801099e-08L,
36125c28e8Piotr Jasiukajtis	+1.2470346689224243164062500000000000000000e-01L,
36225c28e8Piotr Jasiukajtis	+1.1608714804222781515380863268491613205318e-08L,
36325c28e8Piotr Jasiukajtis	+1.3840228319168090820312500000000000000000e-01L,
36425c28e8Piotr Jasiukajtis	+3.9667438227482200873601649187393160823607e-08L,
36525c28e8Piotr Jasiukajtis	+1.5191602706909179687500000000000000000000e-01L,
36625c28e8Piotr Jasiukajtis	+1.4956750178196803424896884511327584958252e-08L,
36725c28e8Piotr Jasiukajtis	+1.6524952650070190429687500000000000000000e-01L,
36825c28e8Piotr Jasiukajtis	+4.6394605258578736449277240313729237989366e-08L,
36925c28e8Piotr Jasiukajtis	+1.7840760946273803710937500000000000000000e-01L,
37025c28e8Piotr Jasiukajtis	+4.8010080260010025241510941968354682199540e-08L,
37125c28e8Piotr Jasiukajtis	+1.9139480590820312500000000000000000000000e-01L,
37225c28e8Piotr Jasiukajtis	+4.7091426329609298807561308873447039132856e-08L,
37325c28e8Piotr Jasiukajtis	+2.0421552658081054687500000000000000000000e-01L,
37425c28e8Piotr Jasiukajtis	+1.4847880344628820386196239272213742113867e-08L,
37525c28e8Piotr Jasiukajtis	+2.1687388420104980468750000000000000000000e-01L,
37625c28e8Piotr Jasiukajtis	+5.4099564554931589525744347498478964801484e-08L,
37725c28e8Piotr Jasiukajtis	+2.2937405109405517578125000000000000000000e-01L,
37825c28e8Piotr Jasiukajtis	+4.9970790654210230725046139871550961365282e-08L,
37925c28e8Piotr Jasiukajtis	+2.4171990156173706054687500000000000000000e-01L,
38025c28e8Piotr Jasiukajtis	+3.5325408107597432515913513900103385655073e-08L,
38125c28e8Piotr Jasiukajtis	+2.5391519069671630859375000000000000000000e-01L,
38225c28e8Piotr Jasiukajtis	+1.9284247135543573297906606667466299224747e-08L,
38325c28e8Piotr Jasiukajtis	+2.6596349477767944335937500000000000000000e-01L,
38425c28e8Piotr Jasiukajtis	+5.3719458497979750926537543389268821141517e-08L,
38525c28e8Piotr Jasiukajtis	+2.7786844968795776367187500000000000000000e-01L,
38625c28e8Piotr Jasiukajtis	+1.3154985425144750329234012330820349974537e-09L,
38725c28e8Piotr Jasiukajtis	+2.8963327407836914062500000000000000000000e-01L,
38825c28e8Piotr Jasiukajtis	+1.8504673536253893055525668970003860369760e-08L,
38925c28e8Piotr Jasiukajtis	+3.0126130580902099609375000000000000000000e-01L,
39025c28e8Piotr Jasiukajtis	+2.4769140784919125538233755492657352680723e-08L,
39125c28e8Piotr Jasiukajtis	+3.1275570392608642578125000000000000000000e-01L,
39225c28e8Piotr Jasiukajtis	+6.0778104626049965596883190321597861455475e-09L,
39325c28e8Piotr Jasiukajtis	+3.2411944866180419921875000000000000000000e-01L,
39425c28e8Piotr Jasiukajtis	+1.9992407776871920760434987352182336158873e-08L,
39525c28e8Piotr Jasiukajtis	+3.3535552024841308593750000000000000000000e-01L,
39625c28e8Piotr Jasiukajtis	+2.1672724744319679579814166199074433006807e-08L,
39725c28e8Piotr Jasiukajtis	+3.4646672010421752929687500000000000000000e-01L,
39825c28e8Piotr Jasiukajtis	+4.7241991051621587188425772950711830538414e-08L,
39925c28e8Piotr Jasiukajtis	+3.5745584964752197265625000000000000000000e-01L,
40025c28e8Piotr Jasiukajtis	+3.9274281801569759490140904474434669956562e-08L,
40125c28e8Piotr Jasiukajtis	+3.6832553148269653320312500000000000000000e-01L,
40225c28e8Piotr Jasiukajtis	+2.9676011119845105154050398826897178765758e-08L,
40325c28e8Piotr Jasiukajtis	+3.7907832860946655273437500000000000000000e-01L,
40425c28e8Piotr Jasiukajtis	+2.4325502905656478345631019858881408009210e-08L,
40525c28e8Piotr Jasiukajtis	+3.8971674442291259765625000000000000000000e-01L,
40625c28e8Piotr Jasiukajtis	+6.7171126157142136040035208670510556529487e-09L,
40725c28e8Piotr Jasiukajtis	+4.0024316310882568359375000000000000000000e-01L,
40825c28e8Piotr Jasiukajtis	+1.0181870233355751019951311700799406124957e-09L,
40925c28e8Piotr Jasiukajtis	+4.1065990924835205078125000000000000000000e-01L,
41025c28e8Piotr Jasiukajtis	+1.5736916335153056203175822787661567534220e-08L,
41125c28e8Piotr Jasiukajtis	+4.2096924781799316406250000000000000000000e-01L,
41225c28e8Piotr Jasiukajtis	+4.6826136472066367161506795972449857268707e-08L,
41325c28e8Piotr Jasiukajtis	+4.3117344379425048828125000000000000000000e-01L,
41425c28e8Piotr Jasiukajtis	+2.1024120852577922478955594998480144051225e-08L,
41525c28e8Piotr Jasiukajtis	+4.4127452373504638671875000000000000000000e-01L,
41625c28e8Piotr Jasiukajtis	+3.7069828842770746441661301225362605528786e-08L,
41725c28e8Piotr Jasiukajtis	+4.5127463340759277343750000000000000000000e-01L,
41825c28e8Piotr Jasiukajtis	+1.0731865811707192383079012478685922879010e-08L,
41925c28e8Piotr Jasiukajtis	+4.6117568016052246093750000000000000000000e-01L,
42025c28e8Piotr Jasiukajtis	+3.4961647705430499925597855358603099030515e-08L,
42125c28e8Piotr Jasiukajtis	+4.7097969055175781250000000000000000000000e-01L,
42225c28e8Piotr Jasiukajtis	+2.4667033200046897856056359251373510964634e-08L,
42325c28e8Piotr Jasiukajtis	+4.8068851232528686523437500000000000000000e-01L,
42425c28e8Piotr Jasiukajtis	+1.7020465042442243455448011551208861216878e-08L,
42525c28e8Piotr Jasiukajtis	+4.9030393362045288085937500000000000000000e-01L,
42625c28e8Piotr Jasiukajtis	+5.4424740957290971159645746860530583309571e-08L,
42725c28e8Piotr Jasiukajtis	+4.9982786178588867187500000000000000000000e-01L,
42825c28e8Piotr Jasiukajtis	+7.7705606579463314152470441415126573566105e-09L,
42925c28e8Piotr Jasiukajtis	+5.0926184654235839843750000000000000000000e-01L,
43025c28e8Piotr Jasiukajtis	+5.5247449548366574919228323824878565745713e-08L,
43125c28e8Piotr Jasiukajtis	+5.1860773563385009765625000000000000000000e-01L,
43225c28e8Piotr Jasiukajtis	+2.8574195534496726996364798698556235730848e-08L,
43325c28e8Piotr Jasiukajtis	+5.2786707878112792968750000000000000000000e-01L,
43425c28e8Piotr Jasiukajtis	+1.0839714455426392217778300963558522088193e-08L,
43525c28e8Piotr Jasiukajtis	+5.3704142570495605468750000000000000000000e-01L,
43625c28e8Piotr Jasiukajtis	+4.0191927599879229244153832299023744345999e-08L,
43725c28e8Piotr Jasiukajtis	+5.4613238573074340820312500000000000000000e-01L,
43825c28e8Piotr Jasiukajtis	+5.1867392242179272209231209163864971792889e-08L,
43925c28e8Piotr Jasiukajtis	+5.5514144897460937500000000000000000000000e-01L,
44025c28e8Piotr Jasiukajtis	+5.8565892217715480359515904050170125743178e-08L,
44125c28e8Piotr Jasiukajtis	+5.6407010555267333984375000000000000000000e-01L,
44225c28e8Piotr Jasiukajtis	+3.2732129626227634290090190711817681692354e-08L,
44325c28e8Piotr Jasiukajtis	+5.7291972637176513671875000000000000000000e-01L,
44425c28e8Piotr Jasiukajtis	+2.7190020372374006726626261068626400393936e-08L,
44525c28e8Piotr Jasiukajtis	+5.8169168233871459960937500000000000000000e-01L,
44625c28e8Piotr Jasiukajtis	+5.7295907882911235753725372340709967597394e-08L,
44725c28e8Piotr Jasiukajtis	+5.9038740396499633789062500000000000000000e-01L,
44825c28e8Piotr Jasiukajtis	+4.2637180036751291708123598757577783615014e-08L,
44925c28e8Piotr Jasiukajtis	+5.9900814294815063476562500000000000000000e-01L,
45025c28e8Piotr Jasiukajtis	+4.6697932764615975024461651502060474048774e-08L,
45125c28e8Piotr Jasiukajtis	+6.0755521059036254882812500000000000000000e-01L,
45225c28e8Piotr Jasiukajtis	+3.9634179246672960152791125371893149820625e-08L,
45325c28e8Piotr Jasiukajtis	+6.1602985858917236328125000000000000000000e-01L,
45425c28e8Piotr Jasiukajtis	+1.8626341656366315928196700650292529688219e-08L,
45525c28e8Piotr Jasiukajtis	+6.2443327903747558593750000000000000000000e-01L,
45625c28e8Piotr Jasiukajtis	+8.9744179151050387440546731199093039879228e-09L,
45725c28e8Piotr Jasiukajtis	+6.3276666402816772460937500000000000000000e-01L,
45825c28e8Piotr Jasiukajtis	+5.5428701049364114685035797584887586099726e-09L,
45925c28e8Piotr Jasiukajtis	+6.4103114604949951171875000000000000000000e-01L,
46025c28e8Piotr Jasiukajtis	+3.3371431779336851334405392546708949047361e-08L,
46125c28e8Piotr Jasiukajtis	+6.4922791719436645507812500000000000000000e-01L,
46225c28e8Piotr Jasiukajtis	+2.9430743363812714969905311122271269100885e-08L,
46325c28e8Piotr Jasiukajtis	+6.5735805034637451171875000000000000000000e-01L,
46425c28e8Piotr Jasiukajtis	+2.2361985518423140023245936165514147093250e-08L,
46525c28e8Piotr Jasiukajtis	+6.6542261838912963867187500000000000000000e-01L,
46625c28e8Piotr Jasiukajtis	+1.4155960810278217610006660181148303091649e-08L,
46725c28e8Piotr Jasiukajtis	+6.7342263460159301757812500000000000000000e-01L,
46825c28e8Piotr Jasiukajtis	+4.0610573702719835388801017264750843477878e-08L,
46925c28e8Piotr Jasiukajtis	+6.8135917186737060546875000000000000000000e-01L,
47025c28e8Piotr Jasiukajtis	+5.2940532463479321559568089441735584156689e-08L,
47125c28e8Piotr Jasiukajtis	+6.8923324346542358398437500000000000000000e-01L,
47225c28e8Piotr Jasiukajtis	+3.7773385396340539337814603903232796216537e-08L,
47325c28e8Piotr Jasiukajtis};
47425c28e8Piotr Jasiukajtis
47525c28e8Piotr Jasiukajtis/*
47625c28e8Piotr Jasiukajtis * S[j],S_trail[j] = 2**(j/32.) for the final computation of exp(t+w)
47725c28e8Piotr Jasiukajtis */
47825c28e8Piotr Jasiukajtisstatic const long double S[] = {
47925c28e8Piotr Jasiukajtis#if defined(__x86)
48025c28e8Piotr Jasiukajtis	+1.0000000000000000000000000e+00L,
48125c28e8Piotr Jasiukajtis	+1.0218971486541166782081522e+00L,
48225c28e8Piotr Jasiukajtis	+1.0442737824274138402382006e+00L,
48325c28e8Piotr Jasiukajtis	+1.0671404006768236181297224e+00L,
48425c28e8Piotr Jasiukajtis	+1.0905077326652576591003302e+00L,
48525c28e8Piotr Jasiukajtis	+1.1143867425958925362894369e+00L,
48625c28e8Piotr Jasiukajtis	+1.1387886347566916536971221e+00L,
48725c28e8Piotr Jasiukajtis	+1.1637248587775775137938619e+00L,
48825c28e8Piotr Jasiukajtis	+1.1892071150027210666875674e+00L,
48925c28e8Piotr Jasiukajtis	+1.2152473599804688780476325e+00L,
49025c28e8Piotr Jasiukajtis	+1.2418578120734840485256747e+00L,
49125c28e8Piotr Jasiukajtis	+1.2690509571917332224885722e+00L,
49225c28e8Piotr Jasiukajtis	+1.2968395546510096659215822e+00L,
49325c28e8Piotr Jasiukajtis	+1.3252366431597412945939118e+00L,
49425c28e8Piotr Jasiukajtis	+1.3542555469368927282668852e+00L,
49525c28e8Piotr Jasiukajtis	+1.3839098819638319548151403e+00L,
49625c28e8Piotr Jasiukajtis	+1.4142135623730950487637881e+00L,
49725c28e8Piotr Jasiukajtis	+1.4451808069770466200253470e+00L,
49825c28e8Piotr Jasiukajtis	+1.4768261459394993113155431e+00L,
49925c28e8Piotr Jasiukajtis	+1.5091644275934227397133885e+00L,
50025c28e8Piotr Jasiukajtis	+1.5422108254079408235859630e+00L,
50125c28e8Piotr Jasiukajtis	+1.5759808451078864864006862e+00L,
50225c28e8Piotr Jasiukajtis	+1.6104903319492543080837174e+00L,
50325c28e8Piotr Jasiukajtis	+1.6457554781539648445110730e+00L,
50425c28e8Piotr Jasiukajtis	+1.6817928305074290860378350e+00L,
50525c28e8Piotr Jasiukajtis	+1.7186192981224779156032914e+00L,
50625c28e8Piotr Jasiukajtis	+1.7562521603732994831094730e+00L,
50725c28e8Piotr Jasiukajtis	+1.7947090750031071864148413e+00L,
50825c28e8Piotr Jasiukajtis	+1.8340080864093424633989166e+00L,
50925c28e8Piotr Jasiukajtis	+1.8741676341102999013002103e+00L,
51025c28e8Piotr Jasiukajtis	+1.9152065613971472938202589e+00L,
51125c28e8Piotr Jasiukajtis	+1.9571441241754002689657438e+00L,
51225c28e8Piotr Jasiukajtis#else
51325c28e8Piotr Jasiukajtis	+1.00000000000000000000000000000000000e+00L,
51425c28e8Piotr Jasiukajtis	+1.02189714865411667823448013478329942e+00L,
51525c28e8Piotr Jasiukajtis	+1.04427378242741384032196647873992910e+00L,
51625c28e8Piotr Jasiukajtis	+1.06714040067682361816952112099280918e+00L,
51725c28e8Piotr Jasiukajtis	+1.09050773266525765920701065576070789e+00L,
51825c28e8Piotr Jasiukajtis	+1.11438674259589253630881295691960313e+00L,
51925c28e8Piotr Jasiukajtis	+1.13878863475669165370383028384151134e+00L,
52025c28e8Piotr Jasiukajtis	+1.16372485877757751381357359909218536e+00L,
52125c28e8Piotr Jasiukajtis	+1.18920711500272106671749997056047593e+00L,
52225c28e8Piotr Jasiukajtis	+1.21524735998046887811652025133879836e+00L,
52325c28e8Piotr Jasiukajtis	+1.24185781207348404859367746872659561e+00L,
52425c28e8Piotr Jasiukajtis	+1.26905095719173322255441908103233805e+00L,
52525c28e8Piotr Jasiukajtis	+1.29683955465100966593375411779245118e+00L,
52625c28e8Piotr Jasiukajtis	+1.32523664315974129462953709549872168e+00L,
52725c28e8Piotr Jasiukajtis	+1.35425554693689272829801474014070273e+00L,
52825c28e8Piotr Jasiukajtis	+1.38390988196383195487265952726519287e+00L,
52925c28e8Piotr Jasiukajtis	+1.41421356237309504880168872420969798e+00L,
53025c28e8Piotr Jasiukajtis	+1.44518080697704662003700624147167095e+00L,
53125c28e8Piotr Jasiukajtis	+1.47682614593949931138690748037404985e+00L,
53225c28e8Piotr Jasiukajtis	+1.50916442759342273976601955103319352e+00L,
53325c28e8Piotr Jasiukajtis	+1.54221082540794082361229186209073479e+00L,
53425c28e8Piotr Jasiukajtis	+1.57598084510788648645527016018190504e+00L,
53525c28e8Piotr Jasiukajtis	+1.61049033194925430817952066735740067e+00L,
53625c28e8Piotr Jasiukajtis	+1.64575547815396484451875672472582254e+00L,
53725c28e8Piotr Jasiukajtis	+1.68179283050742908606225095246642969e+00L,
53825c28e8Piotr Jasiukajtis	+1.71861929812247791562934437645631244e+00L,
53925c28e8Piotr Jasiukajtis	+1.75625216037329948311216061937531314e+00L,
54025c28e8Piotr Jasiukajtis	+1.79470907500310718642770324212778174e+00L,
54125c28e8Piotr Jasiukajtis	+1.83400808640934246348708318958828892e+00L,
54225c28e8Piotr Jasiukajtis	+1.87416763411029990132999894995444645e+00L,
54325c28e8Piotr Jasiukajtis	+1.91520656139714729387261127029583086e+00L,
54425c28e8Piotr Jasiukajtis	+1.95714412417540026901832225162687149e+00L,
54525c28e8Piotr Jasiukajtis#endif
54625c28e8Piotr Jasiukajtis};
54725c28e8Piotr Jasiukajtisstatic const long double S_trail[] = {
54825c28e8Piotr Jasiukajtis#if defined(__x86)
54925c28e8Piotr Jasiukajtis	+0.0000000000000000000000000e+00L,
55025c28e8Piotr Jasiukajtis	+2.6327965667180882569382524e-20L,
55125c28e8Piotr Jasiukajtis	+8.3765863521895191129661899e-20L,
55225c28e8Piotr Jasiukajtis	+3.9798705777454504249209575e-20L,
55325c28e8Piotr Jasiukajtis	+1.0668046596651558640993042e-19L,
55425c28e8Piotr Jasiukajtis	+1.9376009847285360448117114e-20L,
55525c28e8Piotr Jasiukajtis	+6.7081819456112953751277576e-21L,
55625c28e8Piotr Jasiukajtis	+1.9711680502629186462729727e-20L,
55725c28e8Piotr Jasiukajtis	+2.9932584438449523689104569e-20L,
55825c28e8Piotr Jasiukajtis	+6.8887754153039109411061914e-20L,
55925c28e8Piotr Jasiukajtis	+6.8002718741225378942847820e-20L,
56025c28e8Piotr Jasiukajtis	+6.5846917376975403439742349e-20L,
56125c28e8Piotr Jasiukajtis	+1.2171958727511372194876001e-20L,
56225c28e8Piotr Jasiukajtis	+3.5625253228704087115438260e-20L,
56325c28e8Piotr Jasiukajtis	+3.1129551559077560956309179e-20L,
56425c28e8Piotr Jasiukajtis	+5.7519192396164779846216492e-20L,
56525c28e8Piotr Jasiukajtis	+3.7900651177865141593101239e-20L,
56625c28e8Piotr Jasiukajtis	+1.1659262405698741798080115e-20L,
56725c28e8Piotr Jasiukajtis	+7.1364385105284695967172478e-20L,
56825c28e8Piotr Jasiukajtis	+5.2631003710812203588788949e-20L,
56925c28e8Piotr Jasiukajtis	+2.6328853788732632868460580e-20L,
57025c28e8Piotr Jasiukajtis	+5.4583950085438242788190141e-20L,
57125c28e8Piotr Jasiukajtis	+9.5803254376938269960718656e-20L,
57225c28e8Piotr Jasiukajtis	+7.6837733983874245823512279e-21L,
57325c28e8Piotr Jasiukajtis	+2.4415965910835093824202087e-20L,
57425c28e8Piotr Jasiukajtis	+2.6052966871016580981769728e-20L,
57525c28e8Piotr Jasiukajtis	+2.6876456344632553875309579e-21L,
57625c28e8Piotr Jasiukajtis	+1.2861930155613700201703279e-20L,
57725c28e8Piotr Jasiukajtis	+8.8166633394037485606572294e-20L,
57825c28e8Piotr Jasiukajtis	+2.9788615389580190940837037e-20L,
57925c28e8Piotr Jasiukajtis	+5.2352341619805098677422139e-20L,
58025c28e8Piotr Jasiukajtis	+5.2578463064010463732242363e-20L,
58125c28e8Piotr Jasiukajtis#else
58225c28e8Piotr Jasiukajtis	+0.00000000000000000000000000000000000e+00L,
58325c28e8Piotr Jasiukajtis	+1.80506787420330954745573333054573786e-35L,
58425c28e8Piotr Jasiukajtis-9.37452029228042742195756741973083214e-35L,
58525c28e8Piotr Jasiukajtis-1.59696844729275877071290963023149997e-35L,
58625c28e8Piotr Jasiukajtis	+9.11249341012502297851168610167248666e-35L,
58725c28e8Piotr Jasiukajtis-6.50422820697854828723037477525938871e-35L,
58825c28e8Piotr Jasiukajtis-8.14846884452585113732569176748815532e-35L,
58925c28e8Piotr Jasiukajtis-5.06621457672180031337233074514290335e-35L,
59025c28e8Piotr Jasiukajtis-1.35983097468881697374987563824591912e-35L,
59125c28e8Piotr Jasiukajtis	+9.49742763556319647030771056643324660e-35L,
59225c28e8Piotr Jasiukajtis-3.28317052317699860161506596533391526e-36L,
59325c28e8Piotr Jasiukajtis-5.01723570938719041029018653045842895e-35L,
59425c28e8Piotr Jasiukajtis-2.39147479768910917162283430160264014e-35L,
59525c28e8Piotr Jasiukajtis-8.35057135763390881529889073794408385e-36L,
59625c28e8Piotr Jasiukajtis	+7.03675688907326504242173719067187644e-35L,
59725c28e8Piotr Jasiukajtis-5.18248485306464645753689301856695619e-35L,
59825c28e8Piotr Jasiukajtis	+9.42224254862183206569211673639406488e-35L,
59925c28e8Piotr Jasiukajtis-3.96750082539886230916730613021641828e-35L,
60025c28e8Piotr Jasiukajtis	+7.14352899156330061452327361509276724e-35L,
60125c28e8Piotr Jasiukajtis	+1.15987125286798512424651783410044433e-35L,
60225c28e8Piotr Jasiukajtis	+4.69693347835811549530973921320187447e-35L,
60325c28e8Piotr Jasiukajtis-3.38651317599500471079924198499981917e-35L,
60425c28e8Piotr Jasiukajtis-8.58731877429824706886865593510387445e-35L,
60525c28e8Piotr Jasiukajtis-9.60595154874935050318549936224606909e-35L,
60625c28e8Piotr Jasiukajtis	+9.60973393212801278450755869714178581e-35L,
60725c28e8Piotr Jasiukajtis	+6.37839792144002843924476144978084855e-35L,
60825c28e8Piotr Jasiukajtis	+7.79243078569586424945646112516927770e-35L,
60925c28e8Piotr Jasiukajtis	+7.36133776758845652413193083663393220e-35L,
61025c28e8Piotr Jasiukajtis-6.47299514791334723003521457561217053e-35L,
61125c28e8Piotr Jasiukajtis	+8.58747441795369869427879806229522962e-35L,
61225c28e8Piotr Jasiukajtis	+2.37181542282517483569165122830269098e-35L,
61325c28e8Piotr Jasiukajtis-3.02689168209611877300459737342190031e-37L,
61425c28e8Piotr Jasiukajtis#endif
61525c28e8Piotr Jasiukajtis};
61625c28e8Piotr Jasiukajtis/* INDENT ON */
61725c28e8Piotr Jasiukajtis
61825c28e8Piotr Jasiukajtis/* INDENT OFF */
61925c28e8Piotr Jasiukajtis/*
62025c28e8Piotr Jasiukajtis * return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
62125c28e8Piotr Jasiukajtis *     log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
62225c28e8Piotr Jasiukajtis *                = L1 + L2 + L3,
62325c28e8Piotr Jasiukajtis */
62425c28e8Piotr Jasiukajtis/* INDENT ON */
62525c28e8Piotr Jasiukajtisstatic struct LDouble
62625c28e8Piotr Jasiukajtislarge_gam(long double x, int *m) {
62725c28e8Piotr Jasiukajtis	long double z, t1, t2, t3, z2, t5, w, y, u, r, v;
62825c28e8Piotr Jasiukajtis	long double t24 = 16777216.0L, p24 = 1.0L / 16777216.0L;
62925c28e8Piotr Jasiukajtis	int n2, j2, k, ix, j, i;
63025c28e8Piotr Jasiukajtis	struct LDouble zz;
63125c28e8Piotr Jasiukajtis	long double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
63225c28e8Piotr Jasiukajtis
63325c28e8Piotr Jasiukajtis/* INDENT OFF */
63425c28e8Piotr Jasiukajtis/*
63525c28e8Piotr Jasiukajtis * compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
63625c28e8Piotr Jasiukajtis *
63725c28e8Piotr Jasiukajtis *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
63825c28e8Piotr Jasiukajtis *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
63925c28e8Piotr Jasiukajtis *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
64025c28e8Piotr Jasiukajtis *       T2(j) = T2[2j,2j+1] = log(z[j]),
64125c28e8Piotr Jasiukajtis *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + ... + T3[6]s^15
64225c28e8Piotr Jasiukajtis *  Note
64325c28e8Piotr Jasiukajtis *  (1) the leading entries are truncated to 24 binary point.
64425c28e8Piotr Jasiukajtis *  (2) Remez error for T3(s) is bounded by 2**(-72.4)
64525c28e8Piotr Jasiukajtis *                                   2**(-24)
64625c28e8Piotr Jasiukajtis *                           _________V___________________
64725c28e8Piotr Jasiukajtis *               T1(n):     |_________|___________________|
64825c28e8Piotr Jasiukajtis *                             _______ ______________________
64925c28e8Piotr Jasiukajtis *               T2(j):       |_______|______________________|
65025c28e8Piotr Jasiukajtis *                                ____ _______________________
65125c28e8Piotr Jasiukajtis *               2s:             |____|_______________________|
65225c28e8Piotr Jasiukajtis *                                    __________________________
65325c28e8Piotr Jasiukajtis *          +    T3(s)-2s:           |__________________________|
65425c28e8Piotr Jasiukajtis *                       -------------------------------------------
65525c28e8Piotr Jasiukajtis *                          [leading] + [Trailing]
65625c28e8Piotr Jasiukajtis */
65725c28e8Piotr Jasiukajtis	/* INDENT ON */
65825c28e8Piotr Jasiukajtis	ix = H0_WORD(x);
65925c28e8Piotr Jasiukajtis	n2 = (ix >> 16) - 0x3fff;	/* exponent of x, range:3-10 */
66025c28e8Piotr Jasiukajtis	y = scalbnl(x, -n2);	/* y = scale x to [1,2] */
66125c28e8Piotr Jasiukajtis	n2 += n2;		/* 2n */
66225c28e8Piotr Jasiukajtis	j = (ix >> 10) & 0x3f;	/* j */
66325c28e8Piotr Jasiukajtis	z = 1.0078125L + (long double) j * 0.015625L;	/* z[j]=1+j/64+1/128 */
66425c28e8Piotr Jasiukajtis	j2 = j + j;
66525c28e8Piotr Jasiukajtis	t1 = y + z;
66625c28e8Piotr Jasiukajtis	t2 = y - z;
66725c28e8Piotr Jasiukajtis	r = one / t1;
66825c28e8Piotr Jasiukajtis	u = r * t2;		/* u = (y-z)/(y+z) */
66925c28e8Piotr Jasiukajtis	t1 = CHOPPED(t1);
67025c28e8Piotr Jasiukajtis	t4 = T2[j2 + 1] + T1[n2 + 1];
67125c28e8Piotr Jasiukajtis	z2 = u * u;
67225c28e8Piotr Jasiukajtis	k = H0_WORD(u) & 0x7fffffff;
67325c28e8Piotr Jasiukajtis	t3 = T2[j2] + T1[n2];
67425c28e8Piotr Jasiukajtis	for (t5 = T3[6], i = 5; i >= 0; i--)
67525c28e8Piotr Jasiukajtis		t5 = z2 * t5 + T3[i];
67625c28e8Piotr Jasiukajtis	if ((k >> 16) < 0x3fec) {	/* |u|<2**-19 */
67725c28e8Piotr Jasiukajtis		t2 = t4 + u * (two + z2 * t5);
67825c28e8Piotr Jasiukajtis	} else {
67925c28e8Piotr Jasiukajtis		t5 = t4 + (u * z2) * t5;
68025c28e8Piotr Jasiukajtis		u2 = u + u;
68125c28e8Piotr Jasiukajtis		v = (long double) ((int) (u2 * t24)) * p24;
68225c28e8Piotr Jasiukajtis		t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
68325c28e8Piotr Jasiukajtis		t3 += v;
68425c28e8Piotr Jasiukajtis	}
68525c28e8Piotr Jasiukajtis	ss_h = CHOPPED((t2 + t3));
68625c28e8Piotr Jasiukajtis	ss_l = t2 - (ss_h - t3);
68725c28e8Piotr Jasiukajtis/* INDENT OFF */
68825c28e8Piotr Jasiukajtis/*
68925c28e8Piotr Jasiukajtis * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
69025c28e8Piotr Jasiukajtis * where ss = log(x) - 1 in already in extra precision
69125c28e8Piotr Jasiukajtis */
69225c28e8Piotr Jasiukajtis	/* INDENT ON */
69325c28e8Piotr Jasiukajtis	z = one / x;
69425c28e8Piotr Jasiukajtis	r = x - half;
69525c28e8Piotr Jasiukajtis	r_h = CHOPPED((r));
69625c28e8Piotr Jasiukajtis	w_h = r_h * ss_h + hln2pim1_h;
69725c28e8Piotr Jasiukajtis	z2 = z * z;
69825c28e8Piotr Jasiukajtis	w = (r - r_h) * ss_h + r * ss_l;
69925c28e8Piotr Jasiukajtis	t1 = GP[19];
70025c28e8Piotr Jasiukajtis	for (i = 18; i > 0; i--)
70125c28e8Piotr Jasiukajtis		t1 = z2 * t1 + GP[i];
70225c28e8Piotr Jasiukajtis	w += hln2pim1_l;
70325c28e8Piotr Jasiukajtis	w_l = z * (GP[0] + z2 * t1) + w;
70425c28e8Piotr Jasiukajtis	k = (int) ((w_h + w_l) * invln2_32 + half);
70525c28e8Piotr Jasiukajtis
70625c28e8Piotr Jasiukajtis	/* compute the exponential of w_h+w_l */
70725c28e8Piotr Jasiukajtis
70825c28e8Piotr Jasiukajtis	j = k & 0x1f;
70925c28e8Piotr Jasiukajtis	*m = k >> 5;
71025c28e8Piotr Jasiukajtis	t3 = (long double) k;
71125c28e8Piotr Jasiukajtis
71225c28e8Piotr Jasiukajtis	/* perform w - k*ln2_32 (represent as w_h - w_l) */
71325c28e8Piotr Jasiukajtis	t1 = w_h - t3 * ln2_32hi;
71425c28e8Piotr Jasiukajtis	t2 = t3 * ln2_32lo;
71525c28e8Piotr Jasiukajtis	w = t2 - w_l;
71625c28e8Piotr Jasiukajtis	w_h = t1 - w;
71725c28e8Piotr Jasiukajtis	w_l = w - (t1 - w_h);
71825c28e8Piotr Jasiukajtis
71925c28e8Piotr Jasiukajtis	/* compute exp(w_h-w_l) */
72025c28e8Piotr Jasiukajtis	z = w_h - w_l;
72125c28e8Piotr Jasiukajtis	for (t1 = Et[10], i = 9; i >= 0; i--)
72225c28e8Piotr Jasiukajtis		t1 = z * t1 + Et[i];
72325c28e8Piotr Jasiukajtis	t3 = w_h - (w_l - (z * z) * t1);	/* t3 = expm1(z) */
72425c28e8Piotr Jasiukajtis	zz.l = S_trail[j] * (one + t3) + S[j] * t3;
72525c28e8Piotr Jasiukajtis	zz.h = S[j];
72625c28e8Piotr Jasiukajtis	return (zz);
72725c28e8Piotr Jasiukajtis}
72825c28e8Piotr Jasiukajtis
72925c28e8Piotr Jasiukajtis/* INDENT OFF */
73025c28e8Piotr Jasiukajtis/*
73125c28e8Piotr Jasiukajtis * kpsin(x)= sin(pi*x)/pi
73225c28e8Piotr Jasiukajtis *	           3        5        7        9        11                27
73325c28e8Piotr Jasiukajtis *	= x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x  + ... + ks[12]*x
73425c28e8Piotr Jasiukajtis */
73525c28e8Piotr Jasiukajtisstatic const long double ks[] = {
73625c28e8Piotr Jasiukajtis	-1.64493406684822643647241516664602518705158902870e+0000L,
73725c28e8Piotr Jasiukajtis	+8.11742425283353643637002772405874238094995726160e-0001L,
73825c28e8Piotr Jasiukajtis	-1.90751824122084213696472111835337366232282723933e-0001L,
73925c28e8Piotr Jasiukajtis	+2.61478478176548005046532613563241288115395517084e-0002L,
74025c28e8Piotr Jasiukajtis	-2.34608103545582363750893072647117829448016479971e-0003L,
74125c28e8Piotr Jasiukajtis	+1.48428793031071003684606647212534027556262040158e-0004L,
74225c28e8Piotr Jasiukajtis	-6.97587366165638046518462722252768122615952898698e-0006L,
74325c28e8Piotr Jasiukajtis	+2.53121740413702536928659271747187500934840057929e-0007L,
74425c28e8Piotr Jasiukajtis	-7.30471182221385990397683641695766121301933621956e-0009L,
74525c28e8Piotr Jasiukajtis	+1.71653847451163495739958249695549313987973589884e-0010L,
74625c28e8Piotr Jasiukajtis	-3.34813314714560776122245796929054813458341420565e-0012L,
74725c28e8Piotr Jasiukajtis	+5.50724992262622033449487808306969135431411753047e-0014L,
74825c28e8Piotr Jasiukajtis	-7.67678132753577998601234393215802221104236979928e-0016L,
74925c28e8Piotr Jasiukajtis};
75025c28e8Piotr Jasiukajtis/* INDENT ON */
75125c28e8Piotr Jasiukajtis
75225c28e8Piotr Jasiukajtis/*
75325c28e8Piotr Jasiukajtis * assume x is not tiny and positive
75425c28e8Piotr Jasiukajtis */
75525c28e8Piotr Jasiukajtisstatic struct LDouble
75625c28e8Piotr Jasiukajtiskpsin(long double x) {
75725c28e8Piotr Jasiukajtis	long double z, t1, t2;
75825c28e8Piotr Jasiukajtis	struct LDouble xx;
75925c28e8Piotr Jasiukajtis	int i;
76025c28e8Piotr Jasiukajtis
76125c28e8Piotr Jasiukajtis	z = x * x;
76225c28e8Piotr Jasiukajtis	xx.h = x;
76325c28e8Piotr Jasiukajtis	for (t2 = ks[12], i = 11; i > 0; i--)
76425c28e8Piotr Jasiukajtis		t2 = z * t2 + ks[i];
76525c28e8Piotr Jasiukajtis	t1 = z * x;
76625c28e8Piotr Jasiukajtis	t2 *= z * t1;
76725c28e8Piotr Jasiukajtis	xx.l = t1 * ks[0] + t2;
76825c28e8Piotr Jasiukajtis	return (xx);
76925c28e8Piotr Jasiukajtis}
77025c28e8Piotr Jasiukajtis
77125c28e8Piotr Jasiukajtis/* INDENT OFF */
77225c28e8Piotr Jasiukajtis/*
77325c28e8Piotr Jasiukajtis * kpcos(x)= cos(pi*x)/pi
77425c28e8Piotr Jasiukajtis *                     2        4        6        8        10        12
77525c28e8Piotr Jasiukajtis *	= 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x  +kc[5]*x
77625c28e8Piotr Jasiukajtis *
77725c28e8Piotr Jasiukajtis *                     2        4        6        8        10            22
77825c28e8Piotr Jasiukajtis *	= 1/pi - pi/2*x +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x  +...+kc[9]*x
77925c28e8Piotr Jasiukajtis *
78025c28e8Piotr Jasiukajtis * -pi/2*x*x = (npi_2_h + npi_2_l) * (x_f+x_l)*(x_f+x_l)
78125c28e8Piotr Jasiukajtis *	   =  npi_2_h*(x_f+x_l)*(x_f+x_l) + npi_2_l*x*x
78225c28e8Piotr Jasiukajtis *	   =  npi_2_h*x_f*x_f + npi_2_h*(x*x-x_f*x_f) + npi_2_l*x*x
78325c28e8Piotr Jasiukajtis *	   =  npi_2_h*x_f*x_f + npi_2_h*(x+x_f)*(x-x_f) + npi_2_l*x*x
78425c28e8Piotr Jasiukajtis * Here x_f = (long double) (float)x
78525c28e8Piotr Jasiukajtis * Note that pi/2(in hex) =
78625c28e8Piotr Jasiukajtis *  1.921FB54442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29
78725c28e8Piotr Jasiukajtis * npi_2_h = -pi/2 chopped to 25 bits = -1.921FB50000000000000000000000000 =
78825c28e8Piotr Jasiukajtis *  -1.570796310901641845703125000000000 and
78925c28e8Piotr Jasiukajtis * npi_2_l =
79025c28e8Piotr Jasiukajtis *  -0.0000004442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29 =
79125c28e8Piotr Jasiukajtis *  -.0000000158932547735281966916397514420985846996875529104874722961539 =
79225c28e8Piotr Jasiukajtis *  -1.5893254773528196691639751442098584699687552910487472296153e-8
79325c28e8Piotr Jasiukajtis * 1/pi(in hex) =
79425c28e8Piotr Jasiukajtis *  .517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
79525c28e8Piotr Jasiukajtis * will be splitted into:
79625c28e8Piotr Jasiukajtis *  one_pi_h = 1/pi chopped to 48 bits = .517CC1B727220000000000...  and
79725c28e8Piotr Jasiukajtis *  one_pi_l = .0000000000000A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
79825c28e8Piotr Jasiukajtis */
79925c28e8Piotr Jasiukajtis
80025c28e8Piotr Jasiukajtisstatic const long double
80125c28e8Piotr Jasiukajtis#if defined(__x86)
80225c28e8Piotr Jasiukajtisone_pi_h = 0.3183098861481994390487670898437500L,	/* 31 bits */
80325c28e8Piotr Jasiukajtisone_pi_l = 3.559123248900043690127872406891929148e-11L,
80425c28e8Piotr Jasiukajtis#else
80525c28e8Piotr Jasiukajtisone_pi_h = 0.31830988618379052468299050815403461456298828125L,
80625c28e8Piotr Jasiukajtisone_pi_l = 1.46854777018590994109505931010230912897495334688117e-16L,
80725c28e8Piotr Jasiukajtis#endif
80825c28e8Piotr Jasiukajtisnpi_2_h = -1.570796310901641845703125000000000L,
80925c28e8Piotr Jasiukajtisnpi_2_l = -1.5893254773528196691639751442098584699687552910e-8L;
81025c28e8Piotr Jasiukajtis
81125c28e8Piotr Jasiukajtisstatic const long double kc[] = {
81225c28e8Piotr Jasiukajtis	+1.29192819501249250731151312779548918765320728489e+0000L,
81325c28e8Piotr Jasiukajtis	-4.25027339979557573976029596929319207009444090366e-0001L,
81425c28e8Piotr Jasiukajtis	+7.49080661650990096109672954618317623888421628613e-0002L,
81525c28e8Piotr Jasiukajtis	-8.21458866111282287985539464173976555436050215120e-0003L,
81625c28e8Piotr Jasiukajtis	+6.14202578809529228503205255165761204750211603402e-0004L,
81725c28e8Piotr Jasiukajtis	-3.33073432691149607007217330302595267179545908740e-0005L,
81825c28e8Piotr Jasiukajtis	+1.36970959047832085796809745461530865597993680204e-0006L,
81925c28e8Piotr Jasiukajtis	-4.41780774262583514450246512727201806217271097336e-0008L,
82025c28e8Piotr Jasiukajtis	+1.14741409212381858820016567664488123478660705759e-0009L,
82125c28e8Piotr Jasiukajtis	-2.44261236114707374558437500654381006300502749632e-0011L,
82225c28e8Piotr Jasiukajtis};
82325c28e8Piotr Jasiukajtis/* INDENT ON */
82425c28e8Piotr Jasiukajtis
82525c28e8Piotr Jasiukajtis/*
82625c28e8Piotr Jasiukajtis * assume x is not tiny and positive
82725c28e8Piotr Jasiukajtis */
82825c28e8Piotr Jasiukajtisstatic struct LDouble
82925c28e8Piotr Jasiukajtiskpcos(long double x) {
83025c28e8Piotr Jasiukajtis	long double z, t1, t2, t3, t4, x4, x8;
83125c28e8Piotr Jasiukajtis	int i;
83225c28e8Piotr Jasiukajtis	struct LDouble xx;
83325c28e8Piotr Jasiukajtis
83425c28e8Piotr Jasiukajtis	z = x * x;
83525c28e8Piotr Jasiukajtis	xx.h = one_pi_h;
83625c28e8Piotr Jasiukajtis	t1 = (long double) ((float) x);
83725c28e8Piotr Jasiukajtis	x4 = z * z;
83825c28e8Piotr Jasiukajtis	t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);
83925c28e8Piotr Jasiukajtis	for (i = 8, t3 = kc[9]; i >= 0; i--)
84025c28e8Piotr Jasiukajtis		t3 = z * t3 + kc[i];
84125c28e8Piotr Jasiukajtis	t3 = one_pi_l + x4 * t3;
84225c28e8Piotr Jasiukajtis	t4 = t1 * t1 * npi_2_h;
84325c28e8Piotr Jasiukajtis	x8 = t2 + t3;
84425c28e8Piotr Jasiukajtis	xx.l = x8 + t4;
84525c28e8Piotr Jasiukajtis	return (xx);
84625c28e8Piotr Jasiukajtis}
84725c28e8Piotr Jasiukajtis
84825c28e8Piotr Jasiukajtis/* INDENT OFF */
84925c28e8Piotr Jasiukajtisstatic const long double
85025c28e8Piotr Jasiukajtis	/* 0.13486180573279076968979393577465291700642511139552429398233 */
85125c28e8Piotr Jasiukajtis#if defined(__x86)
85225c28e8Piotr Jasiukajtist0z1   =  0.1348618057327907696779385054997035808810L,
85325c28e8Piotr Jasiukajtist0z1_l =  1.1855430274949336125392717150257379614654e-20L,
85425c28e8Piotr Jasiukajtis#else
85525c28e8Piotr Jasiukajtist0z1   =  0.1348618057327907696897939357746529168654L,
85625c28e8Piotr Jasiukajtist0z1_l =  1.4102088588676879418739164486159514674310e-37L,
85725c28e8Piotr Jasiukajtis#endif
85825c28e8Piotr Jasiukajtis	/* 0.46163214496836234126265954232572132846819620400644635129599 */
85925c28e8Piotr Jasiukajtis#if defined(__x86)
86025c28e8Piotr Jasiukajtist0z2   =  0.4616321449683623412538115843295472018326L,
86125c28e8Piotr Jasiukajtist0z2_l =  8.84795799617412663558532305039261747030640e-21L,
86225c28e8Piotr Jasiukajtis#else
86325c28e8Piotr Jasiukajtist0z2   =  0.46163214496836234126265954232572132343318L,
86425c28e8Piotr Jasiukajtist0z2_l =  5.03501162329616380465302666480916271611101e-36L,
86525c28e8Piotr Jasiukajtis#endif
86625c28e8Piotr Jasiukajtis	/* 0.81977310110050060178786870492160699631174407846245179119586 */
86725c28e8Piotr Jasiukajtis#if defined(__x86)
86825c28e8Piotr Jasiukajtist0z3   =  0.81977310110050060178773362329351925836817L,
86925c28e8Piotr Jasiukajtist0z3_l =  1.350816280877379435658077052534574556256230e-22L
87025c28e8Piotr Jasiukajtis#else
87125c28e8Piotr Jasiukajtist0z3   =  0.8197731011005006017878687049216069516957449L,
87225c28e8Piotr Jasiukajtist0z3_l =  4.461599916947014419045492615933551648857380e-35L
87325c28e8Piotr Jasiukajtis#endif
87425c28e8Piotr Jasiukajtis;
87525c28e8Piotr Jasiukajtis/* INDENT ON */
87625c28e8Piotr Jasiukajtis
87725c28e8Piotr Jasiukajtis/*
87825c28e8Piotr Jasiukajtis * gamma(x+i) for 0 <= x < 1
87925c28e8Piotr Jasiukajtis */
88025c28e8Piotr Jasiukajtisstatic struct LDouble
88125c28e8Piotr Jasiukajtisgam_n(int i, long double x) {
88225c28e8Piotr Jasiukajtis	struct LDouble rr = {0.0L, 0.0L}, yy;
88325c28e8Piotr Jasiukajtis	long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
88425c28e8Piotr Jasiukajtis
88525c28e8Piotr Jasiukajtis	/* compute yy = gamma(x+1) */
88625c28e8Piotr Jasiukajtis	if (x > 0.2845L) {
88725c28e8Piotr Jasiukajtis		if (x > 0.6374L) {
88825c28e8Piotr Jasiukajtis			r1 = x - t0z3;
88925c28e8Piotr Jasiukajtis			r2 = CHOPPED((r1 - t0z3_l));
89025c28e8Piotr Jasiukajtis			t2 = r1 - r2;
89125c28e8Piotr Jasiukajtis			yy = GT3(r2, t2 - t0z3_l);
89225c28e8Piotr Jasiukajtis		} else {
89325c28e8Piotr Jasiukajtis			r1 = x - t0z2;
89425c28e8Piotr Jasiukajtis			r2 = CHOPPED((r1 - t0z2_l));
89525c28e8Piotr Jasiukajtis			t2 = r1 - r2;
89625c28e8Piotr Jasiukajtis			yy = GT2(r2, t2 - t0z2_l);
89725c28e8Piotr Jasiukajtis		}
89825c28e8Piotr Jasiukajtis	} else {
89925c28e8Piotr Jasiukajtis		r1 = x - t0z1;
90025c28e8Piotr Jasiukajtis		r2 = CHOPPED((r1 - t0z1_l));
90125c28e8Piotr Jasiukajtis		t2 = r1 - r2;
90225c28e8Piotr Jasiukajtis		yy = GT1(r2, t2 - t0z1_l);
90325c28e8Piotr Jasiukajtis	}
90425c28e8Piotr Jasiukajtis	/* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
90525c28e8Piotr Jasiukajtis	switch (i) {
90625c28e8Piotr Jasiukajtis	case 0:		/* yy/x */
90725c28e8Piotr Jasiukajtis		r1 = one / x;
90825c28e8Piotr Jasiukajtis		xh = CHOPPED((x));	/* x is not tiny */
90925c28e8Piotr Jasiukajtis		rr.h = CHOPPED(((yy.h + yy.l) * r1));
91025c28e8Piotr Jasiukajtis		rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) -
91125c28e8Piotr Jasiukajtis			r1 * yy.l);
91225c28e8Piotr Jasiukajtis		break;
91325c28e8Piotr Jasiukajtis	case 1:		/* yy */
91425c28e8Piotr Jasiukajtis		rr.h = yy.h;
91525c28e8Piotr Jasiukajtis		rr.l = yy.l;
91625c28e8Piotr Jasiukajtis		break;
91725c28e8Piotr Jasiukajtis	case 2:		/* (x+1)*yy */
91825c28e8Piotr Jasiukajtis		z = x + one;	/* may not be exact */
91925c28e8Piotr Jasiukajtis		zh = CHOPPED((z));
92025c28e8Piotr Jasiukajtis		rr.h = zh * yy.h;
92125c28e8Piotr Jasiukajtis		rr.l = z * yy.l + (x - (zh - one)) * yy.h;
92225c28e8Piotr Jasiukajtis		break;
92325c28e8Piotr Jasiukajtis	case 3:		/* (x+2)*(x+1)*yy */
92425c28e8Piotr Jasiukajtis		z1 = x + one;
92525c28e8Piotr Jasiukajtis		z2 = x + 2.0L;
92625c28e8Piotr Jasiukajtis		z = z1 * z2;
92725c28e8Piotr Jasiukajtis		xh = CHOPPED((z));
92825c28e8Piotr Jasiukajtis		zh = CHOPPED((z1));
92925c28e8Piotr Jasiukajtis		xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
93025c28e8Piotr Jasiukajtis
93125c28e8Piotr Jasiukajtis		rr.h = xh * yy.h;
93225c28e8Piotr Jasiukajtis		rr.l = z * yy.l + xl * yy.h;
93325c28e8Piotr Jasiukajtis		break;
93425c28e8Piotr Jasiukajtis
93525c28e8Piotr Jasiukajtis	case 4:		/* (x+1)*(x+3)*(x+2)*yy */
93625c28e8Piotr Jasiukajtis		z1 = x + 2.0L;
93725c28e8Piotr Jasiukajtis		z2 = (x + one) * (x + 3.0L);
93825c28e8Piotr Jasiukajtis		zh = CHOPPED(z1);
93925c28e8Piotr Jasiukajtis		zl = x - (zh - 2.0L);
94025c28e8Piotr Jasiukajtis		xh = CHOPPED(z2);
94125c28e8Piotr Jasiukajtis		xl = zl * (zh + z1) - (xh - (zh * zh - one));
94225c28e8Piotr Jasiukajtis
94325c28e8Piotr Jasiukajtis		/* wh+wl=(x+2)*yy */
94425c28e8Piotr Jasiukajtis		wh = CHOPPED((z1 * (yy.h + yy.l)));
94525c28e8Piotr Jasiukajtis		wl = (zl * yy.h + z1 * yy.l) - (wh - zh * yy.h);
94625c28e8Piotr Jasiukajtis
94725c28e8Piotr Jasiukajtis		rr.h = xh * wh;
948