ArcTanDD.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@
 */
//
//	ArcTanDD.c
//	
//	Double-double Function Library
//	Copyright: © 1995-96 by Apple Computer, Inc., all rights reserved
//	
//	long double atanl( long double x );
//	long double atan2l( long double y, long double x );
//
//	Change History:
//		8/31/96	- PAF - Changed if() to else (see Note 1)
//

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

static const double kBig = 6.755399441055744e+15;				// 0x1.8p52   

static const double kPiBy2 = 1.570796326794896619231322;		// 0x1.921FB54442D18p0   
static const double kPiBy2Tail1  =  6.123233995736766036e-17;	// 0x1.1A62633145C07p-54  
static const double kPiBy2Tail1m  = 6.123233995736764803e-17;	// 0x1.1A62633145C06p-54  
static const double kPiBy2Tail2 = -1.4973849048591698e-33;		// 0x1.F1976B7ED8FBCp-110  
static const double kRecip = 6.869637070016756334588131e-12;	// 0x1.e368613a3b430p-38
static const double kLimit1 = 0.0859375;						// 0x1.6p-4   
static const double kLimit2 = 1.015625;	
static const double kLimit3 = 13.0;	
static const double kLimit4 =  8.11296384146066817e+31;			// 0x1.0p+106 
static const Float  kInfinityu = {0x7f800000};

//***************************************************************************
//    factor -> gives the number of table values between the consecutive    *
//                  integers 1 to 12                                        *
//***************************************************************************

static const float factor[12] = {
	32.0f,16.0f,8.0f,4.0f,4.0f,2.0f,2.0f,1.0f,1.0f,1.0f,1.0f,1.0f
};


//***************************************************************************
//    Adjust     -> index into accurate atan table, for integer arguments   *
//                  from 1 to 12.                                           *
//***************************************************************************

static const float adjust[12] = {
	32.0f,64.0f,88.0f,104.0f,104.0f,116.0f,116.0f,124.0f,124.0f,124.0f,124.0f,124.0f
};


//***************************************************************************
//         atancoeff: coefficients of atan approximating polynomial,        *
//                    multiplied by 145568097675.0                          *
//***************************************************************************

static const double atancoeff[15] = {
	145568097675.0, -48522699225.0,
	 29113619535.0, -20795442525.0,   16174233075.0,
	-13233463425.0,  11197545975.0,   -9704539845.0,
	  8562829275.0,  -7661478825.0,    6931814175.0,
	 -6329047725.0,   5822723907.0,   -5391411025.0,
	  5019589575.0
};

