SinCosTanDD.c   [plain text]

/*
 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
 *
 * @APPLE_LICENSE_HEADER_START@
 * 
 * The contents of this file constitute Original Code as defined in and
 * are subject to the Apple Public Source License Version 1.1 (the
 * "License").  You may not use this file except in compliance with the
 * License.  Please obtain a copy of the License at
 * http://www.apple.com/publicsource and read it before using this file.
 * 
 * This Original Code and all software distributed under the License are
 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT.  Please see the
 * License for the specific language governing rights and limitations
 * under the License.
 * 
 * @APPLE_LICENSE_HEADER_END@
 */
//
//	SinCosTanDD.c
//
//	Double-double Function Library
//	Copyright: © 1995-96 by Apple Computer, Inc., all rights reserved
//	
//	long double sinl( long double x );
//	long double cosl( long double x );
//	long double tanl( long double x );
//
//	Change History:			[7d]
//
//		06/20/96 PAF - Remove /* */ comments (for color editor)

#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"

static const float	k1_8To52	=  6755399441055744.0f;			// 0x1.8p52 
static const double	kPiBy2		=  1.570796326794896619231322;	// 0x1.921FB54442D18p0
static const double	kPiBy2Tail1	=  6.1232339957367660e-17;		// 0x1.1A62633145C07p-54
static const double	kPiBy2Tail2	= -1.4973849048591698e-33;		// 0x1.F1976B7ED8FBCp-110
static const double	kPiBy2Tail3	=  5.5622711043168264e-50;		// 0x1.4CF98E804177Dp-164
static const double	k2ByPi		=  0.636619772367581382;		// 0x1.45f306dc9c883p-1
static const float	k2To52		=  4503599627370496.0f;			// 0x1.0p+52
static const float	k2To54		=  18014398509481984.0f;			// 0x1.0p+54
static const double	k2To45m1	=  35184372088831.0;			// 0x1.fffffffffffffp44

static const uint32_t CosCoeff[] __attribute__ ((aligned(8))) = {
	0xbFE00000,		0x00000000,		0x39685556,		0x00000000,
	0x3FA55555,		0x55555555,		0x3C455555,		0x51E8C72E,
	0xbF56C16C,		0x16C16C17,		0x3BF02698,		0xA05BBBBD,
	0x3EFA01A0,		0x1A0169EB,		0x3B540FCD,		0xE3AE1BE0,
	0xbE927E4D,		0x648A174F,		0xbB361705,		0xD6F2C25B
};
	 
static const uint32_t SinCoeff[] __attribute__ ((aligned(8))) = {
     0xbFC55555,	0x55555555, 	0xbC655555,		0x55555552,
     0x3F811111,	0x11111111, 	0x3C011111,		0x0DBB7C30,
     0xbF2A01A0,	0x1A01A01A, 	0xbB510055,		0x358E4100,
     0x3EC71DE3,	0xA556A5C2, 	0x3B61A976,		0xF8644324,
     0xbE5AE642,	0x52F2B067, 	0xbAD356A3,		0xF6278F6C
};

static const uint32_t CosCoeff16[] __attribute__ ((aligned(8))) = {
     0xBFE00000,	0x00000000, 	0x395DA200,		0x00000000,
     0x3FA55555,	0x55555555, 	0x3C455555,		0x554FC86B,
     0xBF56C16C,	0x16C16C17, 	0x3BEF49FE,		0x788C8E18,
     0x3EFA01A0,	0x1A01A01A, 	0xBB6A022F,		0xCA316020,
     0xBE927E4F,	0xB7789A12, 	0xBB245398,		0xA829E8F5,
     0x3E21EED8,	0xEE008FB4, 	0x3AC49E38,		0xCA6157DE,
     0xBDA938BA,	0xA889516E, 	0x3A46CBC1,		0x4AF6E840
};

static const uint32_t SinCoeff16[] __attribute__ ((aligned(8))) = {
     0xBFC55555,	0x55555555, 	0xBC655555,		0x55555554,
     0x3F811111,	0x11111111, 	0x3C011111,		0x10F8EB9C,
     0xBF2A01A0,	0x1A01A01A, 	0xBB69FF8C,		0xFBC69E00,
     0x3EC71DE3,	0xA556C733, 	0x3B6EFC3D,		0xB602439A,
     0xBE5AE645,	0x67F53972, 	0xBAF653B7,		0xA917A91D,
     0x3DE61246,	0x1048C848, 	0x3A7CFCE5,		0x9E3CEA38,
     0xBD6AE6FF,	0xCFCF788E, 	0x39DA11EB,		0x64C11650
};

