overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22//
23// ArcTanDD.c
24//
25// Double-double Function Library
26// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
27//
28// long double atanl( long double x );
29// long double atan2l( long double y, long double x );
30//
31// Change History:
32// 8/31/96 - PAF - Changed if() to else (see Note 1)
33//
34
35#include "math.h"
36#include "fp_private.h"
37#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
38#include "DD.h"
39
40static const double kBig = 6.755399441055744e+15; // 0x1.8p52
41
42static const double kPiBy2 = 1.570796326794896619231322; // 0x1.921FB54442D18p0
43static const double kPiBy2Tail1 = 6.123233995736766036e-17; // 0x1.1A62633145C07p-54
44static const double kPiBy2Tail1m = 6.123233995736764803e-17; // 0x1.1A62633145C06p-54
45static const double kPiBy2Tail2 = -1.4973849048591698e-33; // 0x1.F1976B7ED8FBCp-110
46static const double kRecip = 6.869637070016756334588131e-12; // 0x1.e368613a3b430p-38
47static const double kLimit1 = 0.0859375; // 0x1.6p-4
48static const double kLimit2 = 1.015625;
49static const double kLimit3 = 13.0;
50static const double kLimit4 = 8.11296384146066817e+31; // 0x1.0p+106
51static const Float kInfinityu = {0x7f800000};
52
53//***************************************************************************
54// factor -> gives the number of table values between the consecutive *
55// integers 1 to 12 *
56//***************************************************************************
57
58static const float factor[12] = {
59 32.0f,16.0f,8.0f,4.0f,4.0f,2.0f,2.0f,1.0f,1.0f,1.0f,1.0f,1.0f
60};
61
62
63//***************************************************************************
64// Adjust -> index into accurate atan table, for integer arguments *
65// from 1 to 12. *
66//***************************************************************************
67
68static const float adjust[12] = {
69 32.0f,64.0f,88.0f,104.0f,104.0f,116.0f,116.0f,124.0f,124.0f,124.0f,124.0f,124.0f
70};
71
72
73//***************************************************************************
74// atancoeff: coefficients of atan approximating polynomial, *
75// multiplied by 145568097675.0 *
76//***************************************************************************
77
78static const double atancoeff[15] = {
79 145568097675.0, -48522699225.0,
80 29113619535.0, -20795442525.0, 16174233075.0,
81 -13233463425.0, 11197545975.0, -9704539845.0,
82 8562829275.0, -7661478825.0, 6931814175.0,
83 -6329047725.0, 5822723907.0, -5391411025.0,
84 5019589575.0
85};
86
87static const uint32_t atanTable[] __attribute__ ((aligned(8))) = {
88
89 //*******************************************************************
90 // Accurate table over range ( .093750, .984375). *
91 //*******************************************************************
92
93 0x3FB80000, 0x96AC4A5A, 0x3FB7EE18, 0xBB5F2BDC, 0x3B1F7151, 0xA60EEDFD,
94 0x3FBC0000, 0xAB2E5990, 0x3FBBE39F, 0x679755D8, 0x3B5AF7C1, 0xBBB2A139,
95 0x3FC00000, 0x7F748304, 0x3FBFD5BB, 0x95A94110, 0x3B5B0A0C, 0x811C880E,
96 0x3FC20000, 0x3DF2157F, 0x3FC1E1FB, 0x37C2C7E5, 0x3B617207, 0xB776AF4F,
97 0x3FC40000, 0x30754446, 0x3FC3D6EF, 0x1814018D, 0xBB63A092, 0x88503B16,
98 0x3FC60000, 0x7996ECD4, 0x3FC5C981, 0x94588C25, 0xBB5B6979, 0x558D8B64,
99 0x3FC80000, 0x76D31B22, 0x3FC7B97B, 0xBE985C07, 0x3B6708AD, 0xCC1C4519,
100 0x3FCA0000, 0xEADCE47A, 0x3FC9A6A9, 0xCAFAF9A8, 0x3B62971C, 0xA562CA82,
101 0x3FCC0000, 0x7E8C3F27, 0x3FCB90D7, 0xCB57348B, 0xBB493C17, 0xA875159C,
102 0x3FCE0000, 0x5B063A6A, 0x3FCD77D6, 0x3569311A, 0xBB5A33AB, 0xA192F16A,
103 0x3FD00000, 0x64EF22FF, 0x3FCF5B76, 0xB72AE095, 0x3B52E368, 0xB50909D6,
104 0x3FD10000, 0x87DE9AAA, 0x3FD09DC6, 0x16C29778, 0xBB55B0DC, 0xF3418D68,
105 0x3FD20000, 0x0AAD7C84, 0x3FD18BF5, 0xACF10E71, 0x3B41EAC1, 0x346154B4,
106 0x3FD30000, 0x1D92088F, 0x3FD27837, 0x3B84D2FF, 0x3B48BBE9, 0xDC9998E5,
107 0x3FD40000, 0xCC4E6285, 0x3FD36277, 0xF12910F6, 0xBB6585F3, 0xFA8B319B,
108 0x3FD50000, 0x766A7592, 0x3FD44AA4, 0xA1AA8D78, 0x3B6C6289, 0x13F656CF,
109 0x3FD60000, 0xCFA33544, 0x3FD530AE, 0x5303BB5D, 0x3B63F533, 0xC1117E20,
110 0x3FD70000, 0xC925EED6, 0x3FD61484, 0xB52DF321, 0xBB70059A, 0x66E93FA9,
111 0x3FD80000, 0x54AA10D9, 0x3FD6F619, 0x8C1ECA84, 0x3B60C60C, 0x08BED883,
112 0x3FD90000, 0x113AC5B7, 0x3FD7D560, 0x5A568BA4, 0xBB541817, 0xB5654BFC,
113 0x3FDA0000, 0x6D101F2E, 0x3FD8B24D, 0x96E71249, 0xBB55A721, 0x1379C827,
114 0x3FDB0000, 0xE246C017, 0x3FD98CD6, 0x05641F92, 0x3B6918FF, 0xCA4E1DF9,
115 0x3FDC0000, 0x39B58346, 0x3FDA64EE, 0xF43C341A, 0x3B6633D0, 0xEECC747A,
116 0x3FDD0000, 0x4C36F71C, 0x3FDB3A91, 0x5CE1B431, 0xBB61E754, 0x0B1A4FB4,
117 0x3FDE0000, 0x66CA921B, 0x3FDC0DB5, 0x1D94F195, 0x3B571DE1, 0xE6981964,
118 0x3FDF0000, 0x3AD3639E, 0x3FDCDE53, 0x72D1AC94, 0xBB736718, 0x370B3962,
119 0x3FE00000, 0x7815CCCD, 0x3FDDAC67, 0xC5849B77, 0x3B52FF7C, 0xC7D84489,
120 0x3FE08000, 0xF40DDDFC, 0x3FDE77ED, 0x00AEBFCD, 0x3B62E958, 0xE33FF707,
121 0x3FE10000, 0xA6C93DAD, 0x3FDF40DE, 0x0F7AA915, 0x3B77641D, 0x15CA0475,
122 0x3FE18000, 0x85EBF4D7, 0x3FE0039C, 0xDAD8BEF9, 0x3B722996, 0x8391BF09,
123 0x3FE20000, 0xC129DA76, 0x3FE0657F, 0x27977717, 0xBB761869, 0xC045B9BD,
124 0x3FE28000, 0x93F751BC, 0x3FE0C614, 0xCA41B0AE, 0xBB7A8A92, 0x3C961154,
125 0x3FE30000, 0x7C1EF6CE, 0x3FE1255D, 0xF7C0A55C, 0x3B591CBF, 0x478B1C38,
126 0x3FE38000, 0x169FA349, 0x3FE1835A, 0x993DD546, 0x3B6FBFAB, 0x698175A0,
127 0x3FE40000, 0x27D6B511, 0x3FE1E00B, 0xC884E585, 0xBB248C8D, 0x9BB39A1C,
128 0x3FE48000, 0x6BE89876, 0x3FE23B72, 0x2F4EF788, 0xBB815690, 0x7DF099CF,
129 0x3FE50000, 0x280F4ABE, 0x3FE2958E, 0x7530C25F, 0xBB72B7E0, 0xEE07EAE7,
130 0x3FE58000, 0x2FEDFD47, 0x3FE2EE62, 0xA50C9A07, 0xBB2B522D, 0xBFFA3A8C,
131 0x3FE60000, 0x875C63BA, 0x3FE345F0, 0x78B8C001, 0xBB408ED6, 0xD056E4C4,
132 0x3FE68000, 0x8A368DBA, 0x3FE39C39, 0x79511750, 0x3B77CB64, 0x0C777681,
133 0x3FE70000, 0xEE7CA285, 0x3FE3F140, 0x55DECB32, 0xBB5AA9DC, 0xBC68CF79,
134 0x3FE78000, 0xA4904A32, 0x3FE44506, 0xC661B510, 0x3B7EBCE3, 0x3EB6C0CD,
135 0x3FE80000, 0x4CA9F6E2, 0x3FE4978F, 0xD4373CAA, 0x3AFF9675, 0xAB73F345,
136 0x3FE88000, 0xD7B95553, 0x3FE4E8DE, 0xE3B770C7, 0xBB715BD3, 0xBFA6C0AD,
137 0x3FE90000, 0x849C7E04, 0x3FE538F5, 0xCDE27023, 0xBB6B45C2, 0x3E5A5CC8,
138 0x3FE98000, 0xC171C662, 0x3FE587D8, 0x95C38C2E, 0x3B7887A1, 0x345B0752,
139 0x3FEA0000, 0x4321104C, 0x3FE5D589, 0xAF86142F, 0xBB8438EB, 0x1AE94AC1,
140 0x3FEA8000, 0x838D8D26, 0x3FE6220D, 0x5F66C783, 0x3B7CA422, 0x8F96E72F,
141 0x3FEB0000, 0xE4DA897E, 0x3FE66D66, 0xBED2B224, 0xBB35FB4D, 0xFD76AD75,
142 0x3FEB8000, 0xD775765C, 0x3FE6B799, 0x0DF9D017, 0x3B7B5383, 0xE0D80102,
143 0x3FEC0000, 0xAA053CA7, 0x3FE700A8, 0x25C3BB63, 0xBB7D5234, 0x29157DFC,
144 0x3FEC8000, 0x773AAC9B, 0x3FE74897, 0xD237C658, 0x3B6ACFB2, 0xF4E46A15,
145 0x3FED0000, 0x3FA5ABCF, 0x3FE78F6B, 0xE04F6C27, 0xBB8680F4, 0x55D1C7DF,
146 0x3FED8000, 0xD31E46AA, 0x3FE7D528, 0x9AC0235F, 0xBB61A1B0, 0x1EABC199,
147 0x3FEE0000, 0xBB676A41, 0x3FE819D1, 0x1AD33A87, 0xBB7726EA, 0x735248A8,
148 0x3FEE8000, 0x76CB5B54, 0x3FE85D69, 0x95ABE3CC, 0x3B772D9C, 0x044A3CEB,
149 0x3FEF0000, 0x26506206, 0x3FE89FF6, 0x131BC92E, 0x3B860889, 0x51E5F051,
150 0x3FEF8000, 0xED69BF50, 0x3FE8E17B, 0x2230274E, 0xBB7FB7C5, 0x5CA2A11A,
151
152 //*******************************************************************
153 // Accurate table over range ( 1.000000, 1.968750). *
154 //*******************************************************************
155
156 0x3FF00000, 0x38A5805B, 0x3FE921FB, 0x8CE9AD0F, 0xBB728518, 0x56659FDA,
157 0x3FF08000, 0xE4CDB94C, 0x3FE9A001, 0x86F9991F, 0xBB7F4F21, 0xC242FD58,
158 0x3FF10000, 0x29DA60EA, 0x3FEA1A26, 0x1A19C30A, 0xBB87F5DE, 0xE73B7B3A,
159 0x3FF18000, 0x9BBB03C5, 0x3FEA908B, 0x882B14EE, 0x3B8A364F, 0xC319894C,
160 0x3FF20000, 0x42A338B5, 0x3FEB034F, 0x7337C8D1, 0xBB768B38, 0xD9BDAC43,
161 0x3FF28000, 0xED0428FC, 0x3FEB7292, 0x7FBAC8B6, 0x3B7B400A, 0x53FD9AAF,
162 0x3FF30000, 0x36946784, 0x3FEBDE71, 0x1A8E3A61, 0x3B6C2B37, 0x13D123E1,
163 0x3FF38000, 0x21F16647, 0x3FEC470A, 0xDA7DBAA7, 0xBB523353, 0xD7A5886A,
164 0x3FF40000, 0xF06227B4, 0x3FECAC7D, 0x132231D2, 0xBB8AEFE3, 0x0558B34B,
165 0x3FF48000, 0xF9AB1CF7, 0x3FED0EE2, 0xE23FB012, 0xBB4465DC, 0xD4A6C05D,
166 0x3FF50000, 0x1DAE2F90, 0x3FED6E57, 0xE51C7E16, 0xBB85A9B9, 0x83503D52,
167 0x3FF58000, 0xD679A51F, 0x3FEDCAF8, 0xC6A4C801, 0x3B8576BE, 0xCEAEFC33,
168 0x3FF60000, 0xCE367C5E, 0x3FEE24DD, 0xD375A188, 0xBB78F5B1, 0x1E4802A0,
169 0x3FF68000, 0xCD05D052, 0x3FEE7C20, 0xCBEB8A43, 0xBB8AD216, 0x43D5F156,
170 0x3FF70000, 0x98AA0697, 0x3FEED0D9, 0xE022B144, 0x3B267908, 0xB6C976B5,
171 0x3FF78000, 0xA3FB1493, 0x3FEF2320, 0xDB8F1029, 0xBB85B081, 0x83128CFA,
172 0x3FF80000, 0x681E0063, 0x3FEF730C, 0x12946C3F, 0xBB8AA69F, 0xBF359098,
173 0x3FF88000, 0x75E901B6, 0x3FEFC0B1, 0xB86DE15F, 0x3B8EAA0C, 0x4663839A,
174 0x3FF90000, 0x736A4590, 0x3FF00613, 0x4FBE5B7E, 0xBB4E3551, 0x42306F18,
175 0x3FF98000, 0xC5B189FC, 0x3FF02ABF, 0xA107BE50, 0x3B63FD1C, 0xCE1D6AEB,
176 0x3FFA0000, 0x41411D33, 0x3FF04E67, 0x3966892F, 0xBB86ED08, 0x00B1EF7B,
177 0x3FFA8000, 0x4DDFE3DA, 0x3FF07113, 0xDB772C28, 0x3B746CEB, 0x807EA687,
178 0x3FFB0000, 0xA7AB4BE2, 0x3FF092CE, 0x72AC060E, 0xBB72F0AF, 0x054AD94E,
179 0x3FFB8000, 0x35EE2E87, 0x3FF0B39F, 0x5C6DBFD7, 0xBB8A699D, 0x14CF8D89,
180 0x3FFC0000, 0x79493CCE, 0x3FF0D38F, 0x4A3683E2, 0xBB7C178E, 0x4E0FB167,
181 0x3FFC8000, 0x7C2C16CD, 0x3FF0F2A5, 0xF7C1C52E, 0xBB77BBEC, 0x0F5FF2D9,
182 0x3FFD0000, 0xF7925B8C, 0x3FF110EB, 0x3A456E0E, 0x3B86B850, 0xB94665FD,
183 0x3FFD8000, 0xD4EDAFAA, 0x3FF12E66, 0x2A98E0EA, 0x3B8FC04D, 0xF88FF47B,
184 0x3FFE0000, 0x1822F5F9, 0x3FF14B1D, 0xDB5166C5, 0x3B400FEE, 0xA21D9DDA,
185 0x3FFE8000, 0x5F375D4F, 0x3FF16719, 0x6EAC8C79, 0xBB45D7BA, 0x99C8508E,
186 0x3FFF0000, 0x983AF36D, 0x3FF1825F, 0x2745DBAE, 0xBB8849C2, 0x328CE199,
187 0x3FFF8000, 0xDF7044EE, 0x3FF19CF5, 0x48D90CBC, 0xBB876AC4, 0x23C76E06,
188
189 //******************************************************************
190 // Accurate table over range ( 2.000000, 2.937500). *
191 //******************************************************************
192
193 0x40000000, 0xC89A22AD, 0x3FF1B6E1, 0xE3296298, 0xBB83E57A, 0x4E33C68C,
194 0x40008000, 0xAE85BCDD, 0x3FF1E8D4, 0xB635433E, 0x3B8FDE86, 0xD38237A1,
195 0x40010000, 0x1E66B4AD, 0x3FF21862, 0xFF00ECE6, 0x3B74E47A, 0x940D38B0,
196 0x40018000, 0x257005B2, 0x3FF245B5, 0x07E38198, 0xBB90341D, 0xF7AEB5A2,
197 0x40020000, 0x55656B58, 0x3FF270EF, 0x71D13D73, 0x3B88C680, 0x13906BD2,
198 0x40028000, 0x466225D9, 0x3FF29A34, 0x0FA8F5E6, 0x3B91F854, 0xF983041D,
199 0x40030000, 0x717EC3D3, 0x3FF2C1A2, 0x640509A5, 0x3B8E1C75, 0xAF6F1BFF,
200 0x40038000, 0x75A7B542, 0x3FF2E757, 0x4A697F72, 0x3B92AB3B, 0x412F53FB,
201 0x40040000, 0x00512EF1, 0x3FF30B6D, 0x7980B2E2, 0x3B8AD702, 0x40329E0C,
202 0x40048000, 0x09ED9B8A, 0x3FF32DFE, 0x0460EC6D, 0x3B6B92BF, 0x54B284E9,
203 0x40050000, 0x1B9F04CA, 0x3FF34F1F, 0xC21A2D1B, 0xBB82C1A6, 0xD9CE49BA,
204 0x40058000, 0x697C1A2B, 0x3FF36EE8, 0x0C4A7DD1, 0x3B73DD4E, 0x5C021394,
205 0x40060000, 0xABB8A398, 0x3FF38D6A, 0x94FD5FEE, 0x3B373349, 0x3B4FC3A6,
206 0x40068000, 0x18374FED, 0x3FF3AAB9, 0x8BB1767A, 0x3B72834E, 0xCBC4CA13,
207 0x40070000, 0xBCD94FDA, 0x3FF3C6E6, 0x7976E8D9, 0xBB931322, 0x7228DCD0,
208 0x40078000, 0x7BA476B6, 0x3FF3E200, 0xC84E84B4, 0xBB87D21A, 0xC44FB391,
209
210 //******************************************************************
211 // Accurate table over range ( 3.000000, 3.875000). *
212 //******************************************************************
213
214 0x40080000, 0xD7180754, 0x3FF3FC17, 0x967F5249, 0xBB64B771, 0x5FA0957B,
215 0x40090000, 0x6DFDCD42, 0x3FF42D70, 0x558EAD29, 0x3B84B923, 0xDD410F01,
216 0x400A0000, 0xE7C427AE, 0x3FF45B54, 0xAB8A2C51, 0x3B91D694, 0xE38B7BFF,
217 0x400B0000, 0x10E3223E, 0x3FF4861B, 0x4FB5B5B4, 0xBB7CB027, 0x50BA194A,
218 0x400C0000, 0x8EA13F5E, 0x3FF4AE11, 0x11ECF83F, 0x3B8B6616, 0xD89A7B16,
219 0x400D0000, 0x06D105F7, 0x3FF4D378, 0xC2906964, 0x3B93659C, 0x0737481A,
220 0x400E0000, 0x727DD5BC, 0x3FF4F68D, 0xF99AE908, 0x3B8514AE, 0x4636F80B,
221 0x400F0000, 0x3FDFD1C5, 0x3FF51785, 0x020F4064, 0xBB65DB78, 0x8577ADC5,
222
223 //******************************************************************
224 // Accurate table over range ( 4.000000, 5.750000). *
225 //******************************************************************
226
227 0x40100000, 0xCF49E240, 0x3FF5368C, 0xC5E4B1C8, 0xBB5E5499, 0xA9C4BF51,
228 0x40110000, 0x94E90733, 0x3FF56F6F, 0x52E31047, 0x3B8F3B21, 0x724AC4D7,
229 0x40120000, 0x46578102, 0x3FF5A250, 0x5F4EF405, 0x3B8657A1, 0x1EBE66F0,
230 0x40130000, 0xCA3D6D0F, 0x3FF5D013, 0xE66FF8CA, 0x3B60B423, 0x26CAA3CA,
231 0x40140000, 0x32964FAE, 0x3FF5F973, 0x1CEDA34B, 0x3B635D66, 0xC3ED1C18,
232 0x40150000, 0xA6AA4915, 0x3FF61F06, 0xDE005085, 0x3B6E8C7F, 0x000B0CDD,
233 0x40160000, 0x8906EAE3, 0x3FF6414D, 0x5593660A, 0x3B95FADD, 0xB8BF7A98,
234 0x40170000, 0x831FA60C, 0x3FF660B0, 0x3BD94D47, 0x3B96491D, 0x5B9C83F6,
235
236 //******************************************************************
237 // Accurate table over range ( 6.000000, 7.500000). *
238 //******************************************************************
239
240 0x40180000, 0x8DCF9463, 0x3FF67D88, 0x73114F7D, 0x3B928119, 0xC1230051,
241 0x401A0000, 0xCBB9EEAD, 0x3FF6B0BA, 0xFB083C0D, 0x3B93280F, 0xAE0D69E3,
242 0x401C0000, 0xDEF559E1, 0x3FF6DCC5, 0x8D8B9598, 0x3B86722D, 0x1022C04A,
243 0x401E0000, 0x744BF1EC, 0x3FF7030D, 0x01605442, 0xBB736E8B, 0xF92A87CA,
244
245 //******************************************************************
246 // Accurate table over range ( 8.000000, 12.000000). *
247 //******************************************************************
248
249 0x40200000, 0x456D4528, 0x3FF7249F, 0xB324E4B7, 0x3B675690, 0xDD79F342,
250 0x40220000, 0x43E2466D, 0x3FF75CBA, 0xD9437CC1, 0x3B46FCCD, 0xEB45ACCD,
251 0x40240000, 0x18A681EA, 0x3FF789BD, 0x2E09D7EA, 0x3B6C324D, 0xB28DD9D9,
252 0x40260000, 0xF3FE273A, 0x3FF7AEA3, 0x9C1AAC21, 0x3B8D6F66, 0xFC2D1461,
253 0x40280000, 0x455A6241, 0x3FF7CD6F, 0x71992B07, 0x3B8B2FE8, 0xBAEBACBC,
254 0x40280000, 0x455A6241, 0x3FF7CD6F, 0x71992B07, 0x3B8B2FE8, 0xBAEBACBC
255};
256
257struct atantblEntry {
258 double One;
259 double Two;
260 double Three;
261} atantblEntry;
262
263long double atanl( long double x )
264{
265 double p1, p2, p3, p2a, p4, p5, p6, p7, p8, p29, p30, p33, p35;
266 double quot, quot2;
267 double rarg, result, restemp;
268 double top, top2, topextra, top3;
269 double bot, bot2;
270 double r, rarg2, rarg2a;
271 double temps;
272 double remainder, residual;
273 double partial;
274 double arg, argl, argsq, argsq2;
275 double temp,temp2;
276 double sum, sumt, suml;
277 double absx;
278 double xHead, xTail;
279 double res, reslow, resmid, restiny, reshi, resbot;
280 double p, q;
281 double prod, prodlow;
282 Double DeN ;
283 double fenv;
284 DblDbl ddx;
285 int i, rflag;
286 uint32_t index;
287
288 struct atantblEntry *atanPtr = (struct atantblEntry *)atanTable - 6;
289
290 FEGETENVD(fenv);
291 FESETENVD(0.0);
292
293 ddx.f = x;
294 xHead = ddx.d[0];
295 xTail = ddx.d[1];
296
297 // x is a NaN or zero
298
299 if ((xHead != xHead)||(xHead == 0.0)) { // x is 0 or a NaN
300 FESETENVD(fenv);
301 return x;
302 }
303
304 // x is not a NaN nor zero
305
306 absx = fabs(xHead);
307
308 //finite result
309
310 //***************************************************************************
311 // This routine breaks the arguments into five ranges: *
312 // Range 1: 0 <= |x| <= .0859375 *
313 // Range 2: .0859375 <= |x| <= 1.015625 *
314 // Range 3: 1.015625 <= |x| < 13 *
315 // Range 4: 13 <= |x| < 2^107 *
316 // Range 5: 2^107 <= |x| <= Infinity *
317 // *
318 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> *
319 // *
320 // In Range 1, a 15th degree polynomial in x^2 is used to evaluate QATAN. *
321 // *
322 // In Range 2 and, *
323 // In Range 3 the accurate table method is used to range-reduce the *
324 // argument to a small value to be used with an interpolating *
325 // polynomial for QATAN. In Ranges 2 and 3, the resultant *
326 // argument is less than 2^-7, so an 8th degree polynomial can *
327 // be used. Range 3 requires additional work to find the *
328 // correct index into the accurate table. *
329 // *
330 // In Range 4, QATAN(x) = Pi/2 - QATAN(1/x). After taking the reciprocal *
331 // of the argument, the quotient is less than 1/13, so that the *
332 // 15th degree polynomial used for Range 1 is used to compute *
333 // QATAN(1/x) *
334 // *
335 // In Range 5, the correctly rounded value differs from Pi/2 by less than *
336 // 1/2 ulp. It is, however, incorrect to let the value become *
337 // correctly rounded, as that correctly rounded value exceeds *
338 // PI/2 and is in the wrong quadrant. Applying the inverse *
339 // function, QTAN, to the correctly rounded value of PI/2 would *
340 // produce a result of the opposite sign to the argument *
341 // originally given to ATAN. *
342 //***************************************************************************
343
344 // Range 5
345
346 if (absx >= kLimit4) { // Filter out unusual arguments
347 ddx.d[0] = kPiBy2;
348 ddx.d[1] = kPiBy2Tail1m; // Intentional misround to
349 if (xHead < 0.0) {
350 ddx.d[0] = -ddx.d[0];
351 ddx.d[1] = -ddx.d[1];
352 }
353 FESETENVD(fenv);
354 return (ddx.f);
355 }
356
357 // Ranges 1 and 4 (DO NOT USE TABLE)
358
359 rflag = 0;
360
361 // Range 4
362
363 if (absx >= kLimit3) { // Range 4
364 rarg = 1.0/xHead;
365 rflag = 1; // indicate large argument
366 p1 = __FNMSUB( xHead, rarg, 1.0 ); // 1.0 - xHead*rarg; // time to complete hi precision
367 p2 = xTail*rarg; // reciprocal of argument, to
368 p3 = p1 - p2; // use as reduced argument.
369 rarg2a = p3*rarg;
370 rarg2 = rarg2a + __FNMSUB( xHead, rarg2a, p3 )*rarg; // rarg2a + (p3 - xHead*rarg2a)*rarg;
371 xHead = rarg;
372 xTail = rarg2;
373 }
374
375 // Ranges 1 and 4 only
376
377 if( (absx < kLimit1) || (absx >= kLimit3) ) { // Code for Range 1 and Range 4
378 temp = 2.0*xHead*xTail; // direct evaluation of series
379 argsq = __FMADD( xHead, xHead, temp ); // xHead*xHead + temp; // compute square of argument
380 argsq2 = __FMADD( xTail, xTail, __FMSUB( xHead, xHead, argsq ) + temp ); // xHead*xHead - argsq + temp + xTail*xTail;
381 sum = atancoeff[14];
382 sum = __FMADD( argsq, sum, atancoeff[13] ); // atancoeff[13] + argsq*sum;
383 sum = __FMADD( argsq, sum, atancoeff[12] ); // atancoeff[12] + argsq*sum;
384 sum = __FMADD( argsq, sum, atancoeff[11] ); // atancoeff[11] + argsq*sum;
385 sumt = __FMADD( argsq, sum, atancoeff[10] ); // atancoeff[10] + argsq*sum;
386 sum = __FMADD( argsq, sum, atancoeff[9] ); // atancoeff[9] + argsq*sumt;
387 suml = __FMADD( argsq, sumt, atancoeff[9] - sum ); // atancoeff[9] - sum + argsq*sumt;
388
389 for (i = 8; i >= 1; i--) {
390 temps = __FMADD( argsq, sum, atancoeff[i] ); // atancoeff[i] + sum*argsq;
391 /* suml = atancoeff[i] - temps+sum*argsq + (sum*argsq2 + suml*argsq); */
392 suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
393 sum=temps;
394 }
395
396 prodlow = __FMADD( suml, argsq, sum*argsq2 ); // suml*argsq + sum*argsq2; // multiply by arg^2
397 prod = __FMADD( sum, argsq, prodlow ); // sum*argsq + prodlow;
398 prodlow = __FMSUB( sum, argsq, prod ) + prodlow; // sum*argsq - prod + prodlow;
399 temp2 = __FMADD( prodlow, xHead, prod*xTail ); // prodlow*xHead + prod*xTail; // multiply by arg
400 temp = __FMADD( prod, xHead, temp2 ); // prod*xHead + temp2;
401 temp2 = __FMSUB( prod, xHead, temp ) + temp2; // prod*xHead - temp + temp2;
402 sum = temp*kRecip; // sum of series ---
403 remainder = __FNMSUB( sum, atancoeff[0], temp ); // (temp - sum*atancoeff[0]);
404 partial = __FADD( remainder, temp2 );
405 residual = remainder - partial + temp2;
406 suml = __FMADD( partial, kRecip, (residual*kRecip) );// partial*kRecip + (residual*kRecip);
407 res = __FADD( xHead, sum ); // except for argument
408 reslow = (xHead - res + sum); // exact
409 resmid = __FADD( xTail, suml );
410 restiny = xTail - resmid + suml;
411 p = __FADD( reslow, resmid ); // sum of middle terms
412 q = reslow - p + resmid; // possible residual
413 reshi = res + p;
414 resbot = (res - reshi + p) + (q + restiny);
415 if (rflag == 1) { // large magnitude argument
416 if (xHead > 0) { // fetch pi/2 with proper sign
417 p1 = kPiBy2;
418 p2 = kPiBy2Tail1;
419 p2a = kPiBy2Tail2;
420 }
421 else {
422 p1 = -kPiBy2;
423 p2 = -kPiBy2Tail1;
424 p2a = -kPiBy2Tail2;
425 }
426 p3 = __FSUB( p1, reshi ); // subtract result from pi/2
427 p4 = p1 - p3 - reshi;
428 p5 = __FSUB( p2, resbot );
429 p6 = p2 - p5 - resbot + p2a;
430 p7 = __FADD( p4, p5 );
431 reshi = __FADD( p3, p7 );
432 if (fabs(p4) > fabs(p5))
433 p8 = p4 - p7 + p5;
434 else
435 p8 = p5 - p7 + p4;
436
437 resbot = p3 - reshi + p7 + (p6 + p8 - restiny);
438 }
439 ddx.d[0] = __FADD( reshi, resbot );
440 ddx.d[1] = (reshi - ddx.d[0]) + resbot;
441
442 FESETENVD(fenv);
443 return (ddx.f);
444 }
445
446 // Ranges 2 and 3 (USE TABLE)
447
448 //if( (absx >= kLimit1) || (absx < kLimit3) ) { // Code for Range 2 and Range 3
449
450 // Note 1: The above statement in the Taligent code always tests true. Logically,
451 // it was probably intended that || be &&. This would make the test logically
452 // equivalent to an "else" clause for the previous if. Replacing with an else
453 // also eliminates the problem of a code path that avoids any return statements
454 // (as detected by MrC as a compiler warning). PAF
455 else {
456
457 DeN.f = __FMADD( 64.0, absx, kBig ); // 64.0*absx + kBig;
458
459 if (absx > kLimit2) { // for large args, use recip. Range 3
460 FESETENVD(kBig+1.0); // use trunc mode momentarily
461 DeN.f = absx + kBig; // determine range of argument
462 FESETENVD(kBig+1.0); // use trunc mode momentarily
463 DeN.f = absx + kBig; // determine range of argument
464 FESETENVD(kBig); // return to round-to-nearest
465 index = DeN.i[1];
466 i = index - 1;
467 DeN.f = __FMADD( factor[i], absx, kBig + adjust[i] ); // factor[i]*absx + kBig + adjust[i];
468 }
469
470 index = DeN.i[1];
471
472 if (xHead < 0) { // use abs. value
473 arg = -xHead;
474 argl = -xTail;
475 }
476 else {
477 arg = xHead;
478 argl = xTail;
479 }
480
481 top = __FSUB( arg, atanPtr[(int32_t)index].One );
482 top2 = __FADD( top, argl );
483 topextra = arg - top - atanPtr[(int32_t)index].One;
484 top3 = top - top2 + argl + topextra; // reduced numerator
485 bot =__FMADD( arg, atanPtr[(int32_t)index].One, 1.0 ); // 1.0 + arg*atanPtr[(int32_t)index].One;
486 /* bot2 = 1.0 - bot + arg*atanPtr[(int32_t)index].One +
487 argl*atanPtr[(int32_t)index].One; */ // denominator
488 bot2 = __FMADD( argl, atanPtr[(int32_t)index].One, __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 - bot ) );
489 r = 1.0/bot;
490 quot = r*top2;
491 quot = __FMADD( __FNMSUB( bot, quot, top2), r, quot ); // quot + (top2 - bot*quot)*r;
492
493 p29 = __FNMSUB( bot, quot, top2); // top2 - quot*bot;
494 p30 = quot*bot2; // negative term
495 p33 = p29 - p30;
496 p35 = p33 + top3;
497 quot2 = p35*r;
498 quot2 = __FMADD( __FNMSUB( quot2, bot, p35 ), r, quot2 ); // quot2 + (p35 - quot2*bot)*r; // second fraction wd.
499 temp = 2.0*quot*quot2; // reduced arg is (quot,quot2)
500 argsq = __FMADD( quot, quot, temp ); // quot*quot + temp; // compute square of argument
501 argsq2 = __FMADD( quot2, quot2, __FMSUB( quot, quot, argsq) + temp ); // quot*quot - argsq + temp + quot2*quot2;
502 sumt = __FMADD( argsq, atancoeff[7], atancoeff[6] );// atancoeff[6] + argsq*atancoeff[7];
503 sum = __FMADD( argsq, sumt, atancoeff[5] ); // atancoeff[5] + argsq*sumt;
504 suml = __FMADD( argsq, sumt, atancoeff[5] - sum ); // atancoeff[5] - sum + argsq*sumt;
505
506 for (i = 4; i >= 1; i--) {
507 temps = __FMADD( sum, argsq, atancoeff[i] ); // atancoeff[i] + sum*argsq;
508 /* suml = atancoeff[i] - temps + sum*argsq + (sum*argsq2 + suml*argsq); */
509 suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
510 sum = temps;
511 }
512
513 prodlow = __FMADD( suml, argsq, sum*argsq2 ); // suml*argsq + sum*argsq2; // mult. by arg^2
514 prod = __FMADD( sum, argsq, prodlow ); // sum*argsq + prodlow;
515 prodlow = __FMSUB( sum, argsq, prod ) + prodlow; // sum*argsq - prod + prodlow;
516 temp2 = __FMADD( prodlow, quot, prod*quot2 ); // prodlow*quot + prod*quot2; // mult. by arg
517 temp = __FMADD( prod, quot, temp2 ); // prod*quot + temp2;
518 temp2 = __FMSUB( prod, quot, temp ) + temp2; // prod*quot - temp + temp2;
519 sum = temp*kRecip; // sum of series ---
520 remainder = __FNMSUB( sum, atancoeff[0], temp ); // (temp - sum*atancoeff[0]);
521 partial = __FADD( remainder, temp2 );
522 residual = remainder - partial + temp2;
523 suml = __FMADD( partial, kRecip, (residual*kRecip) ); // partial*kRecip + (residual*kRecip);
524 res = __FADD( quot, sum ); // except for argument
525 reslow = (quot - res + sum); // exact
526 resmid = __FADD( quot2, suml );
527 restiny = quot2 - resmid + suml;
528 p = __FADD( reslow, resmid ); // sum of middle terms
529 q = reslow - p + resmid; // possible residual
530 reshi = __FADD( res, p );
531 resbot = (res - reshi + p) + (q + restiny);
532 result = __FADD( atanPtr[(int32_t)index].Two, reshi );
533 reslow = __FADD( atanPtr[(int32_t)index].Three, resbot );
534 resmid = atanPtr[(int32_t)index].Two - result + reshi;
535 restiny = resbot - reslow + atanPtr[(int32_t)index].Three;
536 restemp = __FADD( resmid, reslow );
537 reshi = __FADD( result, restemp );
538
539 if (fabs(resmid) > fabs(reslow))
540 restiny = resmid - restemp + reslow + restiny;
541 else
542 restiny = reslow - restemp + resmid + restiny;
543
544 resbot = result - reshi + restemp + restiny;
545
546 if (xHead < 0.0) {
547 reshi = -reshi;
548 resbot = -resbot;
549 result = -result;
550 restemp = -restemp;
551 restiny = -restiny;
552 }
553
554 ddx.d[0] = __FADD( reshi, resbot );
555 ddx.d[1] = (reshi - ddx.d[0]) + resbot;
556
557 FESETENVD(fenv);
558 return (ddx.f); // Normal arguments done
559 }
560}
561
562// End of atanl, start of atan2l:
563
564//*******************************************************************************
565// _QATAN2 expects two quad words as arguments. It returns the arctan *
566// of the quotient of the arguments, but differs from _QATAN in that *
567// the result can be in the interval [-Pi,Pi], whereas _QATAN always *
568// returns a result in the range [-Pi/2, Pi/2]. After argument reduction *
569// and the testing of numerous special cases, _QATAN is used for the bulk *
570// of the work. *
571// *
572// atan2(NaN, anything)---------------------------------> NaN *
573// atan2(anything, NaN)---------------------------------> NaN *
574// atan2(+-0, +(anything but NaN))----------------------> +-0 *
575// atan2(+-0, -(anything but NaN))----------------------> +-Pi *
576// atan2(+-(anything but 0 and NaN), 0 )----------------> +-Pi/2 *
577// atan2(+-(anything but Infinity and NaN), +Infinity)--> +-0 *
578// atan2(+-(anything but Infinity and NaN), -Infinity)--> +-Pi *
579// atan2(+-Infinity, +Infinity)-------------------------> +-Pi/4 *
580// atan2(+-Infinity, -Infinity)-------------------------> +-3Pi/4 *
581// atan2(+-Infinity, (anything but 0, NaN, Infinity))---> +-Pi/2 *
582// *
583//*******************************************************************************
584
585// Floating-point constants
586static const double kTiny = 1.805194375864829576e-276; // 0x1.0p-916 = 0x06B00000
587// kTiny: number that still allows a 107 bit quad to consist of two normal numbers
588
589static const double kRescale = 2.923003274661805836e+48; // 0x1.0p+161
590// Needed to rescale if the denominator is "tiny"
591
592static const double kPi = 3.141592653589793116; // 0x1.921FB54442D18p1
593static const double kPiTail1 = 1.224646799147353207e-16; // 0x1.1A62633145C07p-53
594static const double kPiTail1m = 1.224646799147352961e-16; // 0x1.1A62633145C06p-53
595static const double kPiTail2 = -2.994769809718339666e-33; // 0x1.F1976B7ED8FBCp-109
596
597static const double k3PiBy4 = 2.356194490192344837; // 0x1.2D97C7F3321D2p1
598static const double k3PiBy4Tail1 = 9.184850993605148438e-17;// 0x1.A79394C9E8A0Ap-54
599
600
601long double atan2l( long double y, long double x )
602{
603 double p1, p2, p3, p2a, p4, p5, p6, p7, p8, p29, p30, p33, p35;
604 double p31, p32, p34, p36, p37, p38;
605 double quot, quot2;
606 double rarg, result, restemp;
607 double top, top2, topextra, top3;
608 double extra;
609 double bot, bot2;
610 double r, rarg2, rarg2a;
611 double temps;
612 double remainder, residual;
613 double partial;
614 double arg, argl, argsq, argsq2;
615 double temp,temp2;
616 double sum, sumt, suml;
617 double absx;
618 double xHead, xTail;
619 double res, reslow, resmid, restiny, reshi, resbot;
620 double p, q;
621 double prod, prodlow;
622 double xDen, xDenTail, yNum, yNumTail;
623 double frac3, frac2, frac, fract;
624 double pin, pin2, pin3;
625 double ressmall, restop, resint;
626 Double DeN ;
627 double fpenv;
628 DblDbl ddx, ddy;
629
630 int i, rflag;
631
632 uint32_t index;
633 uint32_t signx, signy, expx, expy, sigx, sigy;
634
635 struct atantblEntry *atanPtr = (struct atantblEntry *)atanTable - 6;
636
637 FEGETENVD(fpenv);
638 FESETENVD(0.0);
639
640 ddx.f = x; // denominator
641 ddy.f = y; // numerator
642
643 xDen = ddx.d[0]; // dden
644 xDenTail = ddx.d[1]; // ddenl
645 yNum = ddy.d[0]; // dnum
646 yNumTail = ddy.d[1]; // dnuml
647
648 expy = (ddy.i[0] >> 20) & 0x7ffu;
649 expx = (ddx.i[0] >> 20) & 0x7ffu; // get biased exponents
650 signy = (ddy.i[0] >> 31) & 1u;
651 signx = (ddx.i[0] >> 31) & 1u; // get signs
652 sigy = (ddy.i[0] & 0xfffffu) | ddy.i[1];
653 sigx = (ddx.i[0] & 0xfffffu) | ddx.i[1];
654
655 if (expy == 0x7ffu) { // y is infinite or NaN
656 if (sigy != 0) {
657 FESETENVD(fpenv);
658 return (y); // NaN y is returned
659 }
660 else if (expx == 0x7ffu) { // infinite y
661 if (sigx != 0) {
662 FESETENVD(fpenv);
663 return (x); // NaN x is returned
664 }
665 else if (signx == 0) { // infinite y and x denominator positive
666 if (signy) { // numerator negative
667 ddx.d[0] = -0.5*kPiBy2 ;
668 ddx.d[1] = -0.5*kPiBy2Tail1;
669 FESETENVD(fpenv);
670 return (ddx.f);
671 }
672 else { // numerator positive
673 ddx.d[0] = 0.5*kPiBy2 ;
674 ddx.d[1] = 0.5*kPiBy2Tail1;
675 FESETENVD(fpenv);
676 return (ddx.f);
677 }
678 }
679 else { // denominator negative
680 if (signy) { // numerator negative
681 ddx.d[0] = -k3PiBy4 ;
682 ddx.d[1] = -k3PiBy4Tail1;
683 FESETENVD(fpenv);
684 return (ddx.f);
685 }
686 else { // numerator positive
687 ddx.d[0] = k3PiBy4 ;
688 ddx.d[1] = k3PiBy4Tail1;
689 FESETENVD(fpenv);
690 return (ddx.f);
691 }
692 }
693 }
694 else { // y is Infinite , x is finite
695 if (signy == 0.0) {
696 if (signx == 0.0) { // first quadrant
697 ddx.d[0] = kPiBy2 ;
698 ddx.d[1] = kPiBy2Tail1m;
699 FESETENVD(fpenv);
700 return (ddx.f);
701 }
702 else { // second quadrant
703 ddx.d[0] = kPiBy2 ;
704 ddx.d[1] = kPiBy2Tail1;
705 FESETENVD(fpenv);
706 return (ddx.f);
707 }
708 }
709 else { // yNum < 0.0
710 if (signx == 0.0) { // fourth quadrant
711 ddx.d[0] = -kPiBy2 ;
712 ddx.d[1] = -kPiBy2Tail1m;
713 FESETENVD(fpenv);
714 return (ddx.f);
715 }
716
717 else { // third quadrant
718 ddx.d[0] = -kPiBy2 ;
719 ddx.d[1] = -kPiBy2Tail1;
720 FESETENVD(fpenv);
721 return (ddx.f);
722 }
723 }
724 }
725 } // [y is infinite or NaN]
726 if (expx == 0x7ffu) { // x is infinite or NaN
727 if (sigx != 0) { // NaN x is returned
728 FESETENVD(fpenv);
729 return (x);
730 }
731 else if (signx) { // x is -INF
732 if (signy == 0) {
733 ddx.d[0] = 2.0*kPiBy2 ;
734 ddx.d[1] = 2.0*kPiBy2Tail1m;
735 FESETENVD(fpenv);
736 return (ddx.f);
737 }
738 else {
739 ddx.d[0] = -2.0*kPiBy2 ;
740 ddx.d[1] = -2.0*kPiBy2Tail1m;
741 FESETENVD(fpenv);
742 return (ddx.f);
743 }
744 }
745 else // x is +INF
746 if (signy == 0) {
747 ddx.d[0] = 0.0 ;
748 ddx.d[1] = 0.0;
749 FESETENVD(fpenv);
750 return (ddx.f);
751 }
752 else {
753 ddx.d[0] = -0.0 ;
754 ddx.d[1] = -0.0;
755 FESETENVD(fpenv);
756 return (ddx.f);
757 }
758 } // x is infinite or NaN
759 if ((sigx | expx) == 0u) { // x is zero
760 if ((sigy | expy) == 0u) { // y is zero
761 if (signx) {
762 if (signy) {
763 ddx.d[0] = -2.0*kPiBy2 ;
764 ddx.d[1] = -2.0*kPiBy2Tail1m;
765 FESETENVD(fpenv);
766 return (ddx.f);
767 }
768 else {
769 ddx.d[0] = 2.0*kPiBy2 ;
770 ddx.d[1] = 2.0*kPiBy2Tail1m;
771 FESETENVD(fpenv);
772 return (ddx.f);
773 }
774 }
775 else {
776 FESETENVD(fpenv);
777 return (y);
778 }
779 }
780 else if (signy) { // y < 0.0
781 if (signx) {
782 ddx.d[0] = -kPiBy2 ;
783 ddx.d[1] = -kPiBy2Tail1;
784 FESETENVD(fpenv);
785 return (ddx.f);
786 }
787 else {
788 ddx.d[0] = -kPiBy2 ;
789 ddx.d[1] = -kPiBy2Tail1m;
790 FESETENVD(fpenv);
791 return (ddx.f);
792 }
793 }
794 else { // y > 0.0
795 if (signx) {
796 ddx.d[0] = kPiBy2 ;
797 ddx.d[1] = kPiBy2Tail1;
798 FESETENVD(fpenv);
799 return (ddx.f);
800 }
801 else {
802 ddx.d[0] = kPiBy2 ;
803 ddx.d[1] = kPiBy2Tail1m;
804 FESETENVD(fpenv);
805 return (ddx.f);
806 }
807 }
808 } // [x is zero]
809 if ((sigy | expy) == 0u) { // [zero y, nonzero x]
810 if (signx == 0.0) { // denominator positive
811 ddx.d[0] = yNum; // correctly signed zero result
812 ddx.d[1] = 0.0;
813 FESETENVD(fpenv);
814 return (ddx.f);
815 }
816 else { // negative denominator
817 if (signy == 0) { // numerator is 0,
818 ddx.d[0] = kPi ;
819 ddx.d[1] = kPiTail1m;
820 FESETENVD(fpenv);
821 return (ddx.f);
822 }
823 else { // numerator is -0
824 ddx.d[0] = -kPi ;
825 ddx.d[1] = -kPiTail1m;
826 FESETENVD(fpenv);
827 return (ddx.f);
828 }
829 }
830 } // [zero y, nonzero x]
831 // END OF SPECIAL CASES
832
833 // The code below is for finite y/x values, we first calculate the ratio y/x
834 // for tiny values of either x and y, and for other finite values with extra precision
835 // producing xHead, xTail and frac3. These quantities and the sign of the numerator and
836 // the denomitor are used then to produce ArcTan (y,x)
837
838 // Beginning of tiny numerator or denominator
839
840 if((((ddy.i[0]) & 0x7fffffff) < 0x06B00000)|| (((ddx.i[0]) & 0x7fffffff) < 0x06B00000)) {
841 xDen = xDen*kRescale; // For a tiny arguments,
842 xDenTail = xDenTail*kRescale; // rescale all the input
843 yNum = yNum*kRescale; // so that the full length
844 yNumTail = yNumTail*kRescale; // division can procede.
845 r = 1.0/xDen;
846 if (fabs(r) < kTiny) { // rescaling has made numerator
847 if (xDen > 0) { // too large. Result is
848 xHead = 0.0; // zero for positive denominator
849 xTail = 0.0;
850 frac3 = 0.0;
851 }
852 else { // for neg. denominator,
853 xHead = yNum; // send numerator to qatan
854 xTail = 0.0; // result will be subtracted
855 frac3 = 0.0; // from properly signed Pi later
856 }
857 }
858 }
859 // end of tiny numerator or denominator
860
861 r = 1.0/xDen; // start computation of num/den
862 // normal non-special values of x and y
863 fract = __FMUL( yNum, r );
864
865 // A detailed division (yHead,yTail)/(xHead,xTail) follows
866 frac = __FMADD( r, __FNMSUB( xDen, fract, yNum ), fract ); // fract + (yNum-xDen*fract)*r; // refine h.o. quotient
867 p29 = __FNMSUB( xDen, frac, yNum ); // yNum - xDen*frac;
868 p30 = __FMUL( frac, xDenTail ); // negative term
869 p31 = __FNMSUB( frac, xDenTail, p30 ); // p30 - frac*xDenTail;
870 p32 = p31;
871 p33 = __FSUB( p29, p30 );
872 p35 = __FADD( p33, yNumTail );
873
874 if (fabs(fract) == kInfinityu.f) {
875 xHead = fract;
876 xTail = 0.0;
877 frac3 = 0.0;
878 }
879 else if (fabs(fract)< kInfinityu.f) {
880 frac2 = __FMUL( p35, r );
881 frac2 = __FMADD( r, __FNMSUB( frac2, xDen, p35 ), frac2 ); // frac2 + (p35-frac2*xDen)*r; // mid. quot. word
882 if (fabs(p29) > fabs(p30))
883 p34 = p29 - p33 - p30 + p32;
884 else
885 p34 = p29 - (p33+p30)+p32;
886 p36 = __FNMSUB( frac2, xDen, p35 ) + __FNMSUB( frac2, xDenTail, p34 ); // p35 - frac2*xDen + p34 - frac2*xDenTail;
887 if (fabs(p33) > fabs(yNumTail))
888 p37 = p33 - p35 + yNumTail;
889 else
890 p37 = yNumTail - p35 + p33;
891 p38 = p37 + p36;
892 frac3 = p38*r/__FMADD( frac, frac, 1.0 ); // p38*r/(1.0+frac*frac); // Fraction done
893 xHead = frac; // Use ATAN routine to do
894 xTail = frac2; // most of the work
895 }
896 // Below we use the just evaluated values of xHead and xTail
897 // to first evaluate the quantity "ressmall" which add an extra precision
898 // contribution to frac3. Then, xHead, xTail, frac3, ressmall and
899 // the sign of the numerator (signy) and denominator (signx)
900 // are used to evaluate ArcTan2(y,x)
901
902 // finite result
903
904 //***************************************************************************
905 // This routine breaks the arguments into five ranges: *
906 // Range 1: 0 <= |x| <= .0859375 *
907 // Range 2: .0859375 <= |x| <= 1.015625 *
908 // Range 3: 1.015625 <= |x| < 13 *
909 // Range 4: 13 <= |x| < 2^107 *
910 // Range 5: 2^107 <= |x| <= Infinity *
911 // *
912 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> *
913 // *
914 // In Range 1, a 15th degree polynomial in x^2 is used to evaluate QATAN. *
915 // *
916 // In Range 2 and, *
917 // In Range 3 the accurate table method is used to range-reduce the *
918 // argument to a small value to be used with an interpolating *
919 // polynomial for QATAN. In Ranges 2 and 3, the resultant *
920 // argument is less than 2^-7, so an 8th degree polynomial can *
921 // be used. Range 3 requires additional work to find the *
922 // correct index into the accurate table. *
923 // *
924 // In Range 4, QATAN(x) = Pi/2 - QATAN(1/x). After taking the reciprocal *
925 // of the argument, the quotient is less than 1/13, so that the *
926 // 15th degree polynomial used for Range 1 is used to compute *
927 // QATAN(1/x) *
928 // *
929 // In Range 5, the correctly rounded value differs from Pi/2 by less than *
930 // 1/2 ulp. It is, however, incorrect to let the value become *
931 // correctly rounded, as that correctly rounded value exceeds *
932 // PI/2 and is in the wrong quadrant. Applying the inverse *
933 // function, QTAN, to the correctly rounded value of PI/2 would *
934 // produce a result of the opposite sign to the argument *
935 // originally given to ATAN. *
936 //***************************************************************************
937
938
939 // THIS CODE ONLY TREATS RANGES 1 TROUGH 5. WE ASSUME THAT xHead and xTail
940 // are such that we are ALWAYS in Range 1-5.
941
942 absx = fabs(xHead);
943
944 // Range 5
945
946 if (absx >= kLimit4) { // Filter out unusual arguments
947 reshi = kPiBy2; // if the second argument is
948 if (signx == 0.0) { // positive, the result is in
949 resbot = kPiBy2Tail1m; // first or fourth quadrant
950 ressmall = 0; // set up appropriate value of
951 if (xHead < 0.0) { // pi/2, truncated
952 reshi = -reshi;
953 resbot = -resbot;
954 }
955 }
956 else { // Otherwise, result is in the
957 reshi = kPiBy2; // second or third quadrant.
958 resbot = kPiBy2Tail1; // set up pi/2, correctly rounded
959 ressmall = kPiBy2Tail2;
960 if (xHead < 0.0) {
961 reshi = -reshi;
962 resbot = -resbot;
963 ressmall = -ressmall;
964 }
965 }
966 }
967
968 // Ranges 1 and 4 (DO NOT USE TABLE)
969
970 rflag = 0;
971
972 // Range 4
973
974 if ((absx >= kLimit3) && (absx < kLimit4) ) { // Range 4
975 rarg = 1.0/xHead;
976 rflag = 1; // indicate large argument
977 p1 = __FNMSUB( xHead, rarg, 1.0 ); // 1.0 - xHead*rarg; // time to complete hi precision
978 p2 = xTail*rarg; // reciprocal of argument, to
979 p3 = p1 - p2; // use as reduced argument.
980 rarg2a = p3*rarg;
981 rarg2 = rarg2a + __FNMSUB( xHead, rarg2a, p3 )*rarg; // rarg2a + (p3 - xHead*rarg2a)*rarg;
982 xHead = rarg;
983 xTail = rarg2;
984 }
985
986 // Ranges 1 and 4 only
987
988 if( (absx < kLimit1) || ((absx >= kLimit3) && (absx < kLimit4) )) { //Code for Range 1 and Range 4
989 temp = 2.0*xHead*xTail; // direct evaluation of series
990 argsq = __FMADD( xHead, xHead, temp ); // xHead*xHead + temp; // compute square of argument
991 argsq2 = __FMADD( xTail, xTail, __FMSUB( xHead, xHead, argsq ) + temp ); // xHead*xHead - argsq + temp + xTail*xTail;
992 sum = atancoeff[14];
993 sum = __FMADD( argsq, sum, atancoeff[13] ); // atancoeff[13] + argsq*sum;
994 sum = __FMADD( argsq, sum, atancoeff[12] ); // atancoeff[12] + argsq*sum;
995 sum = __FMADD( argsq, sum, atancoeff[11] ); // atancoeff[11] + argsq*sum;
996 sumt = __FMADD( argsq, sum, atancoeff[10] ); // atancoeff[10] + argsq*sum;
997 sum = __FMADD( argsq, sum, atancoeff[9] ); // atancoeff[9] + argsq*sumt;
998 suml = __FMADD( argsq, sumt, atancoeff[9] - sum ); // atancoeff[9] - sum + argsq*sumt;
999
1000 for (i = 8; i >= 1; i--) {
1001 temps = __FMADD( argsq, sum, atancoeff[i] ); // atancoeff[i] + sum*argsq;
1002 /* suml = atancoeff[i] - temps+sum*argsq + (sum*argsq2 + suml*argsq); */
1003 suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
1004 sum=temps;
1005 }
1006
1007 prodlow = __FMADD( suml, argsq, sum*argsq2 ); // suml*argsq + sum*argsq2; // multiply by arg^2
1008 prod = __FMADD( sum, argsq, prodlow ); // sum*argsq + prodlow;
1009 prodlow = __FMSUB( sum, argsq, prod ) + prodlow; // sum*argsq - prod + prodlow;
1010 temp2 = __FMADD( prodlow, xHead, prod*xTail ); // prodlow*xHead + prod*xTail; // multiply by arg
1011 temp = __FMADD( prod, xHead, temp2 ); // prod*xHead + temp2;
1012 temp2 = __FMSUB( prod, xHead, temp ) + temp2; // prod*xHead - temp + temp2;
1013 sum = temp*kRecip; // sum of series ---
1014 remainder = __FNMSUB( sum, atancoeff[0], temp ); // (temp - sum*atancoeff[0]);
1015 partial = __FADD( remainder, temp2 );
1016 residual = remainder - partial + temp2;
1017 suml = __FMADD( partial, kRecip, (residual*kRecip) );// partial*kRecip + (residual*kRecip);
1018 res = __FADD( xHead, sum ); // except for argument
1019 reslow = (xHead - res + sum); // exact
1020 resmid = __FADD( xTail, suml );
1021 restiny = xTail - resmid + suml;
1022 p = __FADD( reslow, resmid ); // sum of middle terms
1023 q = reslow - p + resmid; // possible residual
1024 reshi = res + p;
1025 resbot = (res - reshi + p) + (q + restiny);
1026 if (rflag == 1) { // large magnitude argument
1027 if (xHead > 0) { // fetch pi/2 with proper sign
1028 p1 = kPiBy2;
1029 p2 = kPiBy2Tail1;
1030 p2a = kPiBy2Tail2;
1031 }
1032 else {
1033 p1 = -kPiBy2;
1034 p2 = -kPiBy2Tail1;
1035 p2a = -kPiBy2Tail2;
1036 }
1037 p3 = __FSUB( p1, reshi ); // subtract result from pi/2
1038 p4 = p1 - p3 - reshi;
1039 p5 = __FSUB( p2, resbot );
1040 p6 = p2 - p5 - resbot + p2a;
1041 p7 = __FADD( p4, p5 );
1042 reshi = p3 + p7;
1043 if (fabs(p4) > fabs(p5))
1044 p8 = p4 - p7 + p5;
1045 else
1046 p8 = p5 - p7 + p4;
1047 resbot = __FADD( p3 - reshi + p7, (p6+p8-restiny) );
1048 ressmall = (p3 - reshi + p7) - resbot + (p6 + p8 - restiny);
1049 }
1050 else //if rflag = 0 i.e., Range 1
1051 ressmall = (res - reshi + p) - resbot + (q + restiny);
1052 }
1053
1054 // Ranges 2 and 3 (USE TABLE)
1055
1056 if( (absx >= kLimit1) && (absx < kLimit3) ) { // Code for Range 2 and Range 3
1057 DeN.f = __FMADD( 64.0, absx, kBig ); // 64.0 *absx+kBig;
1058 if (absx > kLimit2) { // for large args, use recip. Range 3
1059 FESETENVD(kBig+1.0); // use trunc mode momentarily
1060 DeN.f = absx + kBig; // determine range of argument
1061 FESETENVD(kBig+1.0); // use trunc mode momentarily
1062 DeN.f = absx + kBig; // determine range of argument
1063 FESETENVD(kBig); // return to round-to-nearest
1064 index = DeN.i[1];
1065 i = index - 1;
1066 DeN.f = __FMADD( factor[i], absx, kBig+adjust[i] ); // factor[i]*absx+kBig+adjust[i];
1067 }
1068 index = DeN.i[1];
1069 if (xHead < 0) { // use abs. value
1070 arg = -xHead;
1071 argl = -xTail;
1072 }
1073 else {
1074 arg = xHead;
1075 argl = xTail;
1076 }
1077 top = __FSUB( arg, atanPtr[(int32_t)index].One );
1078 top2 = __FADD( top, argl );
1079 topextra = arg - top - atanPtr[(int32_t)index].One;
1080 top3 = top - top2 + argl + topextra; // reduced numerator
1081 bot = __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 ); // 1.0 + arg*atanPtr[(int32_t)index].One;
1082 /* bot2 = 1.0 - bot + arg*atanPtr[(int32_t)index].One + argl*atanPtr[(int32_t)index].One; */ // denominator
1083 bot2 = __FMADD( argl, atanPtr[(int32_t)index].One, __FMADD( arg, atanPtr[(int32_t)index].One, 1.0 - bot ) );
1084 r = 1.0/bot;
1085 quot = r*top2;
1086 quot = __FMADD( __FNMSUB( bot, quot, top2 ) , r, quot ); // quot + (top2-bot*quot)*r;
1087
1088 p29 = __FNMSUB( bot, quot, top2 ); //top2 - quot*bot;
1089 p30 = quot*bot2; // negative term
1090 p33 = p29-p30;
1091 p35 = p33 + top3;
1092 quot2 = p35*r;
1093 quot2 = __FMADD( __FNMSUB( quot2, bot, p35 ), r, quot2 ); // quot2 + (p35-quot2*bot)*r; // second fraction wd.
1094 temp = 2.0*quot*quot2; // reduced arg is (quot,quot2)
1095 argsq = __FMADD( quot, quot, temp ); // quot*quot + temp; // compute square of argument
1096 argsq2 = __FMADD( quot2, quot2, __FMSUB( quot, quot, argsq) + temp ); // quot*quot - argsq + temp + quot2*quot2;
1097 sumt = __FMADD( argsq, atancoeff[7], atancoeff[6] );// atancoeff[6] + argsq*atancoeff[7];
1098 sum = __FMADD( argsq, sumt, atancoeff[5] ); // atancoeff[5] + argsq*sumt;
1099 suml = __FMADD( argsq, sumt, atancoeff[5] - sum ); // atancoeff[5] - sum + argsq*sumt;
1100
1101 for (i = 4; i >= 1; i--) {
1102 temps = __FMADD( sum, argsq, atancoeff[i] ); // atancoeff[i] + sum*argsq;
1103 /* suml = atancoeff[i] - temps + sum*argsq + (sum*argsq2 + suml*argsq); */
1104 suml = __FMADD( sum, argsq, atancoeff[i] - temps ) + __FMADD( sum, argsq2, suml*argsq );
1105 sum = temps;
1106 }
1107
1108 prodlow = __FMADD( suml, argsq, sum*argsq2 ); // suml*argsq + sum*argsq2; // mult. by arg^2
1109 prod = __FMADD( sum, argsq, prodlow ); // sum*argsq + prodlow;
1110 prodlow = __FMSUB( sum, argsq, prod ) + prodlow; // sum*argsq - prod + prodlow;
1111 temp2 = __FMADD( prodlow, quot, prod*quot2 ); // prodlow*quot + prod*quot2; // mult. by arg
1112 temp = __FMADD( prod, quot, temp2 ); // prod*quot + temp2;
1113 temp2 = __FMSUB( prod, quot, temp ) + temp2; // prod*quot - temp + temp2;
1114 sum = temp*kRecip; // sum of series ---
1115 remainder = __FNMSUB( sum, atancoeff[0], temp ); // (temp - sum*atancoeff[0]);
1116 partial = __FADD( remainder, temp2 );
1117 residual = remainder - partial + temp2;
1118 suml = __FMADD( partial, kRecip, (residual*kRecip) ); // partial*kRecip + (residual*kRecip);
1119 res = __FADD( quot, sum ); // except for argument
1120 reslow = (quot - res + sum); // exact
1121 resmid = __FADD( quot2, suml );
1122 restiny = quot2 - resmid + suml;
1123 p = __FADD( reslow, resmid ); // sum of middle terms
1124 q = reslow - p + resmid; // possible residual
1125 reshi = __FADD( res, p );
1126 resbot = (res - reshi + p) + (q + restiny);
1127 result = __FADD( atanPtr[(int32_t)index].Two, reshi );
1128 reslow = __FADD( atanPtr[(int32_t)index].Three, resbot );
1129 resmid = atanPtr[(int32_t)index].Two - result + reshi;
1130 restiny = resbot - reslow + atanPtr[(int32_t)index].Three;
1131 restemp = __FADD( resmid, reslow );
1132 reshi = __FADD( result, restemp );
1133
1134 if (fabs(resmid) > fabs(reslow))
1135 restiny = resmid - restemp + reslow + restiny;
1136 else
1137 restiny = reslow - restemp + resmid + restiny;
1138 resbot = result - reshi + restemp + restiny;
1139 if (xHead < 0.0) {
1140 reshi = -reshi;
1141 resbot = -resbot;
1142 result = -result;
1143 restemp = -restemp;
1144 restiny = -restiny;
1145 }
1146 ressmall = (result - reshi + restemp) - resbot + restiny;
1147 }
1148 extra = ressmall + frac3;
1149 if (signx == 0) { // For Positive denominator
1150 ddx.d[0] = __FADD( reshi, (resbot + extra) );
1151 ddx.d[1] = (reshi - ddx.d[0]) + (resbot + extra);
1152 }
1153 else {
1154 if (signy == 0) { // For Negative denominator
1155 pin = kPi; // return 2nd|3rd quadrand angle,
1156 pin2 = kPiTail1; // depending on sign of numerator
1157 pin3 = kPiTail2;
1158 }
1159 else {
1160 pin = -kPi; // --> Pi+atan, for dnum < 0,
1161 pin2 = -kPiTail1; // atan-Pi, for dnum > 0
1162 pin3 = -kPiTail2;
1163 }
1164 restop = __FADD( pin, reshi );
1165 resmid = pin - restop + reshi;
1166 reslow = __FADD( pin2, resbot );
1167 restiny = pin2 - reslow + resbot + extra;
1168 resint = __FADD( resmid, reslow );
1169 reshi = __FADD( restop, resint );
1170 resbot = restop - reshi+resint;
1171 if (fabs(resmid) > fabs(reslow))
1172 resbot = resbot + (resmid-resint+reslow+restiny);
1173 else
1174 resbot = resbot + (reslow-resint+resmid+restiny);
1175 ddx.d[0] = __FADD( reshi, resbot );
1176 ddx.d[1] = (reshi - ddx.d[0]) + resbot;
1177 }
1178
1179 FESETENVD(fpenv);
1180 return (ddx.f);
1181}
1182#endif
1183