static const uint32_t atanTable[] __attribute__ ((aligned(8))) = {

	//*******************************************************************
	//	  Accurate table over range (	.093750,   .984375).            *
	//*******************************************************************
	
	0x3FB80000,	0x96AC4A5A,	0x3FB7EE18,	0xBB5F2BDC,	0x3B1F7151,	0xA60EEDFD,
	0x3FBC0000,	0xAB2E5990, 0x3FBBE39F,	0x679755D8, 0x3B5AF7C1,	0xBBB2A139,
	0x3FC00000,	0x7F748304, 0x3FBFD5BB,	0x95A94110, 0x3B5B0A0C,	0x811C880E,
	0x3FC20000,	0x3DF2157F, 0x3FC1E1FB,	0x37C2C7E5, 0x3B617207,	0xB776AF4F,
	0x3FC40000,	0x30754446, 0x3FC3D6EF,	0x1814018D, 0xBB63A092,	0x88503B16,
	0x3FC60000,	0x7996ECD4, 0x3FC5C981,	0x94588C25, 0xBB5B6979,	0x558D8B64,
	0x3FC80000,	0x76D31B22, 0x3FC7B97B,	0xBE985C07, 0x3B6708AD,	0xCC1C4519,
	0x3FCA0000,	0xEADCE47A, 0x3FC9A6A9,	0xCAFAF9A8, 0x3B62971C,	0xA562CA82,
	0x3FCC0000,	0x7E8C3F27, 0x3FCB90D7,	0xCB57348B, 0xBB493C17,	0xA875159C,
	0x3FCE0000,	0x5B063A6A, 0x3FCD77D6,	0x3569311A, 0xBB5A33AB,	0xA192F16A,
	0x3FD00000,	0x64EF22FF, 0x3FCF5B76,	0xB72AE095, 0x3B52E368,	0xB50909D6,
	0x3FD10000,	0x87DE9AAA, 0x3FD09DC6,	0x16C29778, 0xBB55B0DC,	0xF3418D68,
	0x3FD20000,	0x0AAD7C84, 0x3FD18BF5,	0xACF10E71, 0x3B41EAC1,	0x346154B4,
	0x3FD30000,	0x1D92088F, 0x3FD27837,	0x3B84D2FF, 0x3B48BBE9,	0xDC9998E5,
	0x3FD40000,	0xCC4E6285, 0x3FD36277,	0xF12910F6, 0xBB6585F3,	0xFA8B319B,
	0x3FD50000,	0x766A7592, 0x3FD44AA4,	0xA1AA8D78, 0x3B6C6289,	0x13F656CF,
	0x3FD60000,	0xCFA33544, 0x3FD530AE,	0x5303BB5D, 0x3B63F533,	0xC1117E20,
	0x3FD70000,	0xC925EED6, 0x3FD61484,	0xB52DF321, 0xBB70059A,	0x66E93FA9,
	0x3FD80000,	0x54AA10D9, 0x3FD6F619,	0x8C1ECA84, 0x3B60C60C,	0x08BED883,
	0x3FD90000,	0x113AC5B7, 0x3FD7D560,	0x5A568BA4, 0xBB541817,	0xB5654BFC,
	0x3FDA0000,	0x6D101F2E, 0x3FD8B24D,	0x96E71249, 0xBB55A721,	0x1379C827,
	0x3FDB0000,	0xE246C017, 0x3FD98CD6,	0x05641F92, 0x3B6918FF,	0xCA4E1DF9,
	0x3FDC0000,	0x39B58346, 0x3FDA64EE,	0xF43C341A, 0x3B6633D0,	0xEECC747A,
	0x3FDD0000,	0x4C36F71C, 0x3FDB3A91,	0x5CE1B431, 0xBB61E754,	0x0B1A4FB4,
	0x3FDE0000,	0x66CA921B, 0x3FDC0DB5,	0x1D94F195, 0x3B571DE1,	0xE6981964,
	0x3FDF0000,	0x3AD3639E, 0x3FDCDE53,	0x72D1AC94, 0xBB736718,	0x370B3962,
	0x3FE00000,	0x7815CCCD, 0x3FDDAC67,	0xC5849B77, 0x3B52FF7C,	0xC7D84489,
	0x3FE08000,	0xF40DDDFC, 0x3FDE77ED,	0x00AEBFCD, 0x3B62E958,	0xE33FF707,
	0x3FE10000,	0xA6C93DAD, 0x3FDF40DE,	0x0F7AA915, 0x3B77641D,	0x15CA0475,
	0x3FE18000,	0x85EBF4D7, 0x3FE0039C,	0xDAD8BEF9, 0x3B722996,	0x8391BF09,
	0x3FE20000,	0xC129DA76, 0x3FE0657F,	0x27977717, 0xBB761869,	0xC045B9BD,
	0x3FE28000,	0x93F751BC, 0x3FE0C614,	0xCA41B0AE, 0xBB7A8A92,	0x3C961154,
	0x3FE30000,	0x7C1EF6CE, 0x3FE1255D,	0xF7C0A55C, 0x3B591CBF,	0x478B1C38,
	0x3FE38000,	0x169FA349, 0x3FE1835A,	0x993DD546, 0x3B6FBFAB,	0x698175A0,
	0x3FE40000,	0x27D6B511, 0x3FE1E00B,	0xC884E585, 0xBB248C8D,	0x9BB39A1C,
	0x3FE48000,	0x6BE89876, 0x3FE23B72,	0x2F4EF788, 0xBB815690,	0x7DF099CF,
	0x3FE50000,	0x280F4ABE, 0x3FE2958E,	0x7530C25F, 0xBB72B7E0,	0xEE07EAE7,
	0x3FE58000,	0x2FEDFD47, 0x3FE2EE62,	0xA50C9A07, 0xBB2B522D,	0xBFFA3A8C,
	0x3FE60000,	0x875C63BA, 0x3FE345F0,	0x78B8C001, 0xBB408ED6,	0xD056E4C4,
	0x3FE68000,	0x8A368DBA, 0x3FE39C39,	0x79511750, 0x3B77CB64,	0x0C777681,
	0x3FE70000,	0xEE7CA285, 0x3FE3F140,	0x55DECB32, 0xBB5AA9DC,	0xBC68CF79,
	0x3FE78000,	0xA4904A32, 0x3FE44506,	0xC661B510, 0x3B7EBCE3,	0x3EB6C0CD,  
	0x3FE80000,	0x4CA9F6E2, 0x3FE4978F,	0xD4373CAA, 0x3AFF9675,	0xAB73F345,
	0x3FE88000,	0xD7B95553, 0x3FE4E8DE,	0xE3B770C7, 0xBB715BD3,	0xBFA6C0AD,
	0x3FE90000,	0x849C7E04, 0x3FE538F5,	0xCDE27023, 0xBB6B45C2,	0x3E5A5CC8,
	0x3FE98000,	0xC171C662, 0x3FE587D8,	0x95C38C2E, 0x3B7887A1,	0x345B0752,
	0x3FEA0000,	0x4321104C, 0x3FE5D589,	0xAF86142F, 0xBB8438EB,	0x1AE94AC1,
	0x3FEA8000,	0x838D8D26, 0x3FE6220D,	0x5F66C783, 0x3B7CA422,	0x8F96E72F,
	0x3FEB0000,	0xE4DA897E, 0x3FE66D66,	0xBED2B224, 0xBB35FB4D,	0xFD76AD75,
	0x3FEB8000,	0xD775765C, 0x3FE6B799,	0x0DF9D017, 0x3B7B5383,	0xE0D80102,
	0x3FEC0000,	0xAA053CA7, 0x3FE700A8,	0x25C3BB63, 0xBB7D5234,	0x29157DFC,
	0x3FEC8000,	0x773AAC9B, 0x3FE74897,	0xD237C658, 0x3B6ACFB2,	0xF4E46A15,
	0x3FED0000,	0x3FA5ABCF, 0x3FE78F6B,	0xE04F6C27, 0xBB8680F4,	0x55D1C7DF,
	0x3FED8000,	0xD31E46AA, 0x3FE7D528,	0x9AC0235F, 0xBB61A1B0,	0x1EABC199,
	0x3FEE0000,	0xBB676A41, 0x3FE819D1,	0x1AD33A87, 0xBB7726EA,	0x735248A8,
	0x3FEE8000,	0x76CB5B54, 0x3FE85D69,	0x95ABE3CC, 0x3B772D9C,	0x044A3CEB,
	0x3FEF0000,	0x26506206, 0x3FE89FF6,	0x131BC92E, 0x3B860889,	0x51E5F051,
	0x3FEF8000,	0xED69BF50, 0x3FE8E17B,	0x2230274E, 0xBB7FB7C5,	0x5CA2A11A, 
	
	//*******************************************************************
	//	  Accurate table over range (  1.000000,  1.968750).            *
	//*******************************************************************
	
	0x3FF00000,	0x38A5805B, 0x3FE921FB,	0x8CE9AD0F, 0xBB728518,	0x56659FDA,
	0x3FF08000,	0xE4CDB94C, 0x3FE9A001,	0x86F9991F, 0xBB7F4F21,	0xC242FD58,
	0x3FF10000,	0x29DA60EA, 0x3FEA1A26,	0x1A19C30A, 0xBB87F5DE,	0xE73B7B3A,
	0x3FF18000,	0x9BBB03C5, 0x3FEA908B,	0x882B14EE, 0x3B8A364F,	0xC319894C,
	0x3FF20000,	0x42A338B5, 0x3FEB034F,	0x7337C8D1, 0xBB768B38,	0xD9BDAC43,
	0x3FF28000,	0xED0428FC, 0x3FEB7292,	0x7FBAC8B6, 0x3B7B400A,	0x53FD9AAF,
	0x3FF30000,	0x36946784, 0x3FEBDE71,	0x1A8E3A61, 0x3B6C2B37,	0x13D123E1,
	0x3FF38000,	0x21F16647, 0x3FEC470A,	0xDA7DBAA7, 0xBB523353,	0xD7A5886A,
	0x3FF40000,	0xF06227B4, 0x3FECAC7D,	0x132231D2, 0xBB8AEFE3,	0x0558B34B,
	0x3FF48000,	0xF9AB1CF7, 0x3FED0EE2,	0xE23FB012, 0xBB4465DC,	0xD4A6C05D,
	0x3FF50000,	0x1DAE2F90, 0x3FED6E57,	0xE51C7E16, 0xBB85A9B9,	0x83503D52,
	0x3FF58000,	0xD679A51F, 0x3FEDCAF8,	0xC6A4C801, 0x3B8576BE,	0xCEAEFC33,
	0x3FF60000,	0xCE367C5E, 0x3FEE24DD,	0xD375A188, 0xBB78F5B1,	0x1E4802A0,
	0x3FF68000,	0xCD05D052, 0x3FEE7C20,	0xCBEB8A43, 0xBB8AD216,	0x43D5F156,
	0x3FF70000,	0x98AA0697, 0x3FEED0D9,	0xE022B144, 0x3B267908,	0xB6C976B5,
	0x3FF78000,	0xA3FB1493, 0x3FEF2320,	0xDB8F1029, 0xBB85B081,	0x83128CFA,
	0x3FF80000,	0x681E0063, 0x3FEF730C,	0x12946C3F, 0xBB8AA69F,	0xBF359098,
	0x3FF88000,	0x75E901B6, 0x3FEFC0B1,	0xB86DE15F, 0x3B8EAA0C,	0x4663839A,
	0x3FF90000,	0x736A4590, 0x3FF00613,	0x4FBE5B7E, 0xBB4E3551,	0x42306F18,
	0x3FF98000,	0xC5B189FC, 0x3FF02ABF,	0xA107BE50, 0x3B63FD1C,	0xCE1D6AEB,
	0x3FFA0000,	0x41411D33, 0x3FF04E67,	0x3966892F, 0xBB86ED08,	0x00B1EF7B,
	0x3FFA8000,	0x4DDFE3DA, 0x3FF07113,	0xDB772C28, 0x3B746CEB,	0x807EA687,
	0x3FFB0000,	0xA7AB4BE2, 0x3FF092CE,	0x72AC060E, 0xBB72F0AF,	0x054AD94E,
	0x3FFB8000,	0x35EE2E87, 0x3FF0B39F,	0x5C6DBFD7, 0xBB8A699D,	0x14CF8D89,
	0x3FFC0000,	0x79493CCE, 0x3FF0D38F,	0x4A3683E2, 0xBB7C178E,	0x4E0FB167,
	0x3FFC8000,	0x7C2C16CD, 0x3FF0F2A5,	0xF7C1C52E, 0xBB77BBEC,	0x0F5FF2D9,
	0x3FFD0000,	0xF7925B8C, 0x3FF110EB,	0x3A456E0E, 0x3B86B850,	0xB94665FD,
	0x3FFD8000,	0xD4EDAFAA, 0x3FF12E66,	0x2A98E0EA, 0x3B8FC04D,	0xF88FF47B,
	0x3FFE0000,	0x1822F5F9, 0x3FF14B1D,	0xDB5166C5, 0x3B400FEE,	0xA21D9DDA,
	0x3FFE8000,	0x5F375D4F, 0x3FF16719,	0x6EAC8C79, 0xBB45D7BA,	0x99C8508E,
	0x3FFF0000,	0x983AF36D, 0x3FF1825F,	0x2745DBAE, 0xBB8849C2,	0x328CE199,
	0x3FFF8000,	0xDF7044EE, 0x3FF19CF5,	0x48D90CBC, 0xBB876AC4,	0x23C76E06,
	
	//******************************************************************
	//	  Accurate table over range (  2.000000,  2.937500).           *
	//******************************************************************
	
	0x40000000,	0xC89A22AD, 0x3FF1B6E1,	0xE3296298, 0xBB83E57A,	0x4E33C68C,
	0x40008000,	0xAE85BCDD, 0x3FF1E8D4,	0xB635433E, 0x3B8FDE86,	0xD38237A1,
	0x40010000,	0x1E66B4AD, 0x3FF21862,	0xFF00ECE6, 0x3B74E47A,	0x940D38B0,
	0x40018000,	0x257005B2, 0x3FF245B5,	0x07E38198, 0xBB90341D,	0xF7AEB5A2,
	0x40020000,	0x55656B58, 0x3FF270EF,	0x71D13D73, 0x3B88C680,	0x13906BD2,
	0x40028000,	0x466225D9, 0x3FF29A34,	0x0FA8F5E6, 0x3B91F854,	0xF983041D,
	0x40030000,	0x717EC3D3, 0x3FF2C1A2,	0x640509A5, 0x3B8E1C75,	0xAF6F1BFF,
	0x40038000,	0x75A7B542, 0x3FF2E757,	0x4A697F72, 0x3B92AB3B,	0x412F53FB,
	0x40040000,	0x00512EF1, 0x3FF30B6D,	0x7980B2E2, 0x3B8AD702,	0x40329E0C,
	0x40048000,	0x09ED9B8A, 0x3FF32DFE,	0x0460EC6D, 0x3B6B92BF,	0x54B284E9,
	0x40050000,	0x1B9F04CA, 0x3FF34F1F,	0xC21A2D1B, 0xBB82C1A6,	0xD9CE49BA,
	0x40058000,	0x697C1A2B, 0x3FF36EE8,	0x0C4A7DD1, 0x3B73DD4E,	0x5C021394,
	0x40060000,	0xABB8A398, 0x3FF38D6A,	0x94FD5FEE, 0x3B373349,	0x3B4FC3A6,
	0x40068000,	0x18374FED, 0x3FF3AAB9,	0x8BB1767A, 0x3B72834E,	0xCBC4CA13,
	0x40070000,	0xBCD94FDA, 0x3FF3C6E6,	0x7976E8D9, 0xBB931322,	0x7228DCD0,
	0x40078000,	0x7BA476B6, 0x3FF3E200,	0xC84E84B4, 0xBB87D21A,	0xC44FB391,
	
	//******************************************************************
	//	  Accurate table over range (  3.000000,  3.875000).           *
	//******************************************************************
	
	0x40080000,	0xD7180754, 0x3FF3FC17,	0x967F5249, 0xBB64B771,	0x5FA0957B,
	0x40090000,	0x6DFDCD42, 0x3FF42D70,	0x558EAD29, 0x3B84B923,	0xDD410F01,
	0x400A0000,	0xE7C427AE, 0x3FF45B54,	0xAB8A2C51, 0x3B91D694,	0xE38B7BFF,
	0x400B0000,	0x10E3223E, 0x3FF4861B,	0x4FB5B5B4, 0xBB7CB027,	0x50BA194A,
	0x400C0000,	0x8EA13F5E, 0x3FF4AE11,	0x11ECF83F, 0x3B8B6616,	0xD89A7B16,
	0x400D0000,	0x06D105F7, 0x3FF4D378,	0xC2906964, 0x3B93659C,	0x0737481A,
	0x400E0000,	0x727DD5BC, 0x3FF4F68D,	0xF99AE908, 0x3B8514AE,	0x4636F80B,
	0x400F0000,	0x3FDFD1C5, 0x3FF51785,	0x020F4064, 0xBB65DB78,	0x8577ADC5,
	
	//******************************************************************
	//    Accurate table over range (  4.000000,  5.750000).           *
	//******************************************************************
	
	0x40100000,	0xCF49E240, 0x3FF5368C,	0xC5E4B1C8, 0xBB5E5499,	0xA9C4BF51,
	0x40110000,	0x94E90733, 0x3FF56F6F,	0x52E31047, 0x3B8F3B21,	0x724AC4D7,
	0x40120000,	0x46578102, 0x3FF5A250,	0x5F4EF405, 0x3B8657A1,	0x1EBE66F0,
	0x40130000,	0xCA3D6D0F, 0x3FF5D013,	0xE66FF8CA, 0x3B60B423,	0x26CAA3CA,
	0x40140000,	0x32964FAE, 0x3FF5F973,	0x1CEDA34B, 0x3B635D66,	0xC3ED1C18,
	0x40150000,	0xA6AA4915, 0x3FF61F06,	0xDE005085, 0x3B6E8C7F,	0x000B0CDD,
	0x40160000,	0x8906EAE3, 0x3FF6414D,	0x5593660A, 0x3B95FADD,	0xB8BF7A98,
	0x40170000,	0x831FA60C, 0x3FF660B0,	0x3BD94D47, 0x3B96491D,	0x5B9C83F6,
	
	//******************************************************************
	//	  Accurate table over range (  6.000000,  7.500000).           *
	//******************************************************************
	
	0x40180000,	0x8DCF9463, 0x3FF67D88,	0x73114F7D, 0x3B928119,	0xC1230051,
	0x401A0000,	0xCBB9EEAD, 0x3FF6B0BA,	0xFB083C0D, 0x3B93280F,	0xAE0D69E3,
	0x401C0000,	0xDEF559E1, 0x3FF6DCC5,	0x8D8B9598, 0x3B86722D,	0x1022C04A,
	0x401E0000,	0x744BF1EC, 0x3FF7030D,	0x01605442, 0xBB736E8B,	0xF92A87CA,
	
	//******************************************************************
	//    Accurate table over range (  8.000000, 12.000000).           *
	//******************************************************************
	
	0x40200000,	0x456D4528, 0x3FF7249F,	0xB324E4B7, 0x3B675690,	0xDD79F342,
	0x40220000,	0x43E2466D, 0x3FF75CBA,	0xD9437CC1, 0x3B46FCCD,	0xEB45ACCD,
	0x40240000,	0x18A681EA, 0x3FF789BD,	0x2E09D7EA, 0x3B6C324D,	0xB28DD9D9,
	0x40260000,	0xF3FE273A, 0x3FF7AEA3,	0x9C1AAC21, 0x3B8D6F66,	0xFC2D1461,
	0x40280000,	0x455A6241, 0x3FF7CD6F,	0x71992B07, 0x3B8B2FE8,	0xBAEBACBC,
	0x40280000,	0x455A6241, 0x3FF7CD6F,	0x71992B07, 0x3B8B2FE8,	0xBAEBACBC
};