static const uint32_t TrigTable[] __attribute__ ((aligned(8))) = {
	0x3FB80000,	0x00054E6C, 								// angle 	
	0x3FB7F701,	0x032A9959, 	0xBBAD11F0,	0x5495F5BE,		// Sin Head, Sin Tail
	0x3FEFDC06,	0xBF7E5BB6, 	0xBBDCB529,	0x1917F935,		// Cos Head, Cos Tail
	0x3FBC0000,	0x000A5F0A, 
	0x3FBBF1B7,	0x8572884A, 	0x3BBAADFC,	0xF7E52ACA,
	0x3FEFCF0C,	0x800E7577, 	0xBBE4670B,	0xB43517A0,
	0x3FC00000,	0x000F1278, 
	0x3FBFEAAE,	0xEEA4D6F0, 	0xBB9D8153,	0x7343A887,
	0x3FEFC015,	0x527CE390, 	0x3BE16998,	0xA0C1E22D,
	0x3FC20000,	0x00052790, 	
	0x3FC1F0D3,	0xD7B4E938, 	0xBBB94CF6,	0x02B38414,
	0x3FEFAF22,	0x263C1D96, 	0xBBDE4F5F,	0x508C4D7E,
	0x3FC40000,	0x000143F6, 
	0x3FC3EB31,	0x2C5EA6CF, 	0x3BBF0AD1,	0xD6EF035D,
	0x3FEF9C34,	0x0A7CB78E, 	0xBBE87189,	0x89E0691A,
	0x3FC60000,	0x0007E6FF, 
	0x3FC5E44F,	0xCFA8F023, 	0xBBAF148F,	0x61A85B8E,
	0x3FEF874C,	0x2E1E9676, 	0x3BC38F0F,	0xAB6185B5,
	0x3FC80000,	0x000E77F8, 	
	0x3FC7DC10,	0x2FC929C2, 	0xBBC3C05A,	0x114DA314,
	0x3FEF706B,	0xDF9E2180, 	0xBBD8BEA1,	0x93F3F8C7,
	0x3FCA0000,	0x00102F36, 
	0x3FC9D252,	0xD0DE9D1A, 	0x3B8BF04E,	0xADF4B6FE,
	0x3FEF5794,	0x8CFE96A3, 	0xBBE58FFE,	0x7BB60718,
	0x3FCC0000,	0x0009CE78, 	
	0x3FCBC6F8,	0x4EE5F43E, 	0xBBA7E35C,	0xC27E82D6,
	0x3FEF3CC7,	0xC3B3493B, 	0xBBCC7CA2,	0xE7033121,
	0x3FCE0000,	0x00150434, 
	0x3FCDB9E1,	0x5FCA16EB, 	0x3BC107C0,	0xD41B55C5,
	0x3FEF2007,	0x30852C41, 	0x3BE70655,	0x75B2D6FB,
	0x3FD00000,	0x0061B811, 
	0x3FCFAAEE,	0xD5B07238, 	0xBBCB82BA,	0x01ECEAD1,
	0x3FEF0154,	0x9F71D818, 	0x3BE95043,	0xCC0B6426,
	0x3FD10000,	0x0020111B, 	
	0x3FD0CD00,	0xCF125569, 	0xBBC2D015,	0xDA2CFEAE,
	0x3FEEE0B1,	0xFBBCBB9D, 	0xBBE9AA1F,	0xF00E2374,
	0x3FD20000,	0x000C1AD3, 	
	0x3FD1C37D,	0x64D25988, 	0x3BC4B236,	0xE14BCAA2,
	0x3FEEBE21,	0x4F75419B, 	0x3BE52928,	0x2D338F16,
	0x3FD30000,	0x000E7F89, 		
	0x3FD2B8DD,	0xC44C91CC, 	0xBBC5707F,	0x432D5073,
	0x3FEE99A4,	0xC3A5AEA3, 	0x3BD23919,	0x2EA440E3,
	0x3FD40000,	0x0017A801, 	
	0x3FD3AD12,	0x9780568B, 	0xBBC59FD1,	0x16DF09C5,
	0x3FEE733E,	0xA0159A53, 	0xBBB12BF5,	0xDB2DCDCC,
	0x3FD50000,	0x0039F7ED, 	
	0x3FD4A00C,	0x9B461D51, 	0xBBCDFE0E,	0x1445AAA9,
	0x3FEE4AF1,	0x4B20ED63, 	0x3BD374F6,	0x9F86063A,
	0x3FD60000,	0x0000A1AA, 	
	0x3FD591BC,	0x9FA38DCC, 	0xBBCD07DB,	0x418C1954,
	0x3FEE20BF,	0x49ACBB83, 	0x3B864FA4,	0xEB634A4E,
	0x3FD70000,	0x00011BC1, 
	0x3FD68213,	0x8A39E197, 	0xBBBA6968,	0x3C12FB6D,
	0x3FEDF4AB,	0x3EBD5578, 	0xBBD6766E,	0x746C48E7,
	0x3FD80000,	0x0008ED50, 	
	0x3FD77102,	0x557E9094, 	0x3BB488F1,	0x4A8DF658,
	0x3FEDC6B7,	0xEB97B68C, 	0x3BDCA325,	0xE2462D16,
	0x3FD90000,	0x00861E6C, 
	0x3FD85E7A,	0x12FE6D57, 	0x3BA65076,	0x5FCB7841,
	0x3FED96E8,	0x2F58212A, 	0xBBEC65EC,	0x2CEB0E13,
	0x3FDA0000,	0x00021040, 
	0x3FD94A6B,	0xE9F72C06, 	0x3BA04A4B,	0xE0FC13D6,
	0x3FED653F,	0x073DD7E0, 	0x3BEC1A4A,	0xF1A7037E,
	0x3FDB0000,	0x005EF712, 	
	0x3FDA34C9,	0x1D1BB055, 	0x3BC42083,	0xD43A5519,
	0x3FED31BF,	0x8D7A0AAD, 	0xBB9D5F11,	0x5DAECA9B,
	0x3FDC0000,	0x000C075C, 
	0x3FDB1D83,	0x053CFB6A, 	0xBBDA72E9,	0x9A1C2FBE,
	0x3FECFC6C,	0xFA50214C, 	0xBBCB5164,	0xBC69A105,
	0x3FDD0000,	0x0056E33E, 
	0x3FDC048B,	0x17FF5F2B, 	0x3BD8AF58,	0xF8525AD2,
	0x3FECC54A,	0xA29F9263, 	0x3BD22EA3,	0xDBB55D0D,
	0x3FDE0000,	0x005CCC54,
	0x3FDCE9D2,	0xE4276EF1, 	0xBBAAEDF9,	0xDD5DC899,
	0x3FEC8C5B,	0xF8B9244D, 	0xBBC2D9C6,	0xCE4671F4,
	0x3FDF0000,	0x00161118, 
	0x3FDDCD4C,	0x154623DC, 	0xBBD97F5B,	0x8DF7FEBF,
	0x3FEC51A4,	0x8B85F41A, 	0x3BC354E8,	0x86ABC9CE,
	0x3FE00000,	0x0013AF20, 
	0x3FDEAEE8,	0x746D926F, 	0x3BD2B7DF,	0xA4EEC7DD,
	0x3FEC1528,	0x06520D6D, 	0x3BE9EE87,	0x2601172E,
	0x3FE08000,	0x0046DF51, 
	0x3FDF8E99,	0xE7E60D68, 	0x3BAF1859,	0x1DD2BAB4,
	0x3FEBD6EA,	0x30DFA2E0, 	0x3BE95328,	0x8BAB7ACB,
	0x3FE10000,	0x0002BED1, 
	0x3FE03629,	0x39C8F748, 	0x3BCFAF90,	0x9F47E8F6,
	0x3FEB96EE,	0xEF572000, 	0x3BD89C91,	0x79738697,
	0x3FE18000,	0x00404634, 
	0x3FE0A402,	0x1ED4F66C, 	0x3BD4E9FD,	0x2527C4F7,
	0x3FEB553A,	0x40EAA3C8, 	0xBBE9F185,	0x9F0B73D8,
	0x3FE20000,	0x0005ABAA, 
	0x3FE110D0,	0xC4BB683B, 	0xBBB78FD4,	0x4D302311,
	0x3FEB11D0,	0x415F9E99, 	0x3BDAE453,	0xDBCB3704,
	0x3FE28000,	0x002D02D5, 	
	0x3FE17C8E,	0x5F549FEE, 	0x3B6368B3,	0xE72D341D,
	0x3FEACCB5,	0x26DE0531, 	0xBB8CE182,	0xC27CD7B8,
	0x3FE30000,	0x0002B1A8, 
	0x3FE1E734,	0x323892EB, 	0x3BC9A07C,	0xF3B89184,
	0x3FEA85ED,	0x43725E55, 	0x3BE2C483,	0xBE536F66,
	0x3FE38000,	0x000730F9, 
	0x3FE250BB,	0x937E7157, 	0xBBC3FA7C,	0xADF21ED4,
	0x3FEA3D7D,	0x034EA01E, 	0xBBD5786E,	0x86ABB722,
	0x3FE40000,	0x00076EEF, 
	0x3FE2B91D,	0xEA8E4953, 	0x3BE1CB0F,	0x13421234,
	0x3FE9F368,	0xED8CD61E, 	0xBBC0A358,	0x32653793,
	0x3FE48000,	0x00119ECA, 
	0x3FE32054,	0xB156DCB6, 	0xBBBF1C5F,	0x95EB1A26,
	0x3FE9A7B5,	0xA35FDCFF, 	0x3BE6C143,	0x88B5E209,
	0x3FE50000,	0x000045B6, 
	0x3FE38659,	0x74565F66, 	0x3BD5144E,	0x65AA2BBD,
	0x3FE95A67,	0xE00C8774, 	0x3BC1A72A,	0xE0D050C3,
	0x3FE58000,	0x00A4E7B8, 	
	0x3FE3EB25,	0xD3EDE59C, 	0x3BE1F420,	0x5560DD7C,
	0x3FE90B84,	0x77E73599, 	0x3B9250E3,	0x5B17FE76,
	0x3FE60000,	0x000049F4, 	
	0x3FE44EB3,	0x81CF7192, 	0xBBD5B364,	0x7D4D895F,
	0x3FE8BB10,	0x5A5D9A12, 	0xBBD2BB98,	0xF4F2C91C,
	0x3FE68000,	0x0008C197, 	
	0x3FE4B0FC,	0x46B16552, 	0x3BC863B5,	0xB735A540,
	0x3FE86910,	0x8D71FD5A, 	0x3BE5C4E3,	0x402BA3EE,
	0x3FE70000,	0x00340EE5,
	0x3FE511F9,	0xFDA26352,		0x3BCDC369, 0x16E2A8DB,
	0x3FE8158A,	0x316F2658, 	0x3BD69C38, 0x4E2E4085,
	0x3FE78000,	0x005756FD, 
	0x3FE571A6,	0x96AE2DA7, 	0x3BE24CA0,	0xA2E35E02,
	0x3FE7C082,	0x7ECF5E07, 	0xBBCBD891,	0x2A83004E,
	0x3FE80000,	0x004601DC, 
	0x3FE5CFFC,	0x16F2C847, 	0x3BCA6614,	0x963A9B5C,
	0x3FE769FE,	0xC62568E3, 	0x3BE746CD,	0xBABDB543,
	0x3FE88000,	0x002EF948, 
	0x3FE62CF4,	0x99438A17, 	0x3BA3ED99,	0x63704218,
	0x3FE71204,	0x6F86E919, 	0xBBE0671D,	0x4098C629,
	0x3FE90000,	0x0022F620, 
	0x3FE6888A,	0x4E2C1E13, 	0xBBC0D332,	0x374F8A22,
	0x3FE6B898,	0xFA865CFA, 	0xBBDE0573,	0x1D87548D,
	0x3FE98000,	0x00883EB0, 
	0x3FE6E2B7,	0x7C9FF82A, 	0x3BE241BC,	0xA153A6CC,
	0x3FE65DC1,	0xFD8A1C59, 	0xBBDBAD2A,	0x038DFB75
};

static const uint32_t TanCoeff16[] __attribute__ ((aligned(8))) = {
     0x3FD55555,	0x55555555,   0x3C755555,	0x5555554C,
     0x3FC11111,	0x11111111,   0x3C411111,	0x11327FD4,
     0x3FABA1BA,	0x1BA1BA1C,   0xBC47917B,	0x8730ECE8,
     0x3F9664F4,	0x882C10FA,   0xBC08946C,	0x04337BC0,
     0x3F8226E3,	0x55E6C239,   0xBC16F3F8,	0x5EB82A9C,
     0x3F6D6D3D,	0x0E159BD4,   0xBC0BF02E,	0x91A39520,
     0x3F57DA36,	0x44E2BF4A,   0x3BF4CC81,	0xFDF36B54,
     0x3F435582,	0xB4E9E8D9,   0xBBEAAB03,	0xE323ED60,
     0x3F2F5712,	0xDCD1E82C,   0x3BA67F16,	0xB4DC7710,
     0x3F19C9D2,	0x46D3C681,   0x3BA2B1F2,	0xBC709028
};

static const uint32_t CotCoeff16[] __attribute__ ((aligned(8))) = {
     0x3FD55555,	0x55555555,   0x3C755555,	0x55555550,
     0x3F96C16C,	0x16C16C17,   0xBC2F49F4,	0x9F46AA0E,
     0x3F61566A,	0xBC011567,   0xBC050FFB,	0x0703E386,
     0x3F2BBD77,	0x9334EF0B,   0xBBC541A3,	0x5FEBE630,
     0x3EF66A8F,	0x2BF70EDF,   0xBB960EEE,	0x58853E10,
     0x3EC22805,	0xD642B849,   0x3B49A324,	0x92C33F58,
     0x3E8D6DB2,	0xD58EAF82,   0x3B103D9B,	0xCE4CCBB8,
     0x3E57DA1C,	0x0FF6D02F,   0x3AFBDB54,	0xE9F15B6C,
     0x3E2396A6,	0x6CC1B8E2,   0x3AC04FB7,	0x0AEC86D2
};

static const uint32_t TanCoeff[] __attribute__ ((aligned(8))) = {
     0x3ff00000,	0x00000000,   0x00000000,	0x00000000,
     0x3FD55555,	0x55555555,   0x3C755555,	0x55555549,
     0x3FC11111,	0x11111111,   0x3C411111,	0x17771B40,
     0x3FABA1BA,	0x1BA1BA1C,   0xBC47A471,	0x53F3999C,
     0x3F9664F4,	0x882C119E,   0x3C2A931E,	0x61CFDED8,
     0x3F8226E3,	0x5546FD80,   0xBC22700A,	0xDC5E76C2,
     0x3F6D6DCB,	0xC6BCC1FD,   0x3C0EB3DD,	0x1EBC7C6C
};