struct atantblEntry { 
   double One;				
   double Two;
   double Three;
} atantblEntry;
	
long double atanl( long double x )
{
	double p1, p2, p3, p2a, p4, p5, p6, p7, p8, p29, p30, p33, p35;
	double quot, quot2;
	double rarg, result, restemp;
	double top, top2, topextra, top3;
	double bot, bot2;
	double r, rarg2, rarg2a;
	double temps;
	double remainder, residual;
	double partial;  
	double arg, argl, argsq, argsq2;
	double temp,temp2; 
	double sum, sumt, suml;
	double absx;
	double xHead, xTail;
	double res, reslow, resmid, restiny, reshi, resbot;
	double p, q; 
	double prod, prodlow;	
	Double DeN ;
	double fenv;
	DblDbl ddx;
	int i, rflag;
	uint32_t index;
	
	struct atantblEntry *atanPtr = (struct atantblEntry *)atanTable - 6;
	
        FEGETENVD(fenv);
        FESETENVD(0.0);
	
	ddx.f = x;
	xHead = ddx.d[0];
	xTail = ddx.d[1];
	
	// x is a NaN or zero
	
	if ((xHead != xHead)||(xHead == 0.0)) {			// x is 0 or a NaN
                FESETENVD(fenv);
		return x;
	}
	
	// x is not a NaN nor zero
	
	absx = fabs(xHead);
	
	//finite result
	
	//***************************************************************************
	//      This routine breaks the arguments into five ranges:                 *
	//      Range 1:          0 <= |x| <= .0859375                              *
	//      Range 2:   .0859375 <= |x| <= 1.015625                              *
	//      Range 3:   1.015625 <= |x| < 13                                     *
	//      Range 4:         13 <= |x| < 2^107                                  *
	//      Range 5:      2^107 <= |x| <= Infinity                              *
	//                                                                          *
	//         >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>         *
	//                                                                          *
	// In Range 1, a 15th degree polynomial in x^2 is used to evaluate QATAN.   *
	//                                                                          *
	// In Range 2  and,                                                         *
	// In Range 3  the accurate table method is used to range-reduce the        *
	//             argument to a small value to be used with an interpolating   *
	//             polynomial for QATAN.  In Ranges 2 and 3, the resultant      *
	//             argument is less than 2^-7, so an 8th degree polynomial can  *
	//             be used.  Range 3 requires additional work to find the       *
	//             correct index into the accurate table.                       *
	//                                                                          *
	// In Range 4, QATAN(x) = Pi/2 - QATAN(1/x).  After taking the reciprocal   *
	//             of the argument, the quotient is less than 1/13, so that the *
	//             15th degree polynomial used for Range 1 is used to compute   *
	//             QATAN(1/x)                                                   *
	//                                                                          *
	// In Range 5, the correctly rounded value differs from Pi/2 by less than   *
	//             1/2 ulp.  It is, however, incorrect to let the value become  *
	//             correctly rounded, as that correctly rounded value exceeds   *
	//             PI/2 and is in the wrong quadrant.  Applying the inverse     *
	//             function, QTAN, to the correctly rounded value of PI/2 would *
	//             produce a result of the opposite sign to the argument        *
	//             originally given to ATAN.                                    *
	//***************************************************************************
	
	// Range 5
	
	if (absx >= kLimit4) {							// Filter out unusual arguments
		ddx.d[0] = kPiBy2;  
		ddx.d[1] = kPiBy2Tail1m;					// Intentional misround to
		if (xHead < 0.0) {      
			ddx.d[0] = -ddx.d[0];  
			ddx.d[1] = -ddx.d[1];		
		}
                FESETENVD(fenv);
		return (ddx.f);
	}
	
	// Ranges 1 and 4 (DO NOT USE TABLE)
	
	rflag = 0;	
	
	// Range 4
	
	if (absx >= kLimit3) {							// Range 4
		rarg = 1.0/xHead;                         
		rflag = 1;									// indicate large argument
		p1 = __FNMSUB( xHead, rarg, 1.0 );			// 1.0 - xHead*rarg;						// time to complete hi precision
		p2 = xTail*rarg;							// reciprocal of argument, to
		p3 = p1 - p2;								// use as reduced argument.
		rarg2a = p3*rarg;
		rarg2 = rarg2a + __FNMSUB( xHead, rarg2a, p3 )*rarg; // rarg2a + (p3 - xHead*rarg2a)*rarg;
		xHead = rarg;
		xTail = rarg2;
	}
	
	// Ranges 1 and 4 only
	
	if( (absx < kLimit1) || (absx >= kLimit3) ) {	// Code for Range 1 and Range 4
		temp = 2.0*xHead*xTail;						// direct evaluation of series
		argsq = __FMADD( xHead, xHead, temp );		// xHead*xHead + temp;					// compute square of argument
		argsq2 = __FMADD( xTail, xTail, __FMSUB( xHead, xHead, argsq ) + temp ); // xHead*xHead - argsq + temp + xTail*xTail;
		sum = atancoeff[14];
		sum = __FMADD( argsq, sum, atancoeff[13] );			// atancoeff[13] + argsq*sum;
		sum = __FMADD( argsq, sum, atancoeff[12] );			// atancoeff[12] + argsq*sum;
		sum = __FMADD( argsq, sum, atancoeff[11] );			// atancoeff[11] + argsq*sum;
		sumt = __FMADD( argsq, sum, atancoeff[10] );		// atancoeff[10] + argsq*sum;
		sum = __FMADD( argsq, sum, atancoeff[9] );			// atancoeff[9] + argsq*sumt;
		suml = __FMADD( argsq, sumt, atancoeff[9] - sum );	// atancoeff[9] - sum + argsq*sumt;
		
		for (i = 8; i >= 1; i--) {	 		   
			temps = __FMADD( argsq, sum, atancoeff[i] );	// atancoeff[i] + sum*argsq;
			/* suml = atancoeff[i] - temps+sum*argsq + (sum*argsq2 + suml*argsq); */
			suml = __FMADD( sum, argsq, atancoeff[i] -  temps ) + __FMADD( sum, argsq2, suml*argsq );
			sum=temps;
		}
		
		prodlow = __FMADD( suml, argsq, sum*argsq2 );		// suml*argsq + sum*argsq2;			// multiply by arg^2
		prod = __FMADD( sum, argsq, prodlow );				// sum*argsq + prodlow;
		prodlow = __FMSUB( sum, argsq, prod ) + prodlow;	// sum*argsq - prod + prodlow;
		temp2 = __FMADD( prodlow, xHead, prod*xTail );		// prodlow*xHead + prod*xTail;		// multiply by arg
		temp = __FMADD( prod, xHead, temp2 );				// prod*xHead + temp2;
		temp2 = __FMSUB( prod, xHead, temp ) + temp2;		// prod*xHead - temp + temp2;
		sum = temp*kRecip;																		// sum of series ---
		remainder = __FNMSUB( sum, atancoeff[0], temp );	// (temp - sum*atancoeff[0]);
		partial = __FADD( remainder, temp2 );
		residual = remainder - partial + temp2;
		suml = __FMADD( partial, kRecip, (residual*kRecip) );// partial*kRecip + (residual*kRecip);
		res = __FADD( xHead, sum );																// except for argument
		reslow = (xHead - res + sum);															// exact
		resmid = __FADD( xTail, suml );
		restiny = xTail - resmid + suml;
		p = __FADD( reslow, resmid );															// sum of middle terms
		q = reslow - p + resmid;																// possible residual
		reshi = res + p;
		resbot = (res - reshi + p) + (q + restiny);
		if (rflag == 1) {							// large magnitude argument
			if (xHead > 0) {						// fetch pi/2 with proper sign
				p1 = kPiBy2;
				p2 = kPiBy2Tail1;
				p2a = kPiBy2Tail2;
			}
			else {
				p1 = -kPiBy2;
				p2 = -kPiBy2Tail1;
				p2a = -kPiBy2Tail2;
			}
			p3 = __FSUB( p1, reshi );				// subtract result from pi/2
			p4 = p1 - p3 - reshi;
			p5 = __FSUB( p2, resbot );
			p6 = p2 - p5 - resbot + p2a;
			p7 = __FADD( p4, p5 );
			reshi = __FADD( p3, p7 );
			if (fabs(p4) > fabs(p5)) 
				p8 = p4 - p7 + p5;
			else
				p8 = p5 - p7 + p4;
			
			resbot = p3 - reshi + p7 + (p6 + p8 - restiny);
		}
		ddx.d[0] = __FADD( reshi, resbot );  
		ddx.d[1] = (reshi - ddx.d[0]) + resbot;		
		
                FESETENVD(fenv);
		return (ddx.f);
	}
	
	// Ranges 2 and 3  (USE TABLE)
	
	//if( (absx >= kLimit1) || (absx < kLimit3) ) {	// Code for Range 2 and Range 3
	
	// Note 1: The above statement in the Taligent code always tests true.  Logically,
	// it was probably intended that || be &&.  This would make the test logically
	// equivalent to an "else" clause for the previous if.  Replacing with an else
	// also eliminates the problem of a code path that avoids any return statements
	// (as detected by MrC as a compiler warning).  PAF
	else {
		
		DeN.f = __FMADD( 64.0, absx, kBig ); // 64.0*absx + kBig;
		
		if (absx > kLimit2) {						// for large args, use recip.  Range 3
			FESETENVD(kBig+1.0);						// use trunc mode momentarily
			DeN.f = absx + kBig;					// determine range of argument
			FESETENVD(kBig+1.0);						// use trunc mode momentarily
			DeN.f = absx + kBig;					// determine range of argument
			FESETENVD(kBig);							// return to round-to-nearest
			index = DeN.i[1];
			i = index - 1;
			DeN.f = __FMADD( factor[i], absx, kBig + adjust[i] ); // factor[i]*absx + kBig + adjust[i];
		}
		
		index = DeN.i[1];
		
		if (xHead < 0) {							// use abs. value
			arg  = -xHead;
			argl = -xTail;
		}
		else {
			arg  = xHead;
			argl = xTail;
		}
		
		top = __FSUB( arg, atanPtr[(int32_t)index].One );
		top2 = __FADD( top, argl );
		topextra = arg - top - atanPtr[(int32_t)index].One;
		top3 = top - top2 + argl + topextra;												// reduced numerator
		bot =__FMADD( arg, atanPtr[(int32_t)index].One, 1.0 ); // 1.0 + arg*atanPtr[(int32_t)index].One;
		/* bot2 = 1.0 - bot + arg*atanPtr[(int32_t)index].One +
			argl*atanPtr[(int32_t)index].One;	*/											// denominator
		bot2 = __FMADD( argl, atanPtr[(int32_t)index].One, __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 - bot ) );
		r = 1.0/bot;
		quot = r*top2;
		quot = __FMADD( __FNMSUB( bot, quot, top2), r, quot );	// quot + (top2 - bot*quot)*r;
		
		p29 = __FNMSUB( bot, quot, top2);				// top2 - quot*bot;
		p30 = quot*bot2;																	// negative term
		p33 = p29 - p30;
		p35 = p33 + top3;
		quot2 = p35*r;
		quot2 = __FMADD( __FNMSUB( quot2, bot, p35 ), r, quot2 );	// quot2 + (p35 - quot2*bot)*r;		// second fraction wd.
		temp = 2.0*quot*quot2;																// reduced arg is (quot,quot2)
		argsq = __FMADD( quot, quot, temp );			// quot*quot + temp;				// compute square of argument
		argsq2 = __FMADD( quot2, quot2, __FMSUB( quot, quot, argsq) + temp ); // quot*quot - argsq + temp + quot2*quot2;
		sumt = __FMADD( argsq, atancoeff[7], atancoeff[6] );// atancoeff[6] + argsq*atancoeff[7];
		sum = __FMADD( argsq, sumt, atancoeff[5] );		// atancoeff[5] + argsq*sumt;
		suml = __FMADD( argsq, sumt, atancoeff[5] - sum ); // atancoeff[5] - sum + argsq*sumt;
		
		for (i = 4; i >= 1; i--) {	 	
			temps = __FMADD( sum, argsq, atancoeff[i] );	// atancoeff[i] + sum*argsq;
			/* suml = atancoeff[i] - temps + sum*argsq + (sum*argsq2 + suml*argsq); */
			suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
			sum = temps;
		}
		
		prodlow = __FMADD( suml, argsq, sum*argsq2 );		// suml*argsq + sum*argsq2;			// mult. by arg^2
		prod = __FMADD( sum, argsq, prodlow );				// sum*argsq + prodlow;
		prodlow = __FMSUB( sum, argsq, prod ) + prodlow;		// sum*argsq - prod + prodlow;
		temp2 = __FMADD( prodlow, quot, prod*quot2 );		// prodlow*quot + prod*quot2;		// mult. by arg
		temp = __FMADD( prod, quot, temp2 );				// prod*quot + temp2;
		temp2 = __FMSUB( prod, quot, temp ) + temp2;		// prod*quot - temp + temp2;
		sum = temp*kRecip;																		// sum of series ---
		remainder = __FNMSUB( sum, atancoeff[0], temp );		// (temp - sum*atancoeff[0]);
		partial = __FADD( remainder, temp2 );
		residual = remainder - partial + temp2;
		suml = __FMADD( partial, kRecip, (residual*kRecip) ); // partial*kRecip + (residual*kRecip);
		res = __FADD( quot, sum );																// except for argument
		reslow = (quot - res + sum);															// exact
		resmid = __FADD( quot2, suml );
		restiny = quot2 - resmid + suml;
		p = __FADD( reslow, resmid );															// sum of middle terms
		q = reslow - p + resmid;																// possible residual
		reshi = __FADD( res, p );
		resbot = (res - reshi + p) + (q + restiny);
		result = __FADD( atanPtr[(int32_t)index].Two, reshi );
		reslow = __FADD( atanPtr[(int32_t)index].Three, resbot );
		resmid = atanPtr[(int32_t)index].Two - result + reshi;
		restiny = resbot - reslow + atanPtr[(int32_t)index].Three;
		restemp = __FADD( resmid, reslow );
		reshi = __FADD( result, restemp );
		
		if (fabs(resmid) > fabs(reslow)) 
			restiny = resmid - restemp + reslow + restiny;
		else
			restiny = reslow - restemp + resmid + restiny;
		
		resbot = result - reshi + restemp + restiny;
		
		if (xHead < 0.0) {
			reshi = -reshi;
			resbot = -resbot;
			result = -result;
			restemp = -restemp;
			restiny = -restiny;
		}
		
		ddx.d[0] = __FADD( reshi, resbot );  
		ddx.d[1] = (reshi - ddx.d[0]) + resbot;		
		
                FESETENVD(fenv);
		return (ddx.f);								// Normal arguments done
	}
}