static const uint32_t TanTable[] __attribute__ ((aligned(8))) = {
	0x3FB80000,	0x00095B76, 
	0x3FB81210,	0x420B0DDD,   0xBBADBC24,	0xD03EE61B,
	0x40254552,	0xEE62D182,   0xBC23A778,	0xD4D4C82D,
	0x3FBC0000,	0x00F6D8DC,   
	0x3FBC1CB8,	0x85A84FC9,   0xBBA0734F,	0xD6E646D7,
	0x40223676,	0x163ABE75,   0x3C19C932,	0xFC451F09,
	0x3FC00000,	0x002177FC,   
	0x3FC01577,	0xAF3710E9,   0xBBBB1D68,	0x25554EB9,
	0x401FD549,	0xF004A93B,   0x3BF1B140,	0x750F85AD,
	0x3FC20000,	0x0000C2F1,   
	0x3FC21E9E,	0x0175E475,   0x3BBD91B3,	0xA31924A6,
	0x401C41B6,	0xE169DC8D,   0x3C057A49,	0x2ABEE9CF,
	0x3FC40000,	0x004D703E,   
	0x3FC42A13,	0xDFCB15A0,   0xBB94ED3B,	0x3620DA54,
	0x4019642D,	0xFDBA3F8F,   0x3BBB47D8,	0x80CDA3C6,
	0x3FC60000,	0x000CAC24,   
	0x3FC6381F,	0x200F2AED,   0x3BBF3967,	0x0381AC58,
	0x40170B09,	0x205DE76C,   0xBBFEFF3A,	0x14B71AB8,
	0x3FC80000,	0x00473296,   
	0x3FC84906,	0xF15CE818,   0x3BB03E11,	0x36E0F2D3,
	0x4015152E,	0xCDA71D37,   0xBBE77628,	0xDB83D025,
	0x3FCA0000,	0x00B2D7A8,   
	0x3FCA5D14,	0x0081E4D7,   0xBBC0188C,	0x9C6AB2DE,
	0x40136BB4,	0xBA046E82,   0xBC0F80BA,	0xA34DE4A8,
	0x3FCC0000,	0x0066CDE2,   
	0x3FCC7490,	0xA23DC3B4,   0xBBCB362F,	0xA69FE4DA,
	0x4011FE3C,	0xA58F0345,   0xBBFCB31C,	0x2F2152B4,
	0x3FCE0000,	0x00252978,   
	0x3FCE8FC9,	0x01177F3A,   0xBBA5E8BA,	0x10809B66,
	0x4010C0C5,	0xABFB3524,   0xBBFA678B,	0xCA50DEEC,
	0x3FD00000,	0x00638D21,   
	0x3FD05785,	0xA4A65715,   0x3BC29365,	0x7FDA2145,
	0x400F549E,	0x3203FCA2,   0x3BE2971C,	0x11C2A8A3,
	0x3FD10000,	0x002C35A8,   
	0x3FD16953,	0xEACF2DA7,   0x3BC40A91,	0xBA0FE60F,
	0x400D67EC,	0xF327EFBB,   0x3C05D1EC,	0x4B5FCFCB,
	0x3FD20000,	0x001C3AB2,   
	0x3FD27D78,	0xB42A0CE9,   0x3BCCED21, 	0x016F79E8,
	0x400BB0C1,	0xF1319BD2,   0xBC051F68,	0x2CBAE56E,
	0x3FD30000,	0x000F9C12,   
	0x3FD3941E,	0xADA8C533,   0xBBBB6B83,	0x08722CFA,
	0x400A26A8,	0xA4D1FD5D,   0xBBF11692,	0xFACCF74A,
	0x3FD40000,	0x000B6423,   
	0x3FD4AD71,	0xED5E62C2,   0xBBAD5590,	0x3A1538DA,
	0x4008C2DD,	0x5F133FBD,   0xBBDDCA83,	0x1B096361,
	0x3FD50000,	0x0017ECB5,   
	0x3FD5C9A0,	0x105DB3D0,   0x3BCE9A3C,	0xD13B3D17,
	0x40077FE6,	0x3A7531D9,   0x3BE036AE,	0x48AA7EC4,
	0x3FD60000,	0x00785229,   
	0x3FD6E8D8,	0x5AEC50EC,   0xBBC83B79,	0x6A9E4244,
	0x40065948,	0x26E6F5E4,   0xBBFE36C3,	0x52F1803B,
	0x3FD70000,	0x001C8263,   
	0x3FD80B4B,	0xD8D44659,   0xBBB499F2,	0xDC83C2D8,
	0x40054B4F,	0x852854B9,   0x3BFB7522,	0xCB8BEDBE,
	0x3FD80000,	0x0092A4EC,   
	0x3FD9312D,	0x86455645,   0x3B93DC76,	0x9413BABA,
	0x400452E6,	0x95D5E9C0,   0xBBDBF037,	0x35429FA9,
	0x3FD90000,	0x00353578,   
	0x3FDA5AB2,	0x70323F27,   0xBBD5E11A,	0x4F5FE60A,
	0x40036D75,	0xEB84FD4A,   0xBBE3A384,	0xFD52C107,
	0x3FDA0000,	0x002177A0,   
	0x3FDB8811,	0xE4F7B2DE,   0x3BC5B70A,	0x7D15AD3F,
	0x400298CC,	0x1A74BDDF,   0x3BD41900,	0x83E8A76D,
	0x3FDB0000,	0x0071F1F6,   
	0x3FDCB985,	0x9CFB26A0,   0x3BC8BA6E,	0xF094701B,
	0x4001D30A,	0xD7B1694B,   0xBBFECEAB,	0x555C2BC9,
	0x3FDC0000,	0x0071BD74,   
	0x3FDDEF49,	0xEB35D726,   0xBBCCC6A5,	0xDDB5912C,
	0x40011A98,	0x2075779C,   0xBBF95852,	0x2DBF4C5B,
	0x3FDD0000,	0x00A18446,   
	0x3FDF299D,	0xF3CB9E1D,   0xBBDAB47B,	0x6F4E4DBE,
	0x40006E12,	0x72BBF45D,   0x3BE9A288,	0x2C8CD3D9,
	0x3FDE0000,	0x0009F5E0,   
	0x3FE03461,	0xF090AA3D,   0xBBB6411F,	0x2E55E33D,
	0x3FFF988E,	0xC8064412,   0xBBE51DE0,	0x7EC2ADBC,
	0x3FDF0000,	0x00263EAE,   
	0x3FE0D680,	0x92D62104,   0x3BDD6620,	0x8142D974,
	0x3FFE6858,	0x0C88DF44,   0xBBDD8335,	0xD4838950,
	0x3FE00000,	0x002FD5FE,   
	0x3FE17B4F,	0x5C31640E,   0xBBCF7FD9,	0x4387A68E,
	0x3FFD49AD,	0x7DDFB168,   0x3BF71CBF,	0x13808D41,
	0x3FE08000,	0x00760E88,   
	0x3FE222F4,	0xAFFFC590,   0x3BDD3CE8,	0x22236C45,
	0x3FFC3AF4,	0x70DCC003,   0xBBFBC120,	0xA2F987AC,
	0x3FE10000,	0x000F8806,   
	0x3FE2CD98,	0xFEB5905C,   0x3BE169C8,	0x3A579A2C,
	0x3FFB3AC2,	0x7605E88C,   0xBBC11C38,	0xB96450D3,
	0x3FE18000,	0x002B5E5D,   
	0x3FE37B66,	0xF43D000E,   0x3BD2830E,	0xBAC9646C,
	0x3FFA47D6,	0x6F7CCB9D,   0xBBEB41B3,	0x171B37E1,
	0x3FE20000,	0x000A3C10,   
	0x3FE42C8B,	0xA0F7A0E3,   0x3BE38DF3,	0x66A12B85,
	0x3FF96112,	0xDACBFCBC,   0xBBEFEF64,	0xEF1EA7D1,
	0x3FE28000,	0x00595422,   
	0x3FE4E136,	0xB0CFA783,   0xBBBA398C,	0xEE061521,
	0x3FF88578,	0xFAA0A3DE,   0x3BDA1845,	0x4AE391EC,
	0x3FE30000,	0x0026966D,   
	0x3FE5999A,	0x9E477C73,   0xBBCE4AC5,	0x97E84259,
	0x3FF7B424,	0xCEF83049,   0xBBF35B98,	0x9FE75083,
	0x3FE38000,	0x0015B3FD,   
	0x3FE655EC,	0xF397B5BE,   0xBBE51A33,	0x112B4404,
	0x3FF6EC49,	0xA1E2DBE5,   0x3BB86EE4,	0x93E45292,
	0x3FE40000,	0x00146EE3,   
	0x3FE71666,	0x89F330B5,   0xBBE35378,	0x5A4867EE,
	0x3FF62D2F,	0x215A3E9B,   0x3BEA0C9E,	0xF65CC2CF,
	0x3FE48000,	0x000F229A,   
	0x3FE7DB43,	0xD3A2EEDF,   0xBBE0A454,	0xC31E28C0,
	0x3FF5762E,	0xE1176AE6,   0xBBE1C3BD,	0x493469EA,
	0x3FE50000,	0x002E680F,   
	0x3FE8A4C5,	0x2CF1484A,   0xBB9432B2,	0x063946DC,
	0x3FF4C6B2,	0x384B4991,   0x3BD8E9D0,	0xF949DFE6,
	0x3FE58000,	0x00672096,   
	0x3FE9732F,	0x3459A187,   0x3BDDADD7,	0x84F3E87C,
	0x3FF41E30,	0x6C72DBD4,   0xBBEC67AE,	0xAE1EBB1F,
	0x3FE60000,	0x00B6C23F,   
	0x3FEA46CB,	0x2D189DAF,   0x3BB739F4,	0x75A8271C,
	0x3FF37C2D,	0x1B3BE6EB,   0xBBC85203,	0x8CECB4C0,
	0x3FE68000,	0x0028A725,   
	0x3FEB1FE7,	0x6A3CF1FF,   0xBBEA7D29,	0x7EE6B9D1,
	0x3FF2E036,	0xDC04F43A,   0x3BED69F5,	0xAF37688F,
	0x3FE70000,	0x004CA480,   
	0x3FEBFED7,	0xCD2DB8E2,   0xBBCD123B,	0x0537EE2D,
	0x3FF249E6,	0x09D6492E,   0xBBD38032,	0x80A8D415,
	0x3FE78000,	0x000FFDA8,   
	0x3FECE3F6,	0x42FE6158,   0x3BD569F8,	0x793799F7,
	0x3FF1B8DB,	0xBED5846E,   0x3BD9797E,	0x39EAD48F,
	0x3FE80000,	0x001E1B0A,   
	0x3FEDCFA3,	0x61492AAF,   0xBBDA3DEC,	0xB61278D5,
	0x3FF12CC0,	0xE58AA7DC,   0x3BB3088E,	0x1168567F,
	0x3FE88000,	0x000B208A,   
	0x3FEEC247,	0x07D4CEEF,   0xBBE3FF76,	0x1F807DFF,
	0x3FF0A545,	0x6F3C80B0,   0xBBCF3891,	0xE44E8437,
	0x3FE90000,	0x000657F8,   
	0x3FEFBC51,	0x1E0226B6,   0xBBEC6FF5,	0x5F02EDF6,
	0x3FF0221F,	0x9DADF743,   0xBBD42807,	0xFD02BF0B,
	0x3FE98000,	0x00655AE2,   
	0x3FF05F1D,	0x31752EE7,   0xBBD8BE76,	0xB0E002E9,
	0x3FEF4616,	0xC9018897,   0xBBDB4E55,	0xDEAF57E4
};
 