//  End of atanl, start of atan2l:

//*******************************************************************************
//     _QATAN2 expects two quad words as arguments.  It returns the arctan      *
//     of the quotient of the arguments, but differs from _QATAN in that        *
//     the result can be in the interval [-Pi,Pi], whereas _QATAN always        *
//     returns a result in the range [-Pi/2, Pi/2].  After argument reduction   *
//     and the testing of numerous special cases, _QATAN is used for the bulk   *
//     of the work.                                                             *
//                                                                              *
//       atan2(NaN, anything)--------------------------------->   NaN           *
//       atan2(anything, NaN)--------------------------------->   NaN           *
//       atan2(+-0, +(anything but NaN))----------------------> +-0             *
//       atan2(+-0, -(anything but NaN))----------------------> +-Pi            *
//       atan2(+-(anything but 0 and NaN), 0 )----------------> +-Pi/2          *
//       atan2(+-(anything but Infinity and NaN), +Infinity)--> +-0             *
//       atan2(+-(anything but Infinity and NaN), -Infinity)--> +-Pi            *
//       atan2(+-Infinity, +Infinity)-------------------------> +-Pi/4          *
//       atan2(+-Infinity, -Infinity)-------------------------> +-3Pi/4         *
//       atan2(+-Infinity, (anything but 0, NaN, Infinity))---> +-Pi/2          *
//                                                                              *
//*******************************************************************************