#warning  These probably ought to static or typedefs below. They are being exported.

struct TrigTblEntry { 
   double arg;
   double sinHead;
   double sinTail;
   double cosHead;
   double cosTail;
} TrigTblEntry;

struct CosCoeffEntry { 
   double Head;				
   double Tail;
} CosCoeffEntry;

struct SinCoeffEntry { 
   double Head;				
   double Tail;
} SinCoeffEntry;

struct CosCoeff16Entry { 
   double Head;				
   double Tail;
} CosCoeff16Entry;

struct SinCoeff16Entry { 
   double Head;				
   double Tail;
} SinCoeff16Entry;

struct TanTblEntry { 
   double arg;
   double TanHead;
   double TanTail;
   double CotHead;
   double CotTail;
} TanTblEntry;

struct TanCoeffEntry { 
   double Head;				
   double Tail;
} TanCoeffEntry;

struct CotCoeff16Entry { 
   double Head;				
   double Tail;
} CotCoeff16Entry;

struct TanCoeff16Entry { 
   double Head;				
   double Tail;
} TanCoeff16Entry;	  

static void argReduce( double xHead, double xTail, double *dLeft, int *M, double *redHead, double *redTail )   
{
	double amid, alow, carry, high, q, q1;
	double b, blow, b1, c, clow, c1, c2, c3, c4;
	double d, dlow, d1, d2, d3, d4, d5, e1, e2, e3;
	Double DeN;
	amid = 0.0;
	alow = 0.0;
	
	while (fabs(xHead) > k2To45m1)    		// reduce until argument
	{											// isn't larger than k2To45m1
		q = xHead*k2ByPi; 		     
	   	if (fabs(q) < k2To54) {	     		// maybe q(mod 4) != 0
			if (fabs(q) < k2To52) {     		// maybe q isn't even an integer
				if (q > 0)
					q = __FADD( q, k2To52 ) - k2To52;	// remove fraction if any exists
				else
					q = __FSUB( q, k2To52 ) + k2To52;
	      	}			     					// at this point, q is an integer
			q1 = q*0.25; 		     			// and determine q(mod 4)
			if (q > 0)
				q1 = __FADD( q1, k2To52 ) - k2To52;		// giving the quadrant in which the angle resides.
			else
				q1 = __FSUB( q1, k2To52 ) + k2To52;		// right adjust quadrant
			DeN.f = __FNMSUB( 4.0, q1, q ) + k1_8To52; // q - 4.0*q1 + k1_8To52;		     
	      
			*M = *M + DeN.i[1];
	   	}

		//****************************************************************************
		//   Tackling the subtraction of q*(pi/2) from the octuple precision word:   *
		//		      (xHead, xTail, amid, alow)			   						 *
		//****************************************************************************
	   
		high = __FNMSUB( q, kPiBy2, xHead );	// xHead - q*kPiBy2;		     	// exact (Spectacular use of MAF)
		b = q*kPiBy2Tail1;
		blow = __FNMSUB( q, kPiBy2Tail1, b );	// b - q*kPiBy2Tail1;		     	// Yet another use of MAF
		c = q*kPiBy2Tail2;
		clow = __FNMSUB( q, kPiBy2Tail2, c );	// c - q*kPiBy2Tail2;
		d = q*kPiBy2Tail3;
		dlow = __FNMSUB( q, kPiBy2Tail3, d );	// d - q*kPiBy2Tail3;
		b1 = __FADD( xTail, high );
		if (fabs(xTail) > fabs(high)) 		// sum and propagate carries to the right
			c1 = xTail - b1 + high;
		else
			c1 = high - b1 + xTail;
		xHead = __FSUB( b1, b );	
		if (fabs(b1) > fabs(b))
			c2 = b1 - xHead - b + c1;
		else
			c2 = b1 - (xHead+b) + c1;
		c3 = __FADD( amid,  blow);
		if (fabs(amid) > fabs(blow))
			d1 = amid - c3 + blow;
		else
			d1 = blow - c3 + amid;
		c4 = __FSUB( c2, c );
		if (fabs(c2) > fabs(c))
			d2 = c2 - c4 - c + d1;
		else
			d2 = c2 - __FADD( c4, c ) + d1;
		xTail = __FADD( c3, c4 );
		if (fabs(c3) > fabs(c4))
			d3 = c3 - xTail + c4 + d2;
		else
			d3 = c4 - xTail + c3 + d2;
		d4 = __FADD( alow, clow );
		if (fabs(alow) > fabs(clow))
			e1 = alow - d4 + clow;
		else
			e1 = clow - d4 + alow;
		d5 = __FSUB( d3, d );
		if (fabs(d3) > fabs(d))
			e2 = d3 - d5 - d + e1;
		else
			e2 = d3 - (d5+d) + e1;
		amid = __FADD( d4, d5 );
		if (fabs(d4) > fabs(d5))
			e3 = d4 - amid + d5 + e2;
		else
			e3 = d5 - amid + d4 + e2;
		alow = e3 + dlow;
	}
	
	high  = __FADD( xHead, xTail );		 // At long last,
	carry = xHead - high + xTail;		 // A final distillation
	
	xTail  = __FADD( carry, __FADD( amid, alow ) );
	xHead  = high;
	*dLeft = __FSUB( carry, xTail ) + __FADD( amid, alow );
		
	*redHead = xHead;
	*redTail = xTail;
}
	
static void SinCosCore(int M, double xHead, double xTail, double dleft, double *yHead, double *yTail )
{
	double intquo;
	double temp, temps, tempc; 
	double b, c, ca, d; 
	double redarg, redarg1, res, reslow, resmid, resbot, restop, resin, resextra;
	double ressup, restiny, resinf; 
	double p1, p2, p3, p4, p5, p6, p6a,  p7, p8, p9, p10;
	double p11, p12, p13, p14, p15, p16, p17, p18, p19, p20;
	double p21, p22, p23, p23a, p23b, p24, p25, p26, p27, p27a, p28, p29, p30, p31;
	double arg, arg2, arg2a, argsq, argsq2, argmid, argnew, argres;  
	double sin, sin2, sint, sinl, sinargtiny;
	double sintab, sintab2, arg2b;
	double cos, cos2, cosl, cost;
	double prods, prodsl, prodc, prodcl;
	double blow, carry, cerr;
	int i, index;
	Double DeN;
	
	struct TrigTblEntry *TrigPtr = (struct TrigTblEntry *)TrigTable - 6;
	struct CosCoeffEntry *CosPtr = (struct CosCoeffEntry *) CosCoeff;
	struct SinCoeffEntry *SinPtr = (struct SinCoeffEntry *) SinCoeff;
	struct CosCoeff16Entry *Cos16Ptr = (struct CosCoeff16Entry *) CosCoeff16;
	struct SinCoeff16Entry *Sin16Ptr = (struct SinCoeff16Entry *) SinCoeff16;

	DeN.f = __FMADD( xHead, k2ByPi, k1_8To52 );		// xHead*k2ByPi + k1_8To52;				// # of multiples of pi/2 to remove
       	
	intquo = DeN.f - k1_8To52;						// Integer
    redarg = __FNMSUB( intquo, kPiBy2, xHead );		// xHead - intquo*kPiBy2;				// Exact result
	   
	M = M + DeN.i[1];								// M=quadrant in which to comp. sin
       
	b = kPiBy2Tail1*intquo;							// continue argument reduction
	blow = __FMSUB( kPiBy2Tail1, intquo, b );		// kPiBy2Tail1*intquo - b;
	redarg1 = __FSUB( xTail, b );															// second word of reduced argument
	DeN.f = __FMADD( fabs(redarg), 64.0, k1_8To52 );// fabs(redarg)*64.0 + k1_8To52;		// index into acctbl
	index = DeN.i[1];

	ca = kPiBy2Tail2*intquo;

	if (fabs(b) > fabs(xTail))
		carry = xTail - (b+redarg1);    			// tail-b may be inexact
	else
		carry = (xTail-redarg1) - b;
	   
	if (! (index < 6)) {
		if (redarg > 0.0) { 
			redarg = redarg - TrigPtr[index].arg ;	// arg-intquo ln2-tbl(index) ->
			sintab = TrigPtr[index].sinHead ;		// fetch sine of table value
			sintab2 = TrigPtr[index].sinTail ;
		} 
		else {										// much trailing zeroes in redarg
			redarg = redarg + TrigPtr[index].arg;
			sintab =- TrigPtr[index].sinHead ;		// fetch sine of table value
			sintab2 =- TrigPtr[index].sinTail ;
		}
	}
	   
	arg = __FADD( redarg, redarg1 );				// get signif. bits from the right
	arg2a = (redarg - arg + redarg1);
	arg2b = __FSUB( arg2a, blow );             		// fill in with bits from 2nd term
	c = __FSUB( carry, ca );                   		// 3rd word of reduced argument

	argmid = __FADD( arg2b , c );
	argnew = __FADD( arg, argmid );
	argres = arg - argnew + argmid;
	arg = argnew;
	cerr = (carry-c) - ca;

	if (fabs(arg2b) > fabs(c))
		d = arg2b - argmid + c + dleft;
	else
		d = c - argmid + arg2b + dleft;

	// d = d + cerr + (ca-kPiBy2Tail2*intquo) - kPiBy2Tail3*intquo; // 4th word
	d = __FNMSUB( kPiBy2Tail3, intquo, d + cerr + __FNMSUB( kPiBy2Tail2, intquo, ca ) );
	arg2 = __FADD( argres, d );								// reduced arg is (arg, arg2)
	temp = 2.0*arg*arg2;
	argsq = __FMADD( arg, arg, temp );				// arg*arg + temp;
	argsq2 = __FMADD( arg2, arg2, __FMSUB( arg, arg, argsq ) + temp ); // arg*arg - argsq + temp + arg2*arg2;
       
	   
	// Now we have the reduced arg = (arg, arg2)  
	   
	if (index < 6) {								// straight power series
          
		if ((M & 1) == 0) { 						// sin series evaluation
			sin = Sin16Ptr[6].Head ;				// need extra terms when
			sinl = Sin16Ptr[6].Tail ;				// used without table lookup
			sinargtiny = argres - arg2 + d;

			for(i = 5; i > -1; i--) {
				temps = __FMADD( sin, argsq2, sinl*argsq ); // sin*argsq2 + sinl*argsq;
				prods = __FMADD( sin, argsq, temps );		// sin*argsq + temps;
				prodsl = __FMSUB( sin, argsq, prods ) + temps; // sin*argsq - prods + temps;
				sin = __FADD( Sin16Ptr[i].Head, prods );
				sin2 = Sin16Ptr[i].Tail + prodsl;
				sinl = Sin16Ptr[i].Head - sin+prods + sin2;
			}
		} 
		else {                     					// cos series evaluation
			cos = Cos16Ptr[6].Head;
			cosl = Cos16Ptr[6].Tail ;
             
			for(i = 5; i > -1; i--) {
				temps = __FMADD( cos, argsq2, cosl*argsq ); // cos*argsq2 + cosl*argsq;
				prods = __FMADD( cos, argsq, temps );		// cos*argsq + temps;
				prodsl = __FMSUB( cos, argsq, prods) + temps; // cos*argsq - prods + temps;
				cos = __FADD( Cos16Ptr[i].Head, prods );
				cos2 = Cos16Ptr[i].Tail + prodsl;
				cosl = Cos16Ptr[i].Head - cos + prods + cos2;
			}
		}
	} 
	else {											// shorter series in case

		cos = CosPtr[4].Head; 						// where exact table is used.
		sin = SinPtr[4].Head;
		cosl = CosPtr[4].Tail;
		sinl = SinPtr[4].Tail;
		for (i = 3; i > -1; i--) {					// add short series
			tempc = __FMADD( cos, argsq2, cosl*argsq ); // cos*argsq2 + cosl*argsq;
			prodc = __FMADD( cos, argsq, tempc );	   // cos*argsq + tempc;
			prodcl = __FMSUB( cos, argsq, prodc) + tempc; // cos*argsq - prodc + tempc;
			cos = __FADD( CosPtr[i].Head, prodc );
			cos2 = CosPtr[i].Tail + prodcl;
			cosl = CosPtr[i].Head - cos + prodc + cos2;
			temps = __FMADD( sin, argsq2, sinl*argsq ); // sin*argsq2 + sinl*argsq;
			prods = __FMADD( sin, argsq, temps );		// sin*argsq + temps;
			prodsl = __FMSUB( sin, argsq, prods ) + temps; // sin*argsq - prods + temps;
			sin = __FADD( SinPtr[i].Head, prods );
			sin2 = SinPtr[i].Tail + prodsl;
			sinl = SinPtr[i].Head - sin + prods + sin2;
		}
	}
	   
	tempc = __FMADD( cos, argsq2, cosl*argsq );		// cos*argsq2 + cosl*argsq;
	cost = __FMADD( cos, argsq, tempc );			// cos*argsq + tempc;
	cosl = __FMSUB( cos, argsq, cost ) + tempc;		// cos*argsq - cost + tempc;
	cos = cost;                          			// cos(arg) - 1.0 complete
	temps = __FMADD( sin, argsq2, sinl*argsq );		// sin*argsq2 + sinl*argsq;
	sint = __FMADD( sin, argsq, temps );			// sin*argsq + temps;
	sinl = __FMSUB( sin, argsq, sint ) + temps;		// sin*argsq - sint + temps;
	temps = __FMADD( sint, arg2, sinl * arg );		// sint*arg2 + sinl*arg;
	sin = __FMADD( sint, arg, temps );				// sint*arg + temps;
	sinl = __FMSUB( sint, arg, sin ) + temps;		// sint*arg - sin + temps;         			// sin(arg) - arg complete
	   
	if (index < 6) {
		if ((M & 1) == 0) {
			res = __FADD( arg, sin );               // sin of a small angle
			reslow = __FADD( arg2, sinl );          // careful distallation
			resmid = __FADD( arg - res, sin );
			resbot = arg2 - reslow + sinl;
			restop = __FADD( res, (reslow+resmid) );
			resin = __FADD( res - restop, (reslow+resmid) );
			resextra = resmid - (reslow+resmid) + reslow;
			ressup = __FADD( restop, (resin+resbot) );
			restiny = __FSUB( resin, (resin+resbot) + resbot );
			resinf = restop - ressup + (resin+resbot) + (resextra+restiny+sinargtiny);
			*yHead = ressup;
			*yTail = resinf;
		} 
		else {
			res = __FADD( 1.0, cos );               			// cos of a small angle
			reslow = cosl;
			resmid = __FADD( 1.0 - res, cos );
			restop = __FADD( res, (reslow+resmid) );
			resin = __FADD( res - restop, (reslow+resmid) );
			resextra = resmid - (reslow+resmid) + reslow;
			ressup = __FADD( restop, resin );
			resinf = restop - ressup + resin + resextra;
			*yHead = ressup;
			*yTail = resinf;
		}
	}    // end of if (index < 6)
	else {
		if ((M & 1) == 0) {        					// even quadrant-- eval. sine(x)
			p1 = cos*sintab;																// o(-8)  exact
			p2 = __FMSUB( cos, sintab, p1 );		// cos*sintab - p1;          			// o(-62) exact
			p3 = __FMADD( cosl, sintab, cos*sintab2 ); //cosl*sintab + (cos*sintab2); 		// o(-67)
			p4 = arg* TrigPtr[index].cosHead ;												// o(-4)  exact
			p5 = __FMSUB(arg, TrigPtr[index].cosHead, p4 ); // arg* TrigPtr[index].cosHead -p4;	// o(-58) exact
			p6 = arg2* TrigPtr[index].cosHead ;												// o(-58)
			p6a = __FMSUB( arg2, TrigPtr[index].cosHead, p6 ); // arg2* TrigPtr[index].cosHead -p6;
			p7 = sin* TrigPtr[index].cosHead ;												// o(-4)  exact
			p8 = __FMSUB( sin, TrigPtr[index].cosHead, p7 ); // sin* TrigPtr[index].cosHead -p7; // o(-59) exact
			p9 = __FMADD( sinl, TrigPtr[index].cosHead, (sin* TrigPtr[index].cosTail ) ); //sinl* TrigPtr[index].cosHead +(sin* TrigPtr[index].cosTail );
			p10 = __FMADD( arg, TrigPtr[index].cosTail, p9 + p8); // arg* TrigPtr[index].cosTail +p9+p8;  // o(-59)
			p11 = __FADD( p2, p5 );															// o(-58) exact
			p12 = p5 - p11 + p2 + p6a;														// o(-112) exact
			p13 = __FADD( p6, p11 );														// o(-57) exact
			p15 = __FADD( sintab, p4 );														// o(0)   exact
			p16 = sintab - p15 + p4;														// o(-54) exact
			p17 = sintab2 + p10 + p12 + p3;													// o(-58)
			p18 = __FADD( p15, p1 );														// o(0) exact
			p20 = __FADD( p18, p7 );														// o(0) exact
			if (fabs(p6) > fabs(p11))
				p14 = p6 - p13 + p11;														// o(-111)
			else
				p14 = p11 - p13 + p6;
			p19 = p15 - p18 + p1 + p16;														// o(-53) exact
			p21 = p18 - p20 + p7;															// o(-54) exact
			p22 = p13 + p17;																// o(-56)
			p23 = __FADD( p19, p21 );														// o(-52) exact
			p25 = __FADD( p20, p23 );														// o(0) exact
			p26 = p20 - p25 + p23;															// o(-54) exact
			if (fabs(p19) > fabs(p21))
				p24 = p19 - p23 + p21;														// o(-108)
			else
				p24 = p21 - p23 + p19;
			p27 = __FADD( (p14+p24), p22 + p26 );											// o(-54)
			p27a = p26 - p27 + ((p14+p24)+p22);
			p28 = __FADD( p25, p27 );
	
			*yHead = p28;
			*yTail = p25 - p28 + p27 + p27a;
		}
		// end of if ((M & 1) == 0) 
		else {         								// evaluation of cos series.
			p1 = TrigPtr[index].cosHead *cos;													// o(-16) exact
			p2 = __FMSUB( TrigPtr[index].cosHead, cos, p1 ); // TrigPtr[index].cosHead *cos-p1;	// o(-71) exact
			p5 = sintab*arg;																	// o(-8) exact
			p6 = __FMSUB( sintab, arg, p5 );		// sintab*arg - p5 ;						// o(-62) exact
			p7 = sintab*arg2;																	// o(-62)
			p8 = sintab*sin;																	// o(-24) exact
			p9 = __FMSUB( sintab, sin, p8 );		// sintab*sin - p8 ;						// o(-78) exact
			p12 = __FMSUB( TrigPtr[index].cosTail, cos, (sintab2*arg) ); // TrigPtr[index].cosTail *cos - (sintab2*arg); // o(-80) 
			p13 = __FNMSUB( sintab2, sin, p12 );	// p12 - (sintab2*sin);						// o(-72)
			p14 = p13 + p2;																		// o(-70)
			p15 = __FMADD( TrigPtr[index].cosHead, cosl, p14 ); // p14 + TrigPtr[index].cosHead *cosl;	// o(-69)
			p16 = __FNMSUB( sintab, sinl, p15 - p9 );// p15 - p9 - sintab*sinl;					// o(-67)
			p17 = __FADD( TrigPtr[index].cosHead, p1 );											// o(-1) exact
			p18 = TrigPtr[index].cosHead - p17 + p1;											// o(-55) exact-16 bits
			p19 = __FSUB( p17, p5 );															// o(-1) exact
			p20 = p17 - p19 - p5;																// o(-55) exact-8 bits
			p21 = p18 + p20 ;																	// o(-54) exact-17 bits
			p22 = __FSUB( p19, p8 );															// o(-1) exact
			p23 = p19 - p22 - p8;																// o(-53) exact-25 bits
			p24 = TrigPtr[index].cosTail + p16;													// o(-65)
			p23a = __FADD( p23, p21 );															// o(-54)
			p25 = __FSUB( p23a, p6 );															// o(-53) exact
			p26 = p23a - p25 - p6;																// o(-106)
			p27 = __FSUB( p25, p7 );															// o(-53) exact
			if (fabs(p23) > fabs(p21))
				p23b = p23 - p23a + p21;														// o(-106)
			else
				p23b = p21 - p23a + p23;
			p28 = p25 - p27 - p7 + p26 + p23b;													// o(-105)
			p29 = __FADD( p27, p24 );															// o(-53) exact
			p30 = p27 - p29 + p24 + p28;														// o(-104)
			p31 = __FADD( p22, p29 );															// o(-1) exact
	
			*yHead = p31;
			*yTail = p22 - p31 + p29 + p30;
		}   //end of else if ((M & 1) == 0) 

	}	//end of else if (index < 6)
       
	   
	if((M & 2) != 0) {								// change sign for below x axis.
		*yHead = - *yHead;
		*yTail = - *yTail;
	}
}