// Floating-point constants
static const double kTiny = 1.805194375864829576e-276;		// 0x1.0p-916 = 0x06B00000   
// kTiny: number that still allows a 107 bit quad to consist of two normal numbers

static const double kRescale = 2.923003274661805836e+48;	// 0x1.0p+161   
// Needed to rescale if the denominator is "tiny"

static const double kPi = 3.141592653589793116;				// 0x1.921FB54442D18p1   
static const double kPiTail1  = 1.224646799147353207e-16;	// 0x1.1A62633145C07p-53  
static const double kPiTail1m  = 1.224646799147352961e-16;	// 0x1.1A62633145C06p-53  
static const double kPiTail2 = -2.994769809718339666e-33;	// 0x1.F1976B7ED8FBCp-109  

static const double k3PiBy4 = 2.356194490192344837;			// 0x1.2D97C7F3321D2p1   
static const double k3PiBy4Tail1 = 9.184850993605148438e-17;// 0x1.A79394C9E8A0Ap-54  


long double atan2l( long double y, long double x )
{
	double p1, p2, p3, p2a, p4, p5, p6, p7, p8, p29, p30, p33, p35;
	double p31, p32, p34, p36, p37, p38;
	double quot, quot2;
	double rarg, result, restemp;
	double top, top2, topextra, top3;
	double extra;
	double bot, bot2;
	double r, rarg2, rarg2a;
	double temps;
	double remainder, residual;
	double partial;  
	double arg, argl, argsq, argsq2;
	double temp,temp2; 
	double sum, sumt, suml;
	double absx;
	double xHead, xTail;
	double res, reslow, resmid, restiny, reshi, resbot;
	double p, q; 
	double prod, prodlow;	
	double xDen, xDenTail, yNum, yNumTail;
	double frac3, frac2, frac, fract;
	double pin, pin2, pin3;
	double ressmall, restop, resint;
	Double DeN ;
	double fpenv;
	DblDbl ddx, ddy;
	
	int i, rflag;
	
	uint32_t index;
	uint32_t signx, signy, expx, expy, sigx, sigy;
	
	struct atantblEntry *atanPtr = (struct atantblEntry *)atanTable - 6;

        FEGETENVD(fpenv);
        FESETENVD(0.0);
	
	ddx.f = x;				// denominator
	ddy.f = y;				// numerator
	
	xDen     = ddx.d[0];	// dden
	xDenTail = ddx.d[1];	// ddenl    
	yNum     = ddy.d[0];	// dnum
	yNumTail = ddy.d[1];	// dnuml
	
	expy  = (ddy.i[0] >> 20) & 0x7ffu;
	expx  = (ddx.i[0] >> 20) & 0x7ffu;		// get biased exponents
	signy = (ddy.i[0] >> 31) & 1u;
	signx = (ddx.i[0] >> 31) & 1u;			// get signs
	sigy  = (ddy.i[0] & 0xfffffu) | ddy.i[1];
	sigx  = (ddx.i[0] & 0xfffffu) | ddx.i[1];
	
	if (expy == 0x7ffu) {					// y is infinite or NaN
		if (sigy != 0) {
                        FESETENVD(fpenv);
			return (y);						// NaN y is returned
		}
		else if (expx == 0x7ffu) {			// infinite y
			if (sigx != 0) {
                                FESETENVD(fpenv);
				return (x);					// NaN x is returned
			}
			else if (signx == 0) {			// infinite y and x  denominator positive
				if (signy) {				// numerator negative
					ddx.d[0] =  -0.5*kPiBy2 ; 	
					ddx.d[1] =  -0.5*kPiBy2Tail1;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
				else {						// numerator positive
					ddx.d[0] =  0.5*kPiBy2 ; 	
					ddx.d[1] =  0.5*kPiBy2Tail1;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
			}
			else {							// denominator negative
				if (signy) { 				// numerator negative
					ddx.d[0] =  -k3PiBy4 ; 	
					ddx.d[1] =  -k3PiBy4Tail1;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
				else {						// numerator positive
					ddx.d[0] =  k3PiBy4 ; 	
					ddx.d[1] =  k3PiBy4Tail1;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
			}
		}
		else {								// y is Infinite , x is finite
			if (signy == 0.0) {     
				if (signx == 0.0) {			// first quadrant 
					ddx.d[0] = kPiBy2 ; 	
					ddx.d[1] = kPiBy2Tail1m;
                                        FESETENVD(fpenv);
					return (ddx.f);
				} 
				else {						// second quadrant 
					ddx.d[0] = kPiBy2 ; 	
					ddx.d[1] = kPiBy2Tail1;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
			} 
			else  {							// yNum < 0.0
				if (signx == 0.0) {			// fourth quadrant
					ddx.d[0] = -kPiBy2 ; 	
					ddx.d[1] = -kPiBy2Tail1m;
                                        FESETENVD(fpenv);
					return (ddx.f);
				} 
				
				else  {						// third quadrant
					ddx.d[0] = -kPiBy2 ; 	
					ddx.d[1] = -kPiBy2Tail1;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
			}
		}
	}										// [y is infinite or NaN]
	if (expx == 0x7ffu) {					// x is infinite or NaN
		if (sigx != 0)	{					//   NaN x is returned
                        FESETENVD(fpenv);
			return (x);
		}
		else if (signx) {					// x is -INF	
			if (signy == 0) {
				ddx.d[0] = 2.0*kPiBy2 ; 	
				ddx.d[1] = 2.0*kPiBy2Tail1m;
                                FESETENVD(fpenv);
				return (ddx.f);
			}
			else {
				ddx.d[0] = -2.0*kPiBy2 ; 	
				ddx.d[1] = -2.0*kPiBy2Tail1m;
                                FESETENVD(fpenv);
				return (ddx.f);
			}
		}
		else								// x is +INF
		if (signy == 0) {
			ddx.d[0] = 0.0 ; 	
			ddx.d[1] = 0.0;
                        FESETENVD(fpenv);
			return (ddx.f);
		}
		else {
			ddx.d[0] = -0.0 ; 	
			ddx.d[1] = -0.0;
                        FESETENVD(fpenv);
			return (ddx.f);
		}
	}										// x is infinite or NaN
	if ((sigx | expx) == 0u) {				// x is zero
		if ((sigy | expy) == 0u) {			// y is zero
			if (signx) {
				if (signy) {
					ddx.d[0] = -2.0*kPiBy2 ; 	
					ddx.d[1] = -2.0*kPiBy2Tail1m;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
				else {
					ddx.d[0] = 2.0*kPiBy2 ; 	
					ddx.d[1] = 2.0*kPiBy2Tail1m;
                                        FESETENVD(fpenv);
					return (ddx.f);
				}
			}
			else {
                                FESETENVD(fpenv);
				return (y);
			}
		}
		else if (signy)	{					// y < 0.0
			if (signx) {
				ddx.d[0] = -kPiBy2 ; 	
				ddx.d[1] = -kPiBy2Tail1;
                                FESETENVD(fpenv);
				return (ddx.f);
			}
			else {
				ddx.d[0] = -kPiBy2 ; 	
				ddx.d[1] = -kPiBy2Tail1m;
                                FESETENVD(fpenv);
				return (ddx.f);
			}
		}
		else {								// y > 0.0
			if (signx) {
				ddx.d[0] = kPiBy2 ; 	
				ddx.d[1] = kPiBy2Tail1;	
                                FESETENVD(fpenv);
				return (ddx.f);
			}
			else {
				ddx.d[0] = kPiBy2 ; 	
				ddx.d[1] = kPiBy2Tail1m;
                                FESETENVD(fpenv);
				return (ddx.f);
			}
		}
	}										// [x is zero]
	if ((sigy | expy) == 0u) {				// [zero y, nonzero x]
		if (signx == 0.0) {					// denominator positive		
			ddx.d[0] = yNum;				// correctly signed zero result
			ddx.d[1] = 0.0;
                        FESETENVD(fpenv);
			return (ddx.f);
		} 
		else {								// negative denominator
			if (signy == 0) {				// numerator is 0, 
				ddx.d[0] = kPi ; 	
				ddx.d[1] = kPiTail1m;
                                FESETENVD(fpenv);
				return (ddx.f);
			} 
			else {							// numerator is -0  
				ddx.d[0] = -kPi ; 	
				ddx.d[1] = -kPiTail1m;
                                FESETENVD(fpenv);
				return (ddx.f);
			}
		}
	}										// [zero y, nonzero x]
	// END OF SPECIAL CASES
	
	// The code below is for finite y/x values, we first calculate the ratio y/x
	// for tiny values of either x and y, and for other finite values with extra precision
	// producing xHead, xTail and frac3. These quantities and the sign of the numerator and 
	// the denomitor are used then to produce ArcTan (y,x)
	
	//  Beginning of tiny numerator or denominator
	
	if((((ddy.i[0]) & 0x7fffffff) <  0x06B00000)|| (((ddx.i[0]) & 0x7fffffff) < 0x06B00000)) {
		xDen = xDen*kRescale;				// For a tiny arguments,
		xDenTail = xDenTail*kRescale;		// rescale all the input
		yNum = yNum*kRescale;				// so that the full length
		yNumTail = yNumTail*kRescale;		// division can procede.
		r = 1.0/xDen;
		if (fabs(r) < kTiny) {			// rescaling has made numerator
			if (xDen > 0) {					// too large.  Result is
				xHead = 0.0;				// zero for positive denominator
				xTail = 0.0;
				frac3 = 0.0;
			} 
			else {							// for neg. denominator,
				xHead = yNum;				// send numerator to qatan
				xTail = 0.0;				// result will be subtracted
				frac3 = 0.0;				// from properly signed Pi later
			}
		}
	}
	//  end of tiny numerator or denominator
	
	r = 1.0/xDen;							// start computation of num/den
	// normal non-special values of x and y
	fract = __FMUL( yNum, r );
	
	// A detailed division (yHead,yTail)/(xHead,xTail) follows
	frac = __FMADD( r, __FNMSUB( xDen, fract, yNum ), fract );	// fract + (yNum-xDen*fract)*r;		// refine h.o. quotient
	p29 = __FNMSUB( xDen, frac, yNum ); // yNum - xDen*frac;
	p30 = __FMUL( frac, xDenTail );					// negative term
	p31 = __FNMSUB( frac, xDenTail, p30 ); // p30 - frac*xDenTail;
	p32 = p31;
	p33 = __FSUB( p29, p30 );
	p35 = __FADD( p33, yNumTail );
	
	if (fabs(fract) == kInfinityu.f) {
		xHead = fract;
		xTail = 0.0;
		frac3 = 0.0;
	}
	else if (fabs(fract)< kInfinityu.f) {
		frac2 = __FMUL( p35, r );
		frac2 = __FMADD( r, __FNMSUB( frac2, xDen, p35 ), frac2 ); // frac2 + (p35-frac2*xDen)*r;	// mid. quot. word
		if (fabs(p29) > fabs(p30))
			p34 = p29 - p33 - p30 + p32;
		else 
			p34 = p29 - (p33+p30)+p32;
		p36 = __FNMSUB( frac2, xDen, p35 ) + __FNMSUB( frac2, xDenTail, p34 ); // p35 - frac2*xDen + p34 - frac2*xDenTail;
		if (fabs(p33) > fabs(yNumTail))
			p37 = p33 - p35 + yNumTail;
		else  
			p37 = yNumTail - p35 + p33;
		p38 = p37 + p36;
		frac3 = p38*r/__FMADD( frac, frac, 1.0 ); // p38*r/(1.0+frac*frac);		// Fraction done
		xHead = frac;						// Use ATAN routine to do
		xTail = frac2;						// most of the work
	}
	// Below we use the just evaluated values of xHead and xTail
	// to first evaluate the quantity "ressmall" which add an extra precision
	// contribution to frac3. Then, xHead, xTail, frac3, ressmall and 
	// the sign of the numerator (signy) and denominator (signx)
	// are used to evaluate ArcTan2(y,x)
	
	// finite result
	
	//***************************************************************************
	//      This routine breaks the arguments into five ranges:                 *
	//      Range 1:          0 <= |x| <= .0859375                              *
	//      Range 2:   .0859375 <= |x| <= 1.015625                              *
	//      Range 3:   1.015625 <= |x| < 13                                     *
	//      Range 4:         13 <= |x| < 2^107                                  *
	//      Range 5:      2^107 <= |x| <= Infinity                              *
	//                                                                          *
	//         >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>         *
	//                                                                          *
	// In Range 1, a 15th degree polynomial in x^2 is used to evaluate QATAN.   *
	//                                                                          *
	// In Range 2  and,                                                         *
	// In Range 3  the accurate table method is used to range-reduce the        *
	//             argument to a small value to be used with an interpolating   *
	//             polynomial for QATAN.  In Ranges 2 and 3, the resultant      *
	//             argument is less than 2^-7, so an 8th degree polynomial can  *
	//             be used.  Range 3 requires additional work to find the       *
	//             correct index into the accurate table.                       *
	//                                                                          *
	// In Range 4, QATAN(x) = Pi/2 - QATAN(1/x).  After taking the reciprocal   *
	//             of the argument, the quotient is less than 1/13, so that the *
	//             15th degree polynomial used for Range 1 is used to compute   *
	//             QATAN(1/x)                                                   *
	//                                                                          *
	// In Range 5, the correctly rounded value differs from Pi/2 by less than   *
	//             1/2 ulp.  It is, however, incorrect to let the value become  *
	//             correctly rounded, as that correctly rounded value exceeds   *
	//             PI/2 and is in the wrong quadrant.  Applying the inverse     *
	//             function, QTAN, to the correctly rounded value of PI/2 would *
	//             produce a result of the opposite sign to the argument        *
	//             originally given to ATAN.                                    *
	//***************************************************************************
	
	
	//  THIS CODE ONLY TREATS RANGES 1 TROUGH 5. WE ASSUME THAT xHead and xTail
	//	are such that we are ALWAYS in Range 1-5.
	
	absx = fabs(xHead);
	
	// Range 5
	
	if (absx >= kLimit4) {					// Filter out unusual arguments
		reshi = kPiBy2;						// if the second argument is
		if (signx == 0.0) {					// positive, the result is in
			resbot = kPiBy2Tail1m;			// first or fourth quadrant
			ressmall = 0;					// set up appropriate value of
			if (xHead < 0.0) {				// pi/2, truncated
				reshi = -reshi;
				resbot = -resbot;
			}
		}
		else {								// Otherwise, result is in the
			reshi = kPiBy2;					// second or third quadrant.
			resbot = kPiBy2Tail1;			// set up pi/2, correctly rounded
			ressmall = kPiBy2Tail2;
			if (xHead < 0.0) {
				reshi = -reshi;
				resbot = -resbot;
				ressmall = -ressmall;
			}
		}
	}

	// Ranges 1 and 4 (DO NOT USE TABLE)
	
	rflag = 0;	
	
	// Range 4
	
	if ((absx >= kLimit3) && (absx < kLimit4) ) {	// Range 4
		rarg = 1.0/xHead;                         
		rflag = 1;									// indicate large argument
		p1 = __FNMSUB( xHead, rarg, 1.0 );			// 1.0 - xHead*rarg;						// time to complete hi precision
		p2 = xTail*rarg;							// reciprocal of argument, to
		p3 = p1 - p2;								// use as reduced argument.
		rarg2a = p3*rarg;
		rarg2 = rarg2a + __FNMSUB( xHead, rarg2a, p3 )*rarg; // rarg2a + (p3 - xHead*rarg2a)*rarg;
		xHead = rarg;
		xTail = rarg2;
	}
	
	// Ranges 1 and 4 only
	
	if( (absx < kLimit1) || ((absx >= kLimit3) && (absx < kLimit4) )) {	//Code for Range 1 and Range 4
		temp = 2.0*xHead*xTail;						// direct evaluation of series
		argsq = __FMADD( xHead, xHead, temp );		// xHead*xHead + temp;					// compute square of argument
		argsq2 = __FMADD( xTail, xTail, __FMSUB( xHead, xHead, argsq ) + temp ); // xHead*xHead - argsq + temp + xTail*xTail;
		sum = atancoeff[14];
		sum = __FMADD( argsq, sum, atancoeff[13] );			// atancoeff[13] + argsq*sum;
		sum = __FMADD( argsq, sum, atancoeff[12] );			// atancoeff[12] + argsq*sum;
		sum = __FMADD( argsq, sum, atancoeff[11] );			// atancoeff[11] + argsq*sum;
		sumt = __FMADD( argsq, sum, atancoeff[10] );		// atancoeff[10] + argsq*sum;
		sum = __FMADD( argsq, sum, atancoeff[9] );			// atancoeff[9] + argsq*sumt;
		suml = __FMADD( argsq, sumt, atancoeff[9] - sum );	// atancoeff[9] - sum + argsq*sumt;
		
		for (i = 8; i >= 1; i--) {	 		   
			temps = __FMADD( argsq, sum, atancoeff[i] );	// atancoeff[i] + sum*argsq;
			/* suml = atancoeff[i] - temps+sum*argsq + (sum*argsq2 + suml*argsq); */
			suml = __FMADD( sum, argsq, atancoeff[i] -  temps ) + __FMADD( sum, argsq2, suml*argsq );
			sum=temps;
		}
		
		prodlow = __FMADD( suml, argsq, sum*argsq2 );		// suml*argsq + sum*argsq2;			// multiply by arg^2
		prod = __FMADD( sum, argsq, prodlow );				// sum*argsq + prodlow;
		prodlow = __FMSUB( sum, argsq, prod ) + prodlow;	// sum*argsq - prod + prodlow;
		temp2 = __FMADD( prodlow, xHead, prod*xTail );		// prodlow*xHead + prod*xTail;		// multiply by arg
		temp = __FMADD( prod, xHead, temp2 );				// prod*xHead + temp2;
		temp2 = __FMSUB( prod, xHead, temp ) + temp2;		// prod*xHead - temp + temp2;
		sum = temp*kRecip;																		// sum of series ---
		remainder = __FNMSUB( sum, atancoeff[0], temp );	// (temp - sum*atancoeff[0]);
		partial = __FADD( remainder, temp2 );
		residual = remainder - partial + temp2;
		suml = __FMADD( partial, kRecip, (residual*kRecip) );// partial*kRecip + (residual*kRecip);
		res = __FADD( xHead, sum );																// except for argument
		reslow = (xHead - res + sum);															// exact
		resmid = __FADD( xTail, suml );
		restiny = xTail - resmid + suml;
		p = __FADD( reslow, resmid );															// sum of middle terms
		q = reslow - p + resmid;																// possible residual
		reshi = res + p;
		resbot = (res - reshi + p) + (q + restiny);
		if (rflag == 1) {							// large magnitude argument
			if (xHead > 0) {						//  fetch pi/2 with proper sign
				p1 = kPiBy2;
				p2 = kPiBy2Tail1;
				p2a = kPiBy2Tail2;
			}
			else {
				p1 = -kPiBy2;
				p2 = -kPiBy2Tail1;
				p2a = -kPiBy2Tail2;
			}
			p3 = __FSUB( p1, reshi );				// subtract result from pi/2
			p4 = p1 - p3 - reshi;
			p5 = __FSUB( p2, resbot );
			p6 = p2 - p5 - resbot + p2a;
			p7 = __FADD( p4, p5 );
			reshi = p3 + p7;
			if (fabs(p4) > fabs(p5)) 
				p8 = p4 - p7 + p5;
			else
				p8 = p5 - p7 + p4;
			resbot = __FADD( p3 - reshi + p7, (p6+p8-restiny) );
			ressmall = (p3 - reshi + p7) - resbot + (p6 + p8 - restiny);
		}
		else   //if rflag = 0   i.e., Range 1
			ressmall = (res - reshi + p) - resbot + (q + restiny);
	}
	
	// Ranges 2 and 3  (USE TABLE)
	
	if( (absx >= kLimit1) && (absx < kLimit3) ) {	// Code for Range 2 and Range 3
		DeN.f = __FMADD( 64.0, absx, kBig );		// 64.0 *absx+kBig;
		if (absx > kLimit2) {						// for large args, use recip.  Range 3
			FESETENVD(kBig+1.0);						// use trunc mode momentarily
			DeN.f = absx + kBig;					// determine range of argument
			FESETENVD(kBig+1.0);						// use trunc mode momentarily
			DeN.f = absx + kBig;					// determine range of argument
			FESETENVD(kBig);							// return to round-to-nearest
			index = DeN.i[1];
			i = index - 1;
			DeN.f = __FMADD( factor[i], absx, kBig+adjust[i] ); // factor[i]*absx+kBig+adjust[i];
		}
		index = DeN.i[1];
		if (xHead < 0) {							// use abs. value
			arg = -xHead;
			argl = -xTail;
		}
		else {
			arg = xHead;
			argl = xTail;
		}
		top = __FSUB( arg, atanPtr[(int32_t)index].One );
		top2 = __FADD( top, argl );
		topextra = arg - top - atanPtr[(int32_t)index].One;
		top3 = top - top2 + argl + topextra;		// reduced numerator
		bot = __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 );	 // 1.0 + arg*atanPtr[(int32_t)index].One;
		/* bot2 = 1.0 - bot + arg*atanPtr[(int32_t)index].One + argl*atanPtr[(int32_t)index].One; */ // denominator
		bot2 = __FMADD( argl, atanPtr[(int32_t)index].One, __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 - bot ) );
		r = 1.0/bot;
		quot = r*top2;
		quot = __FMADD( __FNMSUB( bot, quot, top2 ) , r, quot ); // quot + (top2-bot*quot)*r;
		
		p29 = __FNMSUB( bot, quot, top2 );						//top2 - quot*bot;
		p30 = quot*bot2;							// negative term
		p33 = p29-p30;
		p35 = p33 + top3;
		quot2 = p35*r;
		quot2 = __FMADD( __FNMSUB( quot2, bot, p35 ), r, quot2 ); // quot2 + (p35-quot2*bot)*r;			// second fraction wd.
		temp = 2.0*quot*quot2;						// reduced arg is (quot,quot2)
		argsq = __FMADD( quot, quot, temp );			// quot*quot + temp;				// compute square of argument
		argsq2 = __FMADD( quot2, quot2, __FMSUB( quot, quot, argsq) + temp ); // quot*quot - argsq + temp + quot2*quot2;
		sumt = __FMADD( argsq, atancoeff[7], atancoeff[6] );// atancoeff[6] + argsq*atancoeff[7];
		sum = __FMADD( argsq, sumt, atancoeff[5] );		// atancoeff[5] + argsq*sumt;
		suml = __FMADD( argsq, sumt, atancoeff[5] - sum ); // atancoeff[5] - sum + argsq*sumt;
		
		for (i = 4; i >= 1; i--) {	 	
			temps = __FMADD( sum, argsq, atancoeff[i] );	// atancoeff[i] + sum*argsq;
			/* suml = atancoeff[i] - temps + sum*argsq + (sum*argsq2 + suml*argsq); */
			suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
			sum = temps;
		}
		
		prodlow = __FMADD( suml, argsq, sum*argsq2 );		// suml*argsq + sum*argsq2;			// mult. by arg^2
		prod = __FMADD( sum, argsq, prodlow );				// sum*argsq + prodlow;
		prodlow = __FMSUB( sum, argsq, prod ) + prodlow;		// sum*argsq - prod + prodlow;
		temp2 = __FMADD( prodlow, quot, prod*quot2 );		// prodlow*quot + prod*quot2;		// mult. by arg
		temp = __FMADD( prod, quot, temp2 );				// prod*quot + temp2;
		temp2 = __FMSUB( prod, quot, temp ) + temp2;		// prod*quot - temp + temp2;
		sum = temp*kRecip;																		// sum of series ---
		remainder = __FNMSUB( sum, atancoeff[0], temp );		// (temp - sum*atancoeff[0]);
		partial = __FADD( remainder, temp2 );
		residual = remainder - partial + temp2;
		suml = __FMADD( partial, kRecip, (residual*kRecip) ); // partial*kRecip + (residual*kRecip);
		res = __FADD( quot, sum );																// except for argument
		reslow = (quot - res + sum);															// exact
		resmid = __FADD( quot2, suml );
		restiny = quot2 - resmid + suml;
		p = __FADD( reslow, resmid );															// sum of middle terms
		q = reslow - p + resmid;																// possible residual
		reshi = __FADD( res, p );
		resbot = (res - reshi + p) + (q + restiny);
		result = __FADD( atanPtr[(int32_t)index].Two, reshi );
		reslow = __FADD( atanPtr[(int32_t)index].Three, resbot );
		resmid = atanPtr[(int32_t)index].Two - result + reshi;
		restiny = resbot - reslow + atanPtr[(int32_t)index].Three;
		restemp = __FADD( resmid, reslow );
		reshi = __FADD( result, restemp );
		
		if (fabs(resmid) > fabs(reslow)) 
			restiny = resmid - restemp + reslow + restiny;
		else
			restiny = reslow - restemp + resmid + restiny;
		resbot = result - reshi + restemp + restiny;
		if (xHead < 0.0) {
			reshi = -reshi;
			resbot = -resbot;
			result = -result;
			restemp = -restemp;
			restiny = -restiny;
		}
		ressmall = (result - reshi + restemp) - resbot + restiny;
	}
	extra = ressmall + frac3;
	if (signx == 0) {								// For Positive denominator
		ddx.d[0] = __FADD( reshi, (resbot + extra) );  
		ddx.d[1] = (reshi - ddx.d[0]) + (resbot + extra);					  
	}
	else {
		if (signy == 0) {							// For Negative denominator
			pin = kPi;								// return 2nd|3rd quadrand angle,
			pin2 = kPiTail1;						// depending on sign of numerator
			pin3 = kPiTail2;
		}
		else {
			pin = -kPi;								// -->  Pi+atan, for dnum < 0,  
			pin2 = -kPiTail1;						//      atan-Pi, for dnum > 0
			pin3 = -kPiTail2;
		}
		restop = __FADD( pin, reshi );
		resmid = pin - restop + reshi;
		reslow = __FADD( pin2, resbot );
		restiny = pin2 - reslow + resbot + extra;
		resint = __FADD( resmid, reslow );
		reshi = __FADD( restop, resint );
		resbot = restop - reshi+resint;
		if (fabs(resmid) > fabs(reslow)) 
			resbot = resbot + (resmid-resint+reslow+restiny);
		else 
			resbot = resbot + (reslow-resint+resmid+restiny);
		ddx.d[0] = __FADD( reshi, resbot );  
		ddx.d[1] = (reshi - ddx.d[0]) + resbot;		
	}

        FESETENVD(fpenv);
	return (ddx.f);
}
#endif