long double sinl( long double x )
{
	uint32_t hibits;
	double xHead, xTail, dleft, fenv, yHead, yTail;
	DblDbl ddx;
	int M;

        FEGETENVD(fenv);
        FESETENVD(0.0);
	
	ddx.f = x;
	
	xHead = ddx.d[0];
	xTail = ddx.d[1];

	hibits = ddx.i[0] & 0x7FFFFFFFu;				// highest 32 bits as uint32_t

	// NaNs and infinities
	
	if (hibits >= 0x7ff00000u) {					// special cases: NaN and infinity
		if (xHead != xHead) {						// x is a NaN
			FESETENVD(fenv);
			return x;
		}
		else {
			ddx.d[0] = xHead - xHead;				// inf - inf gives NaN
			ddx.d[1] = ddx.d[0];
			FESETENVD(fenv);
			return (ddx.f);
		}
	}

	// x = zero hence xHead is zero
	if ((hibits | ddx.i[1]) == 0) {
                FESETENVD(fenv);
		return x;
	}
		
	// x is not infinity, NaN, or zero
    
	//for tiny x, sin x = x
	if (hibits < 0x29900000u) { 					// 2^(-358);  358 = (1022 + 52 )/3
                FESETENVD(fenv);
		return x;
	}
		
	M = 0;	// M = 0 for Sine, M = 1 for Cosine
	dleft = 0.0;
	if (hibits >= 0x42d00000u) {					// |x| > = 2 ^ 46
		argReduce(xHead, xTail, &dleft, &M, &xHead, &xTail);
	}   
	
	//  Call core algorithm
	SinCosCore( M, xHead, xTail, dleft, &yHead, &yTail );
	
	ddx.d[0] = yHead + yTail; 						// Cannonize result
	ddx.d[1] = (yHead - ddx.d[0]) + yTail; 
	
        FESETENVD(fenv);
	return (ddx.f);
}

long double cosl( long double x )
{
	uint32_t hibits;
	double xHead, xTail, dleft, fenv, yHead, yTail;
	DblDbl ddx;
	int M;

        FEGETENVD(fenv);
        FESETENVD(0.0);
	
	ddx.f = x;
	
	xHead = ddx.d[0];
	xTail = ddx.d[1];

	hibits = ddx.i[0] & 0x7FFFFFFFu;				// highest 32 bits as uint32_t

	// NaNs and infinities
	
	if (hibits >= 0x7ff00000u) {					// special cases: NaN and infinity
		if (xHead != xHead) {						// x is a NaN
			FESETENVD(fenv);
			return x;
		}
		else {
			ddx.d[0] = xHead - xHead;				// inf - inf gives NaN
			ddx.d[1] = ddx.d[0];
			FESETENVD(fenv);
			return (ddx.f);
		}
	}

	// x = zero hence xHead is zero
	if ((hibits | ddx.i[1]) == 0) {
                FESETENVD(fenv);
		return (long double) 1.0;
	}
		
	// x is not infinity, NaN, or zero
    
	//for tiny x, cos x = 1.0
	if (hibits < 0x1e600000u) { 					// 2^(-537);  537 = (1022 + 52 )/2
                FESETENVD(fenv);
		return 1.0;
	}
		
	M = 1;	// M = 0 for Sine, M = 1 for Cosine
	dleft = 0.0;
	if (hibits >= 0x42d00000u) {					// |x| > = 2 ^ 46
		argReduce(xHead, xTail, &dleft, &M, &xHead, &xTail);
	}   
	
	//  Call core algorithm
	SinCosCore( M, xHead, xTail, dleft, &yHead, &yTail );
	
	ddx.d[0] = __FADD( yHead, yTail ); 						// Cannonize result
	ddx.d[1] = (yHead - ddx.d[0]) + yTail; 
	
        FESETENVD(fenv);
	return (ddx.f);
}

long double tanl( long double x )
{
    double temp, temps; 
	double b, c, c1, c2,  ca, cerr, d; 
	double frac, frac2, frac3, t1, t2;
	double rec, res, reslow, resmid, resbot, restop, resin, resextra;
	double ressup, restiny, resinf; 
	double rarg, rarg0, rarg2, rarg3, rarg2a, redarg, redarg1, redargold;
	double p1, p2, p3, p4, p5, p6,  p7, p8, p9, p10;
	double p11, p12, p13, p14, p15, p16, p17, p18, p19, p20;
	double p21, p22, p23, p24, p25, p26, p27, p28, p29, p30;
	double p31, p32, p33, p34, p35, p36, p37, p38, p39, p40;
	double p41, p42, p43, p44, p45, p46;
	double arg, arg2, arg2a, arg2b, argsq, argsq2, argres, argmid, argnew;  
	double cot, cot2, cotl, cott,  tan, tan2, tant, tanl;
	double cotargtiny, tanargtiny; 
	double prods, prodsl, q;
	double blow, carry, dleft;
	int i, index, M;
	
	uint32_t hibits;
    Double DeN;
	double fenv;

	double xHead, xTail;
	DblDbl ddx;

	struct TanTblEntry *TanTablePtr = (struct TanTblEntry *)TanTable - 6;  
	struct TanCoeffEntry *TanPtr = (struct TanCoeffEntry *) TanCoeff - 1;  
	struct TanCoeff16Entry *Tan16Ptr = (struct TanCoeff16Entry *) TanCoeff16 - 1;  
	struct CotCoeff16Entry *Cot16Ptr = (struct CotCoeff16Entry *) CotCoeff16 - 1;  

        FEGETENVD(fenv);
        FESETENVD(0.0);

	ddx.f = x;
	
	xHead = ddx.d[0];
	xTail = ddx.d[1];
	
	hibits = ddx.i[0] & 0x7FFFFFFFu; 				// highest 32 bits as uint32_t
	
	if (hibits >= 0x7ff00000u){					// special cases: NaN and infinity
		if (xHead != xHead) {						// x is a NaN
			FESETENVD(fenv);
			return x;
		}
		else {
			ddx.d[0] = xHead - xHead;				// inf - inf gives NaN & raises invalid
			ddx.d[1] = ddx.d[0];
			FESETENVD(fenv);
			return (ddx.f);
		}
	}
		
	// x is not infinity nor a NaN
    
	// x = zero
	if ((hibits | ddx.i[1]) == 0) {
                FESETENVD(fenv);
		return x;
	}


	// finite x 

	//for tiny x, tan x = x
	if (hibits < 0x29900000u) {					// 2^(-358);  358 = (1022 + 52 )/3
                FESETENVD(fenv);
		return x;
	}

	// finite result but not tiny

	M=0;											// determine quadrant
	dleft=0.0;

    if (hibits >= 0x42b00000u){					// |x| > = 2 ^ 44
		argReduce (xHead, xTail, &dleft, &M, &xHead, &xTail);
	}   
	
	// now finite result < 2 ^ 44	

	DeN.f = __FMADD( xHead, k2ByPi, k1_8To52 );		// xHead*k2ByPi + k1_8To52;				// # of multiples of pi/2 to remove
	q = DeN.f - k1_8To52;							// to remove integer
    redarg = __FNMSUB( q, kPiBy2, xHead );			// xHead - intquo*kPiBy2;				// Exact result
	
	M = M + DeN.i[1];								// M=quadrant in which to compute sin
	b = kPiBy2Tail1*q;								// more argument reduction
	blow = __FMSUB( kPiBy2Tail1, q, b );			// kPiBy2Tail1*q - b;
	redarg1 = __FSUB( xTail, b );															// second word of reduced argument
	DeN.f = __FMADD( fabs(redarg), 64.0, k1_8To52 );// fabs(redarg)*64.0 + k1_8To52;		// index into acctbl
	ca = kPiBy2Tail2*q;
	
	if (fabs(b) > fabs(xTail))
		carry = xTail - (b+redarg1);				// tail-b may be inexact
	else
		carry = (xTail-redarg1) - b;
	
	redargold = redarg;
	index = DeN.i[1];
	  
	if (!(index < 6)) {
		if (redargold > 0.0)
			redarg = redarg - TanTablePtr[index].arg;	// arg-q ln2-tbl(index) ->
														// mucho trailing zeroes in redarg
		else
			redarg = redarg + TanTablePtr[index].arg;
	}
	
	arg = __FADD( redarg, redarg1 );				// get signif. bits from the right
	arg2a = (redarg - arg + redarg1);
	arg2b = __FSUB( arg2a, blow );             		// fill in with bits from 2nd term
	c = __FSUB( carry, ca );                   		// 3rd word of reduced argument

	argmid = __FADD( arg2b , c );
	argnew = __FADD( arg, argmid );
	argres = arg - argnew + argmid;
	arg = argnew;
	cerr = (carry-c) - ca;
	
	if (fabs(arg2b) > fabs(c))
		d = arg2b - argmid + c + dleft;
	else
		d = c - argmid + arg2b + dleft;
	
	// d = d + cerr + (ca-kPiBy2Tail2*q) - kPiBy2Tail3*q; // 4th word
	d = __FNMSUB( kPiBy2Tail3, q, d + cerr + __FNMSUB( kPiBy2Tail2, q, ca ) );
	arg2 = __FADD( argres, d );								// reduced arg is (arg, arg2)
	temp = 2.0*arg*arg2;
	argsq = __FMADD( arg, arg, temp );				// arg*arg + temp;
	argsq2 = __FMADD( arg2, arg2, __FMSUB( arg, arg, argsq ) + temp ); // arg*arg - argsq + temp + arg2*arg2;

	if (index < 6) {								// straight power series
		if ((M & 1) == 0) { 						// 2n multiple of pi/2 -> use TAN
			tan = Tan16Ptr[10].Head;				// need extra terms when
			tanl = Tan16Ptr[10].Tail;				// used without table lookup
			tanargtiny = argres - arg2 + d;
			
			for(i = 9; i > 0; i--) {
				temps = __FMADD( tan, argsq2, tanl*argsq );		// tan*argsq2 + tanl*argsq;
				prods = __FMADD( tan, argsq, temps );			// tan*argsq + temps;
				prodsl = __FMSUB( tan, argsq, prods ) + temps;	// tan*argsq - prods + temps;
				tan = __FADD( Tan16Ptr[i].Head, prods );
				tan2 = Tan16Ptr[i].Tail + prodsl;
				tanl = Tan16Ptr[i].Head - tan + prods + tan2;
			}
			temps = __FMADD( tan, argsq2, tanl*argsq );			// tan*argsq2 + tanl*argsq;
			tant = __FMADD( tan, argsq, temps );				// tan*argsq + temps;
			tanl = __FMSUB( tan, argsq, tant ) + temps;			// tan*argsq - tant + temps;
			temps = __FMADD( tant, arg2, tanl*arg );			// tant*arg2 + tanl*arg;
			tan = __FMADD( tant, arg, temps );					// tant*arg + temps;
			tanl = __FMSUB( tant, arg, tan ) + temps;			// tant*arg - tan + temps;			// tan(arg)-arg complete
			res = __FADD( arg, tan );																// tan of a small angle
			reslow = __FADD( arg2, tanl );															// distill VERY CAREFULLY
			resmid = arg - res + tan;
			resbot = arg2 - reslow + tanl;
			restop = __FADD( res, (reslow+resmid) );
			resin = res - restop + (reslow+resmid);
			resextra = resmid - (reslow+resmid) + reslow;
			ressup = __FADD( restop, (resin+resbot) );
			restiny = resin - (resin+resbot) + resbot;
			resinf = restop - ressup + (resin+resbot) + ( resextra+restiny+tanargtiny);
			
			//getting a cannonical result
			
			ddx.d[0] = __FADD( ressup, resinf ); 			
			ddx.d[1] = (ressup - ddx.d[0]) + resinf; 
			
			FESETENVD(fenv);
			return (ddx.f);
		}
		else {										// near an 2k+1 multiple of pi/2
													// tan (x)=-cot(x+(2k+1) pi/2)
			cot = Cot16Ptr[9].Head;					// need extra terms when
			cotl = Cot16Ptr[9].Tail;				// used without table lookup
			cotargtiny = argres - arg2 + d;
			rarg0 = 1.0/arg;						// We use the cotangent
			
			for(i = 8; i > 0; i--){ 				// series, which is one term shorter
				temps = __FMADD( cot, argsq2, cotl*argsq );		// cot*argsq2 + cotl*argsq;
				prods = __FMADD( cot, argsq, temps );			// cot*argsq + temps;
				prodsl = __FMSUB( cot, argsq, prods ) + temps;	// cot*argsq - prods + temps;
				cot = __FADD( Cot16Ptr[i].Head, prods );
				cot2 = Cot16Ptr[i].Tail + prodsl;
				cotl = Cot16Ptr[i].Head - cot + prods + cot2;
			}
			cott = cot;
			rarg = rarg0;
			temps = __FMADD( cott, arg2, cotl*arg );		// cott*arg2 + cotl*arg;
			cot = __FMADD( cott, arg, temps );				// cott*arg + temps;
			cotl = __FMSUB( cott, arg, cot ) + temps;		// cott*arg - cot + temps;			// 1/arg-cot(arg) completed
			p1 = __FNMSUB( arg, rarg, 1.0 );				// 1.0 - arg*rarg;					// next complete hi precision
			p2 = arg2*rarg;																		// reciprocal of argument, to subtract
			p3 = __FSUB( p1, p2 );																		// from series.  The reciprocal will
			rarg2a = p3*rarg;																	// be the dominant term.
			rarg2 = rarg2a + __FNMSUB( arg, rarg2a, p3 )*rarg;	// rarg2a + (p3-arg*rarg2a)*rarg
			
			if (fabs(p1) > fabs(p2))
				p4 = p1 - p3 - p2;
			else
				p4 = p1 - (p3+p2);
			
			p5 = p4 + __FNMSUB( arg2, rarg, p2 );			// p4 + (p2-arg2*rarg)
			p6 = __FNMSUB( cotargtiny, rarg, p5 );			// p5 - cotargtiny*rarg;
			p7 = __FNMSUB( rarg2, arg, p3 );				// p3 - rarg2*arg;
			rarg3 = __FNMSUB( arg2, rarg2, p6+p7 )*rarg;	// ((p6+p7) - arg2*rarg2)*rarg;   	// low order recip. bits
			
			//**********************************************************************
			//     BEGIN the subtraction of 1/arg from the series                  *
			//              |cot| much greater than  |rarg|                        *
			//**********************************************************************
			
			p8 = __FSUB( cot, rarg );
			p9 = cot - (rarg+p8);
			p10 = __FSUB( cotl, rarg2 );
			p11 = cotl - (rarg2+p10);
			p12 = p11 - rarg3;
			p13 = __FADD( p9, p10 );
			
			if (fabs(p9) > fabs(p10))
				p14 = p9 - p13 + p10;
			else
				p14 = p10 - p13 + p9;
			
			p15 = p12 + p14;
			p16 = __FADD( p8, p13 );         		// HO term of result
			p17 = p8 - p16 + p13;					// almost LO term
			p18 = p17 + p15;						// assimilate any other LO bits
			
			// getting a cannonical result
			
			ddx.d[0] = __FADD( p16, p18 ); 			
			ddx.d[1] = (p16 - ddx.d[0]) + p18; 
			
			FESETENVD(fenv);
			return (ddx.f);
		}
	}
	else {											// Serious cases
													// Accurate tables will be used.
		tant = TanPtr[7].Head;						// Need fewer terms when
		tant = __FMADD( tant, argsq, TanPtr[6].Head );		// TanPtr[6].Head + tant*argsq;
		tan = __FMADD( tant, argsq, TanPtr[5].Head );		// TanPtr[5].Head + tant*argsq;
		tanl = __FMADD( tant, argsq, TanPtr[5].Head - tan );// TanPtr[5].Head - tan + tant*argsq;
		tanargtiny = argres - arg2 + d;
		
		for(i = 4; i > 0; i--) {
			temps = __FMADD( tan, argsq2, tanl*argsq );		// tan*argsq2 + tanl*argsq;
			prods = __FMADD( tan, argsq, temps );			// tan*argsq + temps;
			prodsl = __FMSUB( tan, argsq, prods ) + temps;	// tan*argsq - prods + temps;
			tan = __FADD( TanPtr[i].Head, prods );
			tan2 = TanPtr[i].Tail + prodsl;
			tanl = TanPtr[i].Head -tan + prods + tan2;
		}
		tant = tan;
		temps = __FMADD( tant, arg2, tanl*arg );			// tant*arg2 + tanl*arg;
		tan = __FMADD( tant, arg, temps );					// tant*arg + temps;
		tanl = __FMSUB( tant, arg, tan ) + temps;			// tant*arg - tan + temps;				// tan(arg) complete
		if (redargold > 0) {																		// read out table values
			t1 = TanTablePtr[index].TanHead;														// (t1,t2) are tan(x)
			t2 = TanTablePtr[index].TanTail;
			c1 = TanTablePtr[index].CotHead;														// (c1,c2) are cot(x)
			c2 = TanTablePtr[index].CotTail;
		}
		else {
			t1 = -TanTablePtr[index].TanHead;
			t2 = -TanTablePtr[index].TanTail;
			c1 = -TanTablePtr[index].CotHead;
			c2 = -TanTablePtr[index].CotTail;
		}
		
		p1 = __FADD( t1, c1 );						// Tan(x)+Cot(x) needed
		p2 = __FADD( t2, c2 );						// for all quadrants
		p3 = c1 - p1 + t1;
		p4 = c2 - p2 + t2;
		p5 = __FADD( p3, p2 );
		p6 = p3 - p5 + p2 + p4;						// p1+p5+p6=Tan(x)+Cot(x) exactly
													// where x = table lookup argument
													// Next, this sum needs to be multiplied
													// by the tangent of the reduced argument
		p7 = p1*tan;
		p8 = __FMSUB( p1, tan, p7 );				// p1*tan - p7;
		p9 = p5*tan;
		p10 = __FMSUB( p5, tan, p9 );				// p5*tan - p9;
		p11 = p1*tanl;
		p12 = __FMSUB( p1, tanl, p11 );				// p1*tanl - p11;
		p13 = __FMADD( p5, tanl, p10 + p12 );		// p10 + p12 + p5*tanl;
		p14 = __FMADD( p6, tan, p13 );				// p13 + p6*tan;
		p15 = __FADD( p8, p9 );
		p17 = __FADD( p15, p11 );
		
		if (fabs(p8) > fabs(p9))
			p16 = p8 - p15 + p9;
		else
			p16 = p9 - p15 + p8;
		
		if (fabs(p11) > fabs(p15))
			p18 = p11 - p17 + p15 + (p14+p16);
		else
			p18 = p15 - p17 + p11 + (p14+p16);
		
		p19 = __FADD( p7, p17 );
		p20 = p7 - p19 + p17;
		p21 = __FADD( p20, p18 );
		p22 = p20 - p21 + p18;
		
		//***********************************************************************
		//         At this point, different paths are taken for arguments       *
		//         "neighborhood" of even and odd multiples of pi/2.            *
		//***********************************************************************
		
		if ((M & 1) == 0) {
			
			//***********************************************************************
			//         Let the argument A be written as A=x+y, where                *
			//         x=TanTablePtr[index].arg.  Then if T=tan(x) and C=cot(x)     *
			//         as read from tantbl, and t=tan(y),  (computed above as       *
			//         (tan+tanl), then compute tan(A) as:                          *
			//                                                                      *
			//                                t ( C+T)                              *
			//                tan(A)= T   +  -----------                            *
			//                                  C-t                                 *
			//                                                                      *
			//         numerator represented as (p19,p21,p22) in the code below.    *
			//                                                                      *
			//***********************************************************************
			p23 = __FSUB( c1, tan );							// High order part of C-t
			p24 = __FSUB( c2, tanl );
			rec = 1.0/p23;
			p25 = c1 - p23 - tan;
			frac = p19*rec;							// start computation of fraction
			
			if (fabs(c2) > fabs(tanl))
				p26 = c2 - p24 - tanl;
			else
				p26 = c2 - (p24+tanl);
			
			p27 = __FADD( p25, p24 );
			p28 = p25 - p27 + p24 + p26;													// C-t=(p23,p27,p28) done
			frac = frac + __FNMSUB( p23, frac, p19 )*rec; // frac + (p19-p23*frac)*rec;		// refine lead fraction wd.
			p29 = __FNMSUB( p23, frac, p19 );			//p19 - p23*frac;
			p30 = frac*p27;																	// p30 a negative term
			p31 = __FNMSUB( frac, p27, p30 ) + p22;		// p30 - frac*p27 + p22;
			p32 = __FNMSUB( frac, p28, p31 );			// p31 - frac*p28;
			p33 = __FSUB( p29, p30 );
			p35 = p33 + p21;
			frac2 = p35*rec;
			frac2 = frac2 + __FNMSUB( frac2, p23, p35)*rec; // frac2 + (p35-frac2*p23)*rec;	// second fraction wd.
			
			if (fabs(p29) > fabs(p30))
				p34 = p29 - p33 - p30 + p32;
			else
				p34 = p29 - (p33+p30) + p32;
			
			p36 = __FNMSUB( frac2, p23, p35 ) + __FNMSUB( frac2, p27, p34 );	// p35 - frac2*p23 + p34 - frac2*p27;
			
			if (fabs(p33) > fabs(p21))
				p37 = p33 - p35 + p21;
			else
				p37 = p21 - p35 + p33;
			
			p38 = p37 + p36;
			frac3 = p38*rec;						// Fraction of above formula done
			
			//***********************************************************************
			//        Finally, add (tan,tanl) and (frac,frac2,frac3) for result     *
			//***********************************************************************
			
			p39 = __FADD( t1, frac );
			p40 = t1 - p39 + frac;
			p41 = __FADD( t2, frac2 );
			
			if (fabs(t2) > fabs(frac2))
				p42 = t2 - p41 + frac2 + frac3;
			else
				p42 = frac2 - p41 + t2 + frac3;
			
			p43 = __FADD( p40, p41 );
			p44 = p40 - p43 + p41 + p42;
			p45 = __FADD( p39, p43 );
			p46 = p39 - p45 + p43 + p44;
			
			// getting a cannonical result
			
			ddx.d[0] = __FADD( p45, p46 ); 			
			ddx.d[1] = (p45 - ddx.d[0]) + p46; 
			
			FESETENVD(fenv);
			return (ddx.f);
		}
		else {
			
			//***********************************************************************
			//      The remaining case is for argments in the neighborhood of       *
			//      even multiples of pi/2.                                         *
			//      Use the identity                                                *
			//                          tan(A+(2k+1)pi/2)=-cot(A)                   *
			//                                                                      *
			//      Let the reduced argument A be written as                        *
			//                     A=x+y, where  x= TanTablePtr[index].arg          *
			//      Then if T=tan(x) and C=cot(x) as read from tantbl, and t=tan(y),*
			//      (computed above as (tan+tanl), compute cot(A) as:               *
			//                                                                      *
			//                                  t(C+T)                              *
			//                      cot(A)=C - ----------                           *
			//                                  T+t                                 *
			//***********************************************************************
			p23 = __FADD( t1, tan );							// High order part of T+t
			p24 = __FADD( t2, tanl );
			rec = 1.0/p23;
			p25 = t1 - p23 + tan;
			frac = p19*rec;							// start computation of fraction
			
			if (fabs(t2) > fabs(tanl))
				p26 = t2 - p24 + tanl;
			else
				p26 = tanl - p24 + t2;
			
			p27 = __FADD( p25, p24 );
			p28 = p25 - p27 + p24 + p26;													// T+t= (p23,p27,p28) done
			frac = frac + __FNMSUB( p23, frac, p19 )*rec; // frac + (p19-p23*frac)*rec;		// refine lead fraction wd.
			p29 = __FNMSUB( p23, frac, p19 );		// p19 - p23*frac;
			p30 = frac*p27;																	// p30 a negative term
			p31 = __FNMSUB( frac, p27, p30 ) + p22; // p30 - frac*p27 + p22;
			p32 = __FNMSUB( frac, p28, p31 );		// p31 - frac*p28;
			p33 = __FSUB( p29, p30 );
			p35 = __FADD( p33, p21 );
			frac2 = p35*rec;
			frac2 = frac2, __FNMSUB( frac2, p23, p35 )*rec; // frac2 + (p35-frac2*p23)*rec;	// second fraction wd.
			
			if (fabs(p29) > fabs(p30))
				p34 = p29 - p33 - p30 + p32;
			else
				p34 = p29 - (p33+p30) + p32;
			
			p36 = __FNMSUB( frac2, p23, p35 ) + __FNMSUB( frac2, p27, p34 ); // p35 - frac2*p23 + p34 - frac2*p27;
			
			if (fabs(p33) > fabs(p21))
				p37 = p33 - p35 + p21;
			else
				p37 = p21 - p35 + p33;
			
			p38 = p37 + p36;
			frac3 = p38*rec;						// Fraction done
			
			//*************************************************************************
			//    At last, subtract (cot,cotl) from (frac,frac2,frac3) for the result *
			//*************************************************************************
			p39 = __FSUB( frac, c1 );
			p40 = frac - (c1+p39);
			p41 = __FSUB( frac2, c2 );
			if (fabs(c2) > fabs(frac2))
				p42 = frac2 - (c2+p41) + frac3;
			else
				p42 = frac2 - p41 - c2 + frac3;
			
			p43 = __FADD( p40, p41 );
			p44 = p40 - p43 + p41 + p42;
			p45 = __FADD( p39, p43 );
			p46 = p39 - p45 + p43 + p44;
			
			//getting a cannonical result
			
			ddx.d[0] = __FADD( p45, p46 ); 			
			ddx.d[1] = (p45 - ddx.d[0]) + p46; 
			
			FESETENVD(fenv);
			return (ddx.f);
		}
	}
}
#endif