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// SinCosTanDD.c
24//
25// Double-double Function Library
26// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
27//
28// long double sinl( long double x );
29// long double cosl( long double x );
30// long double tanl( long double x );
31//
32// Change History: [7d]
33//
34// 06/20/96 PAF - Remove /* */ comments (for color editor)
35
36#include "math.h"
37#include "fp_private.h"
38#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
39#include "DD.h"
40
41static const float k1_8To52 = 6755399441055744.0f; // 0x1.8p52
42static const double kPiBy2 = 1.570796326794896619231322; // 0x1.921FB54442D18p0
43static const double kPiBy2Tail1 = 6.1232339957367660e-17; // 0x1.1A62633145C07p-54
44static const double kPiBy2Tail2 = -1.4973849048591698e-33; // 0x1.F1976B7ED8FBCp-110
45static const double kPiBy2Tail3 = 5.5622711043168264e-50; // 0x1.4CF98E804177Dp-164
46static const double k2ByPi = 0.636619772367581382; // 0x1.45f306dc9c883p-1
47static const float k2To52 = 4503599627370496.0f; // 0x1.0p+52
48static const float k2To54 = 18014398509481984.0f; // 0x1.0p+54
49static const double k2To45m1 = 35184372088831.0; // 0x1.fffffffffffffp44
50
51static const uint32_t CosCoeff[] __attribute__ ((aligned(8))) = {
52 0xbFE00000, 0x00000000, 0x39685556, 0x00000000,
53 0x3FA55555, 0x55555555, 0x3C455555, 0x51E8C72E,
54 0xbF56C16C, 0x16C16C17, 0x3BF02698, 0xA05BBBBD,
55 0x3EFA01A0, 0x1A0169EB, 0x3B540FCD, 0xE3AE1BE0,
56 0xbE927E4D, 0x648A174F, 0xbB361705, 0xD6F2C25B
57};
58
59static const uint32_t SinCoeff[] __attribute__ ((aligned(8))) = {
60 0xbFC55555, 0x55555555, 0xbC655555, 0x55555552,
61 0x3F811111, 0x11111111, 0x3C011111, 0x0DBB7C30,
62 0xbF2A01A0, 0x1A01A01A, 0xbB510055, 0x358E4100,
63 0x3EC71DE3, 0xA556A5C2, 0x3B61A976, 0xF8644324,
64 0xbE5AE642, 0x52F2B067, 0xbAD356A3, 0xF6278F6C
65};
66
67static const uint32_t CosCoeff16[] __attribute__ ((aligned(8))) = {
68 0xBFE00000, 0x00000000, 0x395DA200, 0x00000000,
69 0x3FA55555, 0x55555555, 0x3C455555, 0x554FC86B,
70 0xBF56C16C, 0x16C16C17, 0x3BEF49FE, 0x788C8E18,
71 0x3EFA01A0, 0x1A01A01A, 0xBB6A022F, 0xCA316020,
72 0xBE927E4F, 0xB7789A12, 0xBB245398, 0xA829E8F5,
73 0x3E21EED8, 0xEE008FB4, 0x3AC49E38, 0xCA6157DE,
74 0xBDA938BA, 0xA889516E, 0x3A46CBC1, 0x4AF6E840
75};
76
77static const uint32_t SinCoeff16[] __attribute__ ((aligned(8))) = {
78 0xBFC55555, 0x55555555, 0xBC655555, 0x55555554,
79 0x3F811111, 0x11111111, 0x3C011111, 0x10F8EB9C,
80 0xBF2A01A0, 0x1A01A01A, 0xBB69FF8C, 0xFBC69E00,
81 0x3EC71DE3, 0xA556C733, 0x3B6EFC3D, 0xB602439A,
82 0xBE5AE645, 0x67F53972, 0xBAF653B7, 0xA917A91D,
83 0x3DE61246, 0x1048C848, 0x3A7CFCE5, 0x9E3CEA38,
84 0xBD6AE6FF, 0xCFCF788E, 0x39DA11EB, 0x64C11650
85};
86
87static const uint32_t TrigTable[] __attribute__ ((aligned(8))) = {
88 0x3FB80000, 0x00054E6C, // angle
89 0x3FB7F701, 0x032A9959, 0xBBAD11F0, 0x5495F5BE, // Sin Head, Sin Tail
90 0x3FEFDC06, 0xBF7E5BB6, 0xBBDCB529, 0x1917F935, // Cos Head, Cos Tail
91 0x3FBC0000, 0x000A5F0A,
92 0x3FBBF1B7, 0x8572884A, 0x3BBAADFC, 0xF7E52ACA,
93 0x3FEFCF0C, 0x800E7577, 0xBBE4670B, 0xB43517A0,
94 0x3FC00000, 0x000F1278,
95 0x3FBFEAAE, 0xEEA4D6F0, 0xBB9D8153, 0x7343A887,
96 0x3FEFC015, 0x527CE390, 0x3BE16998, 0xA0C1E22D,
97 0x3FC20000, 0x00052790,
98 0x3FC1F0D3, 0xD7B4E938, 0xBBB94CF6, 0x02B38414,
99 0x3FEFAF22, 0x263C1D96, 0xBBDE4F5F, 0x508C4D7E,
100 0x3FC40000, 0x000143F6,
101 0x3FC3EB31, 0x2C5EA6CF, 0x3BBF0AD1, 0xD6EF035D,
102 0x3FEF9C34, 0x0A7CB78E, 0xBBE87189, 0x89E0691A,
103 0x3FC60000, 0x0007E6FF,
104 0x3FC5E44F, 0xCFA8F023, 0xBBAF148F, 0x61A85B8E,
105 0x3FEF874C, 0x2E1E9676, 0x3BC38F0F, 0xAB6185B5,
106 0x3FC80000, 0x000E77F8,
107 0x3FC7DC10, 0x2FC929C2, 0xBBC3C05A, 0x114DA314,
108 0x3FEF706B, 0xDF9E2180, 0xBBD8BEA1, 0x93F3F8C7,
109 0x3FCA0000, 0x00102F36,
110 0x3FC9D252, 0xD0DE9D1A, 0x3B8BF04E, 0xADF4B6FE,
111 0x3FEF5794, 0x8CFE96A3, 0xBBE58FFE, 0x7BB60718,
112 0x3FCC0000, 0x0009CE78,
113 0x3FCBC6F8, 0x4EE5F43E, 0xBBA7E35C, 0xC27E82D6,
114 0x3FEF3CC7, 0xC3B3493B, 0xBBCC7CA2, 0xE7033121,
115 0x3FCE0000, 0x00150434,
116 0x3FCDB9E1, 0x5FCA16EB, 0x3BC107C0, 0xD41B55C5,
117 0x3FEF2007, 0x30852C41, 0x3BE70655, 0x75B2D6FB,
118 0x3FD00000, 0x0061B811,
119 0x3FCFAAEE, 0xD5B07238, 0xBBCB82BA, 0x01ECEAD1,
120 0x3FEF0154, 0x9F71D818, 0x3BE95043, 0xCC0B6426,
121 0x3FD10000, 0x0020111B,
122 0x3FD0CD00, 0xCF125569, 0xBBC2D015, 0xDA2CFEAE,
123 0x3FEEE0B1, 0xFBBCBB9D, 0xBBE9AA1F, 0xF00E2374,
124 0x3FD20000, 0x000C1AD3,
125 0x3FD1C37D, 0x64D25988, 0x3BC4B236, 0xE14BCAA2,
126 0x3FEEBE21, 0x4F75419B, 0x3BE52928, 0x2D338F16,
127 0x3FD30000, 0x000E7F89,
128 0x3FD2B8DD, 0xC44C91CC, 0xBBC5707F, 0x432D5073,
129 0x3FEE99A4, 0xC3A5AEA3, 0x3BD23919, 0x2EA440E3,
130 0x3FD40000, 0x0017A801,
131 0x3FD3AD12, 0x9780568B, 0xBBC59FD1, 0x16DF09C5,
132 0x3FEE733E, 0xA0159A53, 0xBBB12BF5, 0xDB2DCDCC,
133 0x3FD50000, 0x0039F7ED,
134 0x3FD4A00C, 0x9B461D51, 0xBBCDFE0E, 0x1445AAA9,
135 0x3FEE4AF1, 0x4B20ED63, 0x3BD374F6, 0x9F86063A,
136 0x3FD60000, 0x0000A1AA,
137 0x3FD591BC, 0x9FA38DCC, 0xBBCD07DB, 0x418C1954,
138 0x3FEE20BF, 0x49ACBB83, 0x3B864FA4, 0xEB634A4E,
139 0x3FD70000, 0x00011BC1,
140 0x3FD68213, 0x8A39E197, 0xBBBA6968, 0x3C12FB6D,
141 0x3FEDF4AB, 0x3EBD5578, 0xBBD6766E, 0x746C48E7,
142 0x3FD80000, 0x0008ED50,
143 0x3FD77102, 0x557E9094, 0x3BB488F1, 0x4A8DF658,
144 0x3FEDC6B7, 0xEB97B68C, 0x3BDCA325, 0xE2462D16,
145 0x3FD90000, 0x00861E6C,
146 0x3FD85E7A, 0x12FE6D57, 0x3BA65076, 0x5FCB7841,
147 0x3FED96E8, 0x2F58212A, 0xBBEC65EC, 0x2CEB0E13,
148 0x3FDA0000, 0x00021040,
149 0x3FD94A6B, 0xE9F72C06, 0x3BA04A4B, 0xE0FC13D6,
150 0x3FED653F, 0x073DD7E0, 0x3BEC1A4A, 0xF1A7037E,
151 0x3FDB0000, 0x005EF712,
152 0x3FDA34C9, 0x1D1BB055, 0x3BC42083, 0xD43A5519,
153 0x3FED31BF, 0x8D7A0AAD, 0xBB9D5F11, 0x5DAECA9B,
154 0x3FDC0000, 0x000C075C,
155 0x3FDB1D83, 0x053CFB6A, 0xBBDA72E9, 0x9A1C2FBE,
156 0x3FECFC6C, 0xFA50214C, 0xBBCB5164, 0xBC69A105,
157 0x3FDD0000, 0x0056E33E,
158 0x3FDC048B, 0x17FF5F2B, 0x3BD8AF58, 0xF8525AD2,
159 0x3FECC54A, 0xA29F9263, 0x3BD22EA3, 0xDBB55D0D,
160 0x3FDE0000, 0x005CCC54,
161 0x3FDCE9D2, 0xE4276EF1, 0xBBAAEDF9, 0xDD5DC899,
162 0x3FEC8C5B, 0xF8B9244D, 0xBBC2D9C6, 0xCE4671F4,
163 0x3FDF0000, 0x00161118,
164 0x3FDDCD4C, 0x154623DC, 0xBBD97F5B, 0x8DF7FEBF,
165 0x3FEC51A4, 0x8B85F41A, 0x3BC354E8, 0x86ABC9CE,
166 0x3FE00000, 0x0013AF20,
167 0x3FDEAEE8, 0x746D926F, 0x3BD2B7DF, 0xA4EEC7DD,
168 0x3FEC1528, 0x06520D6D, 0x3BE9EE87, 0x2601172E,
169 0x3FE08000, 0x0046DF51,
170 0x3FDF8E99, 0xE7E60D68, 0x3BAF1859, 0x1DD2BAB4,
171 0x3FEBD6EA, 0x30DFA2E0, 0x3BE95328, 0x8BAB7ACB,
172 0x3FE10000, 0x0002BED1,
173 0x3FE03629, 0x39C8F748, 0x3BCFAF90, 0x9F47E8F6,
174 0x3FEB96EE, 0xEF572000, 0x3BD89C91, 0x79738697,
175 0x3FE18000, 0x00404634,
176 0x3FE0A402, 0x1ED4F66C, 0x3BD4E9FD, 0x2527C4F7,
177 0x3FEB553A, 0x40EAA3C8, 0xBBE9F185, 0x9F0B73D8,
178 0x3FE20000, 0x0005ABAA,
179 0x3FE110D0, 0xC4BB683B, 0xBBB78FD4, 0x4D302311,
180 0x3FEB11D0, 0x415F9E99, 0x3BDAE453, 0xDBCB3704,
181 0x3FE28000, 0x002D02D5,
182 0x3FE17C8E, 0x5F549FEE, 0x3B6368B3, 0xE72D341D,
183 0x3FEACCB5, 0x26DE0531, 0xBB8CE182, 0xC27CD7B8,
184 0x3FE30000, 0x0002B1A8,
185 0x3FE1E734, 0x323892EB, 0x3BC9A07C, 0xF3B89184,
186 0x3FEA85ED, 0x43725E55, 0x3BE2C483, 0xBE536F66,
187 0x3FE38000, 0x000730F9,
188 0x3FE250BB, 0x937E7157, 0xBBC3FA7C, 0xADF21ED4,
189 0x3FEA3D7D, 0x034EA01E, 0xBBD5786E, 0x86ABB722,
190 0x3FE40000, 0x00076EEF,
191 0x3FE2B91D, 0xEA8E4953, 0x3BE1CB0F, 0x13421234,
192 0x3FE9F368, 0xED8CD61E, 0xBBC0A358, 0x32653793,
193 0x3FE48000, 0x00119ECA,
194 0x3FE32054, 0xB156DCB6, 0xBBBF1C5F, 0x95EB1A26,
195 0x3FE9A7B5, 0xA35FDCFF, 0x3BE6C143, 0x88B5E209,
196 0x3FE50000, 0x000045B6,
197 0x3FE38659, 0x74565F66, 0x3BD5144E, 0x65AA2BBD,
198 0x3FE95A67, 0xE00C8774, 0x3BC1A72A, 0xE0D050C3,
199 0x3FE58000, 0x00A4E7B8,
200 0x3FE3EB25, 0xD3EDE59C, 0x3BE1F420, 0x5560DD7C,
201 0x3FE90B84, 0x77E73599, 0x3B9250E3, 0x5B17FE76,
202 0x3FE60000, 0x000049F4,
203 0x3FE44EB3, 0x81CF7192, 0xBBD5B364, 0x7D4D895F,
204 0x3FE8BB10, 0x5A5D9A12, 0xBBD2BB98, 0xF4F2C91C,
205 0x3FE68000, 0x0008C197,
206 0x3FE4B0FC, 0x46B16552, 0x3BC863B5, 0xB735A540,
207 0x3FE86910, 0x8D71FD5A, 0x3BE5C4E3, 0x402BA3EE,
208 0x3FE70000, 0x00340EE5,
209 0x3FE511F9, 0xFDA26352, 0x3BCDC369, 0x16E2A8DB,
210 0x3FE8158A, 0x316F2658, 0x3BD69C38, 0x4E2E4085,
211 0x3FE78000, 0x005756FD,
212 0x3FE571A6, 0x96AE2DA7, 0x3BE24CA0, 0xA2E35E02,
213 0x3FE7C082, 0x7ECF5E07, 0xBBCBD891, 0x2A83004E,
214 0x3FE80000, 0x004601DC,
215 0x3FE5CFFC, 0x16F2C847, 0x3BCA6614, 0x963A9B5C,
216 0x3FE769FE, 0xC62568E3, 0x3BE746CD, 0xBABDB543,
217 0x3FE88000, 0x002EF948,
218 0x3FE62CF4, 0x99438A17, 0x3BA3ED99, 0x63704218,
219 0x3FE71204, 0x6F86E919, 0xBBE0671D, 0x4098C629,
220 0x3FE90000, 0x0022F620,
221 0x3FE6888A, 0x4E2C1E13, 0xBBC0D332, 0x374F8A22,
222 0x3FE6B898, 0xFA865CFA, 0xBBDE0573, 0x1D87548D,
223 0x3FE98000, 0x00883EB0,
224 0x3FE6E2B7, 0x7C9FF82A, 0x3BE241BC, 0xA153A6CC,
225 0x3FE65DC1, 0xFD8A1C59, 0xBBDBAD2A, 0x038DFB75
226};
227
228static const uint32_t TanCoeff16[] __attribute__ ((aligned(8))) = {
229 0x3FD55555, 0x55555555, 0x3C755555, 0x5555554C,
230 0x3FC11111, 0x11111111, 0x3C411111, 0x11327FD4,
231 0x3FABA1BA, 0x1BA1BA1C, 0xBC47917B, 0x8730ECE8,
232 0x3F9664F4, 0x882C10FA, 0xBC08946C, 0x04337BC0,
233 0x3F8226E3, 0x55E6C239, 0xBC16F3F8, 0x5EB82A9C,
234 0x3F6D6D3D, 0x0E159BD4, 0xBC0BF02E, 0x91A39520,
235 0x3F57DA36, 0x44E2BF4A, 0x3BF4CC81, 0xFDF36B54,
236 0x3F435582, 0xB4E9E8D9, 0xBBEAAB03, 0xE323ED60,
237 0x3F2F5712, 0xDCD1E82C, 0x3BA67F16, 0xB4DC7710,
238 0x3F19C9D2, 0x46D3C681, 0x3BA2B1F2, 0xBC709028
239};
240
241static const uint32_t CotCoeff16[] __attribute__ ((aligned(8))) = {
242 0x3FD55555, 0x55555555, 0x3C755555, 0x55555550,
243 0x3F96C16C, 0x16C16C17, 0xBC2F49F4, 0x9F46AA0E,
244 0x3F61566A, 0xBC011567, 0xBC050FFB, 0x0703E386,
245 0x3F2BBD77, 0x9334EF0B, 0xBBC541A3, 0x5FEBE630,
246 0x3EF66A8F, 0x2BF70EDF, 0xBB960EEE, 0x58853E10,
247 0x3EC22805, 0xD642B849, 0x3B49A324, 0x92C33F58,
248 0x3E8D6DB2, 0xD58EAF82, 0x3B103D9B, 0xCE4CCBB8,
249 0x3E57DA1C, 0x0FF6D02F, 0x3AFBDB54, 0xE9F15B6C,
250 0x3E2396A6, 0x6CC1B8E2, 0x3AC04FB7, 0x0AEC86D2
251};
252
253static const uint32_t TanCoeff[] __attribute__ ((aligned(8))) = {
254 0x3ff00000, 0x00000000, 0x00000000, 0x00000000,
255 0x3FD55555, 0x55555555, 0x3C755555, 0x55555549,
256 0x3FC11111, 0x11111111, 0x3C411111, 0x17771B40,
257 0x3FABA1BA, 0x1BA1BA1C, 0xBC47A471, 0x53F3999C,
258 0x3F9664F4, 0x882C119E, 0x3C2A931E, 0x61CFDED8,
259 0x3F8226E3, 0x5546FD80, 0xBC22700A, 0xDC5E76C2,
260 0x3F6D6DCB, 0xC6BCC1FD, 0x3C0EB3DD, 0x1EBC7C6C
261};
262
263static const uint32_t TanTable[] __attribute__ ((aligned(8))) = {
264 0x3FB80000, 0x00095B76,
265 0x3FB81210, 0x420B0DDD, 0xBBADBC24, 0xD03EE61B,
266 0x40254552, 0xEE62D182, 0xBC23A778, 0xD4D4C82D,
267 0x3FBC0000, 0x00F6D8DC,
268 0x3FBC1CB8, 0x85A84FC9, 0xBBA0734F, 0xD6E646D7,
269 0x40223676, 0x163ABE75, 0x3C19C932, 0xFC451F09,
270 0x3FC00000, 0x002177FC,
271 0x3FC01577, 0xAF3710E9, 0xBBBB1D68, 0x25554EB9,
272 0x401FD549, 0xF004A93B, 0x3BF1B140, 0x750F85AD,
273 0x3FC20000, 0x0000C2F1,
274 0x3FC21E9E, 0x0175E475, 0x3BBD91B3, 0xA31924A6,
275 0x401C41B6, 0xE169DC8D, 0x3C057A49, 0x2ABEE9CF,
276 0x3FC40000, 0x004D703E,
277 0x3FC42A13, 0xDFCB15A0, 0xBB94ED3B, 0x3620DA54,
278 0x4019642D, 0xFDBA3F8F, 0x3BBB47D8, 0x80CDA3C6,
279 0x3FC60000, 0x000CAC24,
280 0x3FC6381F, 0x200F2AED, 0x3BBF3967, 0x0381AC58,
281 0x40170B09, 0x205DE76C, 0xBBFEFF3A, 0x14B71AB8,
282 0x3FC80000, 0x00473296,
283 0x3FC84906, 0xF15CE818, 0x3BB03E11, 0x36E0F2D3,
284 0x4015152E, 0xCDA71D37, 0xBBE77628, 0xDB83D025,
285 0x3FCA0000, 0x00B2D7A8,
286 0x3FCA5D14, 0x0081E4D7, 0xBBC0188C, 0x9C6AB2DE,
287 0x40136BB4, 0xBA046E82, 0xBC0F80BA, 0xA34DE4A8,
288 0x3FCC0000, 0x0066CDE2,
289 0x3FCC7490, 0xA23DC3B4, 0xBBCB362F, 0xA69FE4DA,
290 0x4011FE3C, 0xA58F0345, 0xBBFCB31C, 0x2F2152B4,
291 0x3FCE0000, 0x00252978,
292 0x3FCE8FC9, 0x01177F3A, 0xBBA5E8BA, 0x10809B66,
293 0x4010C0C5, 0xABFB3524, 0xBBFA678B, 0xCA50DEEC,
294 0x3FD00000, 0x00638D21,
295 0x3FD05785, 0xA4A65715, 0x3BC29365, 0x7FDA2145,
296 0x400F549E, 0x3203FCA2, 0x3BE2971C, 0x11C2A8A3,
297 0x3FD10000, 0x002C35A8,
298 0x3FD16953, 0xEACF2DA7, 0x3BC40A91, 0xBA0FE60F,
299 0x400D67EC, 0xF327EFBB, 0x3C05D1EC, 0x4B5FCFCB,
300 0x3FD20000, 0x001C3AB2,
301 0x3FD27D78, 0xB42A0CE9, 0x3BCCED21, 0x016F79E8,
302 0x400BB0C1, 0xF1319BD2, 0xBC051F68, 0x2CBAE56E,
303 0x3FD30000, 0x000F9C12,
304 0x3FD3941E, 0xADA8C533, 0xBBBB6B83, 0x08722CFA,
305 0x400A26A8, 0xA4D1FD5D, 0xBBF11692, 0xFACCF74A,
306 0x3FD40000, 0x000B6423,
307 0x3FD4AD71, 0xED5E62C2, 0xBBAD5590, 0x3A1538DA,
308 0x4008C2DD, 0x5F133FBD, 0xBBDDCA83, 0x1B096361,
309 0x3FD50000, 0x0017ECB5,
310 0x3FD5C9A0, 0x105DB3D0, 0x3BCE9A3C, 0xD13B3D17,
311 0x40077FE6, 0x3A7531D9, 0x3BE036AE, 0x48AA7EC4,
312 0x3FD60000, 0x00785229,
313 0x3FD6E8D8, 0x5AEC50EC, 0xBBC83B79, 0x6A9E4244,
314 0x40065948, 0x26E6F5E4, 0xBBFE36C3, 0x52F1803B,
315 0x3FD70000, 0x001C8263,
316 0x3FD80B4B, 0xD8D44659, 0xBBB499F2, 0xDC83C2D8,
317 0x40054B4F, 0x852854B9, 0x3BFB7522, 0xCB8BEDBE,
318 0x3FD80000, 0x0092A4EC,
319 0x3FD9312D, 0x86455645, 0x3B93DC76, 0x9413BABA,
320 0x400452E6, 0x95D5E9C0, 0xBBDBF037, 0x35429FA9,
321 0x3FD90000, 0x00353578,
322 0x3FDA5AB2, 0x70323F27, 0xBBD5E11A, 0x4F5FE60A,
323 0x40036D75, 0xEB84FD4A, 0xBBE3A384, 0xFD52C107,
324 0x3FDA0000, 0x002177A0,
325 0x3FDB8811, 0xE4F7B2DE, 0x3BC5B70A, 0x7D15AD3F,
326 0x400298CC, 0x1A74BDDF, 0x3BD41900, 0x83E8A76D,
327 0x3FDB0000, 0x0071F1F6,
328 0x3FDCB985, 0x9CFB26A0, 0x3BC8BA6E, 0xF094701B,
329 0x4001D30A, 0xD7B1694B, 0xBBFECEAB, 0x555C2BC9,
330 0x3FDC0000, 0x0071BD74,
331 0x3FDDEF49, 0xEB35D726, 0xBBCCC6A5, 0xDDB5912C,
332 0x40011A98, 0x2075779C, 0xBBF95852, 0x2DBF4C5B,
333 0x3FDD0000, 0x00A18446,
334 0x3FDF299D, 0xF3CB9E1D, 0xBBDAB47B, 0x6F4E4DBE,
335 0x40006E12, 0x72BBF45D, 0x3BE9A288, 0x2C8CD3D9,
336 0x3FDE0000, 0x0009F5E0,
337 0x3FE03461, 0xF090AA3D, 0xBBB6411F, 0x2E55E33D,
338 0x3FFF988E, 0xC8064412, 0xBBE51DE0, 0x7EC2ADBC,
339 0x3FDF0000, 0x00263EAE,
340 0x3FE0D680, 0x92D62104, 0x3BDD6620, 0x8142D974,
341 0x3FFE6858, 0x0C88DF44, 0xBBDD8335, 0xD4838950,
342 0x3FE00000, 0x002FD5FE,
343 0x3FE17B4F, 0x5C31640E, 0xBBCF7FD9, 0x4387A68E,
344 0x3FFD49AD, 0x7DDFB168, 0x3BF71CBF, 0x13808D41,
345 0x3FE08000, 0x00760E88,
346 0x3FE222F4, 0xAFFFC590, 0x3BDD3CE8, 0x22236C45,
347 0x3FFC3AF4, 0x70DCC003, 0xBBFBC120, 0xA2F987AC,
348 0x3FE10000, 0x000F8806,
349 0x3FE2CD98, 0xFEB5905C, 0x3BE169C8, 0x3A579A2C,
350 0x3FFB3AC2, 0x7605E88C, 0xBBC11C38, 0xB96450D3,
351 0x3FE18000, 0x002B5E5D,
352 0x3FE37B66, 0xF43D000E, 0x3BD2830E, 0xBAC9646C,
353 0x3FFA47D6, 0x6F7CCB9D, 0xBBEB41B3, 0x171B37E1,
354 0x3FE20000, 0x000A3C10,
355 0x3FE42C8B, 0xA0F7A0E3, 0x3BE38DF3, 0x66A12B85,
356 0x3FF96112, 0xDACBFCBC, 0xBBEFEF64, 0xEF1EA7D1,
357 0x3FE28000, 0x00595422,
358 0x3FE4E136, 0xB0CFA783, 0xBBBA398C, 0xEE061521,
359 0x3FF88578, 0xFAA0A3DE, 0x3BDA1845, 0x4AE391EC,
360 0x3FE30000, 0x0026966D,
361 0x3FE5999A, 0x9E477C73, 0xBBCE4AC5, 0x97E84259,
362 0x3FF7B424, 0xCEF83049, 0xBBF35B98, 0x9FE75083,
363 0x3FE38000, 0x0015B3FD,
364 0x3FE655EC, 0xF397B5BE, 0xBBE51A33, 0x112B4404,
365 0x3FF6EC49, 0xA1E2DBE5, 0x3BB86EE4, 0x93E45292,
366 0x3FE40000, 0x00146EE3,
367 0x3FE71666, 0x89F330B5, 0xBBE35378, 0x5A4867EE,
368 0x3FF62D2F, 0x215A3E9B, 0x3BEA0C9E, 0xF65CC2CF,
369 0x3FE48000, 0x000F229A,
370 0x3FE7DB43, 0xD3A2EEDF, 0xBBE0A454, 0xC31E28C0,
371 0x3FF5762E, 0xE1176AE6, 0xBBE1C3BD, 0x493469EA,
372 0x3FE50000, 0x002E680F,
373 0x3FE8A4C5, 0x2CF1484A, 0xBB9432B2, 0x063946DC,
374 0x3FF4C6B2, 0x384B4991, 0x3BD8E9D0, 0xF949DFE6,
375 0x3FE58000, 0x00672096,
376 0x3FE9732F, 0x3459A187, 0x3BDDADD7, 0x84F3E87C,
377 0x3FF41E30, 0x6C72DBD4, 0xBBEC67AE, 0xAE1EBB1F,
378 0x3FE60000, 0x00B6C23F,
379 0x3FEA46CB, 0x2D189DAF, 0x3BB739F4, 0x75A8271C,
380 0x3FF37C2D, 0x1B3BE6EB, 0xBBC85203, 0x8CECB4C0,
381 0x3FE68000, 0x0028A725,
382 0x3FEB1FE7, 0x6A3CF1FF, 0xBBEA7D29, 0x7EE6B9D1,
383 0x3FF2E036, 0xDC04F43A, 0x3BED69F5, 0xAF37688F,
384 0x3FE70000, 0x004CA480,
385 0x3FEBFED7, 0xCD2DB8E2, 0xBBCD123B, 0x0537EE2D,
386 0x3FF249E6, 0x09D6492E, 0xBBD38032, 0x80A8D415,
387 0x3FE78000, 0x000FFDA8,
388 0x3FECE3F6, 0x42FE6158, 0x3BD569F8, 0x793799F7,
389 0x3FF1B8DB, 0xBED5846E, 0x3BD9797E, 0x39EAD48F,
390 0x3FE80000, 0x001E1B0A,
391 0x3FEDCFA3, 0x61492AAF, 0xBBDA3DEC, 0xB61278D5,
392 0x3FF12CC0, 0xE58AA7DC, 0x3BB3088E, 0x1168567F,
393 0x3FE88000, 0x000B208A,
394 0x3FEEC247, 0x07D4CEEF, 0xBBE3FF76, 0x1F807DFF,
395 0x3FF0A545, 0x6F3C80B0, 0xBBCF3891, 0xE44E8437,
396 0x3FE90000, 0x000657F8,
397 0x3FEFBC51, 0x1E0226B6, 0xBBEC6FF5, 0x5F02EDF6,
398 0x3FF0221F, 0x9DADF743, 0xBBD42807, 0xFD02BF0B,
399 0x3FE98000, 0x00655AE2,
400 0x3FF05F1D, 0x31752EE7, 0xBBD8BE76, 0xB0E002E9,
401 0x3FEF4616, 0xC9018897, 0xBBDB4E55, 0xDEAF57E4
402};
403
404#warning These probably ought to static or typedefs below. They are being exported.
405
406struct TrigTblEntry {
407 double arg;
408 double sinHead;
409 double sinTail;
410 double cosHead;
411 double cosTail;
412} TrigTblEntry;
413
414struct CosCoeffEntry {
415 double Head;
416 double Tail;
417} CosCoeffEntry;
418
419struct SinCoeffEntry {
420 double Head;
421 double Tail;
422} SinCoeffEntry;
423
424struct CosCoeff16Entry {
425 double Head;
426 double Tail;
427} CosCoeff16Entry;
428
429struct SinCoeff16Entry {
430 double Head;
431 double Tail;
432} SinCoeff16Entry;
433
434struct TanTblEntry {
435 double arg;
436 double TanHead;
437 double TanTail;
438 double CotHead;
439 double CotTail;
440} TanTblEntry;
441
442struct TanCoeffEntry {
443 double Head;
444 double Tail;
445} TanCoeffEntry;
446
447struct CotCoeff16Entry {
448 double Head;
449 double Tail;
450} CotCoeff16Entry;
451
452struct TanCoeff16Entry {
453 double Head;
454 double Tail;
455} TanCoeff16Entry;
456
457static void argReduce( double xHead, double xTail, double *dLeft, int *M, double *redHead, double *redTail )
458{
459 double amid, alow, carry, high, q, q1;
460 double b, blow, b1, c, clow, c1, c2, c3, c4;
461 double d, dlow, d1, d2, d3, d4, d5, e1, e2, e3;
462 Double DeN;
463 amid = 0.0;
464 alow = 0.0;
465
466 while (fabs(xHead) > k2To45m1) // reduce until argument
467 { // isn't larger than k2To45m1
468 q = xHead*k2ByPi;
469 if (fabs(q) < k2To54) { // maybe q(mod 4) != 0
470 if (fabs(q) < k2To52) { // maybe q isn't even an integer
471 if (q > 0)
472 q = __FADD( q, k2To52 ) - k2To52; // remove fraction if any exists
473 else
474 q = __FSUB( q, k2To52 ) + k2To52;
475 } // at this point, q is an integer
476 q1 = q*0.25; // and determine q(mod 4)
477 if (q > 0)
478 q1 = __FADD( q1, k2To52 ) - k2To52; // giving the quadrant in which the angle resides.
479 else
480 q1 = __FSUB( q1, k2To52 ) + k2To52; // right adjust quadrant
481 DeN.f = __FNMSUB( 4.0, q1, q ) + k1_8To52; // q - 4.0*q1 + k1_8To52;
482
483 *M = *M + DeN.i[1];
484 }
485
486 //****************************************************************************
487 // Tackling the subtraction of q*(pi/2) from the octuple precision word: *
488 // (xHead, xTail, amid, alow) *
489 //****************************************************************************
490
491 high = __FNMSUB( q, kPiBy2, xHead ); // xHead - q*kPiBy2; // exact (Spectacular use of MAF)
492 b = q*kPiBy2Tail1;
493 blow = __FNMSUB( q, kPiBy2Tail1, b ); // b - q*kPiBy2Tail1; // Yet another use of MAF
494 c = q*kPiBy2Tail2;
495 clow = __FNMSUB( q, kPiBy2Tail2, c ); // c - q*kPiBy2Tail2;
496 d = q*kPiBy2Tail3;
497 dlow = __FNMSUB( q, kPiBy2Tail3, d ); // d - q*kPiBy2Tail3;
498 b1 = __FADD( xTail, high );
499 if (fabs(xTail) > fabs(high)) // sum and propagate carries to the right
500 c1 = xTail - b1 + high;
501 else
502 c1 = high - b1 + xTail;
503 xHead = __FSUB( b1, b );
504 if (fabs(b1) > fabs(b))
505 c2 = b1 - xHead - b + c1;
506 else
507 c2 = b1 - (xHead+b) + c1;
508 c3 = __FADD( amid, blow);
509 if (fabs(amid) > fabs(blow))
510 d1 = amid - c3 + blow;
511 else
512 d1 = blow - c3 + amid;
513 c4 = __FSUB( c2, c );
514 if (fabs(c2) > fabs(c))
515 d2 = c2 - c4 - c + d1;
516 else
517 d2 = c2 - __FADD( c4, c ) + d1;
518 xTail = __FADD( c3, c4 );
519 if (fabs(c3) > fabs(c4))
520 d3 = c3 - xTail + c4 + d2;
521 else
522 d3 = c4 - xTail + c3 + d2;
523 d4 = __FADD( alow, clow );
524 if (fabs(alow) > fabs(clow))
525 e1 = alow - d4 + clow;
526 else
527 e1 = clow - d4 + alow;
528 d5 = __FSUB( d3, d );
529 if (fabs(d3) > fabs(d))
530 e2 = d3 - d5 - d + e1;
531 else
532 e2 = d3 - (d5+d) + e1;
533 amid = __FADD( d4, d5 );
534 if (fabs(d4) > fabs(d5))
535 e3 = d4 - amid + d5 + e2;
536 else
537 e3 = d5 - amid + d4 + e2;
538 alow = e3 + dlow;
539 }
540
541 high = __FADD( xHead, xTail ); // At long last,
542 carry = xHead - high + xTail; // A final distillation
543
544 xTail = __FADD( carry, __FADD( amid, alow ) );
545 xHead = high;
546 *dLeft = __FSUB( carry, xTail ) + __FADD( amid, alow );
547
548 *redHead = xHead;
549 *redTail = xTail;
550}
551
552static void SinCosCore(int M, double xHead, double xTail, double dleft, double *yHead, double *yTail )
553{
554 double intquo;
555 double temp, temps, tempc;
556 double b, c, ca, d;
557 double redarg, redarg1, res, reslow, resmid, resbot, restop, resin, resextra;
558 double ressup, restiny, resinf;
559 double p1, p2, p3, p4, p5, p6, p6a, p7, p8, p9, p10;
560 double p11, p12, p13, p14, p15, p16, p17, p18, p19, p20;
561 double p21, p22, p23, p23a, p23b, p24, p25, p26, p27, p27a, p28, p29, p30, p31;
562 double arg, arg2, arg2a, argsq, argsq2, argmid, argnew, argres;
563 double sin, sin2, sint, sinl, sinargtiny;
564 double sintab, sintab2, arg2b;
565 double cos, cos2, cosl, cost;
566 double prods, prodsl, prodc, prodcl;
567 double blow, carry, cerr;
568 int i, index;
569 Double DeN;
570
571 struct TrigTblEntry *TrigPtr = (struct TrigTblEntry *)TrigTable - 6;
572 struct CosCoeffEntry *CosPtr = (struct CosCoeffEntry *) CosCoeff;
573 struct SinCoeffEntry *SinPtr = (struct SinCoeffEntry *) SinCoeff;
574 struct CosCoeff16Entry *Cos16Ptr = (struct CosCoeff16Entry *) CosCoeff16;
575 struct SinCoeff16Entry *Sin16Ptr = (struct SinCoeff16Entry *) SinCoeff16;
576
577 DeN.f = __FMADD( xHead, k2ByPi, k1_8To52 ); // xHead*k2ByPi + k1_8To52; // # of multiples of pi/2 to remove
578
579 intquo = DeN.f - k1_8To52; // Integer
580 redarg = __FNMSUB( intquo, kPiBy2, xHead ); // xHead - intquo*kPiBy2; // Exact result
581
582 M = M + DeN.i[1]; // M=quadrant in which to comp. sin
583
584 b = kPiBy2Tail1*intquo; // continue argument reduction
585 blow = __FMSUB( kPiBy2Tail1, intquo, b ); // kPiBy2Tail1*intquo - b;
586 redarg1 = __FSUB( xTail, b ); // second word of reduced argument
587 DeN.f = __FMADD( fabs(redarg), 64.0, k1_8To52 );// fabs(redarg)*64.0 + k1_8To52; // index into acctbl
588 index = DeN.i[1];
589
590 ca = kPiBy2Tail2*intquo;
591
592 if (fabs(b) > fabs(xTail))
593 carry = xTail - (b+redarg1); // tail-b may be inexact
594 else
595 carry = (xTail-redarg1) - b;
596
597 if (! (index < 6)) {
598 if (redarg > 0.0) {
599 redarg = redarg - TrigPtr[index].arg ; // arg-intquo ln2-tbl(index) ->
600 sintab = TrigPtr[index].sinHead ; // fetch sine of table value
601 sintab2 = TrigPtr[index].sinTail ;
602 }
603 else { // much trailing zeroes in redarg
604 redarg = redarg + TrigPtr[index].arg;
605 sintab =- TrigPtr[index].sinHead ; // fetch sine of table value
606 sintab2 =- TrigPtr[index].sinTail ;
607 }
608 }
609
610 arg = __FADD( redarg, redarg1 ); // get signif. bits from the right
611 arg2a = (redarg - arg + redarg1);
612 arg2b = __FSUB( arg2a, blow ); // fill in with bits from 2nd term
613 c = __FSUB( carry, ca ); // 3rd word of reduced argument
614
615 argmid = __FADD( arg2b , c );
616 argnew = __FADD( arg, argmid );
617 argres = arg - argnew + argmid;
618 arg = argnew;
619 cerr = (carry-c) - ca;
620
621 if (fabs(arg2b) > fabs(c))
622 d = arg2b - argmid + c + dleft;
623 else
624 d = c - argmid + arg2b + dleft;
625
626 // d = d + cerr + (ca-kPiBy2Tail2*intquo) - kPiBy2Tail3*intquo; // 4th word
627 d = __FNMSUB( kPiBy2Tail3, intquo, d + cerr + __FNMSUB( kPiBy2Tail2, intquo, ca ) );
628 arg2 = __FADD( argres, d ); // reduced arg is (arg, arg2)
629 temp = 2.0*arg*arg2;
630 argsq = __FMADD( arg, arg, temp ); // arg*arg + temp;
631 argsq2 = __FMADD( arg2, arg2, __FMSUB( arg, arg, argsq ) + temp ); // arg*arg - argsq + temp + arg2*arg2;
632
633
634 // Now we have the reduced arg = (arg, arg2)
635
636 if (index < 6) { // straight power series
637
638 if ((M & 1) == 0) { // sin series evaluation
639 sin = Sin16Ptr[6].Head ; // need extra terms when
640 sinl = Sin16Ptr[6].Tail ; // used without table lookup
641 sinargtiny = argres - arg2 + d;
642
643 for(i = 5; i > -1; i--) {
644 temps = __FMADD( sin, argsq2, sinl*argsq ); // sin*argsq2 + sinl*argsq;
645 prods = __FMADD( sin, argsq, temps ); // sin*argsq + temps;
646 prodsl = __FMSUB( sin, argsq, prods ) + temps; // sin*argsq - prods + temps;
647 sin = __FADD( Sin16Ptr[i].Head, prods );
648 sin2 = Sin16Ptr[i].Tail + prodsl;
649 sinl = Sin16Ptr[i].Head - sin+prods + sin2;
650 }
651 }
652 else { // cos series evaluation
653 cos = Cos16Ptr[6].Head;
654 cosl = Cos16Ptr[6].Tail ;
655
656 for(i = 5; i > -1; i--) {
657 temps = __FMADD( cos, argsq2, cosl*argsq ); // cos*argsq2 + cosl*argsq;
658 prods = __FMADD( cos, argsq, temps ); // cos*argsq + temps;
659 prodsl = __FMSUB( cos, argsq, prods) + temps; // cos*argsq - prods + temps;
660 cos = __FADD( Cos16Ptr[i].Head, prods );
661 cos2 = Cos16Ptr[i].Tail + prodsl;
662 cosl = Cos16Ptr[i].Head - cos + prods + cos2;
663 }
664 }
665 }
666 else { // shorter series in case
667
668 cos = CosPtr[4].Head; // where exact table is used.
669 sin = SinPtr[4].Head;
670 cosl = CosPtr[4].Tail;
671 sinl = SinPtr[4].Tail;
672 for (i = 3; i > -1; i--) { // add short series
673 tempc = __FMADD( cos, argsq2, cosl*argsq ); // cos*argsq2 + cosl*argsq;
674 prodc = __FMADD( cos, argsq, tempc ); // cos*argsq + tempc;
675 prodcl = __FMSUB( cos, argsq, prodc) + tempc; // cos*argsq - prodc + tempc;
676 cos = __FADD( CosPtr[i].Head, prodc );
677 cos2 = CosPtr[i].Tail + prodcl;
678 cosl = CosPtr[i].Head - cos + prodc + cos2;
679 temps = __FMADD( sin, argsq2, sinl*argsq ); // sin*argsq2 + sinl*argsq;
680 prods = __FMADD( sin, argsq, temps ); // sin*argsq + temps;
681 prodsl = __FMSUB( sin, argsq, prods ) + temps; // sin*argsq - prods + temps;
682 sin = __FADD( SinPtr[i].Head, prods );
683 sin2 = SinPtr[i].Tail + prodsl;
684 sinl = SinPtr[i].Head - sin + prods + sin2;
685 }
686 }
687
688 tempc = __FMADD( cos, argsq2, cosl*argsq ); // cos*argsq2 + cosl*argsq;
689 cost = __FMADD( cos, argsq, tempc ); // cos*argsq + tempc;
690 cosl = __FMSUB( cos, argsq, cost ) + tempc; // cos*argsq - cost + tempc;
691 cos = cost; // cos(arg) - 1.0 complete
692 temps = __FMADD( sin, argsq2, sinl*argsq ); // sin*argsq2 + sinl*argsq;
693 sint = __FMADD( sin, argsq, temps ); // sin*argsq + temps;
694 sinl = __FMSUB( sin, argsq, sint ) + temps; // sin*argsq - sint + temps;
695 temps = __FMADD( sint, arg2, sinl * arg ); // sint*arg2 + sinl*arg;
696 sin = __FMADD( sint, arg, temps ); // sint*arg + temps;
697 sinl = __FMSUB( sint, arg, sin ) + temps; // sint*arg - sin + temps; // sin(arg) - arg complete
698
699 if (index < 6) {
700 if ((M & 1) == 0) {
701 res = __FADD( arg, sin ); // sin of a small angle
702 reslow = __FADD( arg2, sinl ); // careful distallation
703 resmid = __FADD( arg - res, sin );
704 resbot = arg2 - reslow + sinl;
705 restop = __FADD( res, (reslow+resmid) );
706 resin = __FADD( res - restop, (reslow+resmid) );
707 resextra = resmid - (reslow+resmid) + reslow;
708 ressup = __FADD( restop, (resin+resbot) );
709 restiny = __FSUB( resin, (resin+resbot) + resbot );
710 resinf = restop - ressup + (resin+resbot) + (resextra+restiny+sinargtiny);
711 *yHead = ressup;
712 *yTail = resinf;
713 }
714 else {
715 res = __FADD( 1.0, cos ); // cos of a small angle
716 reslow = cosl;
717 resmid = __FADD( 1.0 - res, cos );
718 restop = __FADD( res, (reslow+resmid) );
719 resin = __FADD( res - restop, (reslow+resmid) );
720 resextra = resmid - (reslow+resmid) + reslow;
721 ressup = __FADD( restop, resin );
722 resinf = restop - ressup + resin + resextra;
723 *yHead = ressup;
724 *yTail = resinf;
725 }
726 } // end of if (index < 6)
727 else {
728 if ((M & 1) == 0) { // even quadrant-- eval. sine(x)
729 p1 = cos*sintab; // o(-8) exact
730 p2 = __FMSUB( cos, sintab, p1 ); // cos*sintab - p1; // o(-62) exact
731 p3 = __FMADD( cosl, sintab, cos*sintab2 ); //cosl*sintab + (cos*sintab2); // o(-67)
732 p4 = arg* TrigPtr[index].cosHead ; // o(-4) exact
733 p5 = __FMSUB(arg, TrigPtr[index].cosHead, p4 ); // arg* TrigPtr[index].cosHead -p4; // o(-58) exact
734 p6 = arg2* TrigPtr[index].cosHead ; // o(-58)
735 p6a = __FMSUB( arg2, TrigPtr[index].cosHead, p6 ); // arg2* TrigPtr[index].cosHead -p6;
736 p7 = sin* TrigPtr[index].cosHead ; // o(-4) exact
737 p8 = __FMSUB( sin, TrigPtr[index].cosHead, p7 ); // sin* TrigPtr[index].cosHead -p7; // o(-59) exact
738 p9 = __FMADD( sinl, TrigPtr[index].cosHead, (sin* TrigPtr[index].cosTail ) ); //sinl* TrigPtr[index].cosHead +(sin* TrigPtr[index].cosTail );
739 p10 = __FMADD( arg, TrigPtr[index].cosTail, p9 + p8); // arg* TrigPtr[index].cosTail +p9+p8; // o(-59)
740 p11 = __FADD( p2, p5 ); // o(-58) exact
741 p12 = p5 - p11 + p2 + p6a; // o(-112) exact
742 p13 = __FADD( p6, p11 ); // o(-57) exact
743 p15 = __FADD( sintab, p4 ); // o(0) exact
744 p16 = sintab - p15 + p4; // o(-54) exact
745 p17 = sintab2 + p10 + p12 + p3; // o(-58)
746 p18 = __FADD( p15, p1 ); // o(0) exact
747 p20 = __FADD( p18, p7 ); // o(0) exact
748 if (fabs(p6) > fabs(p11))
749 p14 = p6 - p13 + p11; // o(-111)
750 else
751 p14 = p11 - p13 + p6;
752 p19 = p15 - p18 + p1 + p16; // o(-53) exact
753 p21 = p18 - p20 + p7; // o(-54) exact
754 p22 = p13 + p17; // o(-56)
755 p23 = __FADD( p19, p21 ); // o(-52) exact
756 p25 = __FADD( p20, p23 ); // o(0) exact
757 p26 = p20 - p25 + p23; // o(-54) exact
758 if (fabs(p19) > fabs(p21))
759 p24 = p19 - p23 + p21; // o(-108)
760 else
761 p24 = p21 - p23 + p19;
762 p27 = __FADD( (p14+p24), p22 + p26 ); // o(-54)
763 p27a = p26 - p27 + ((p14+p24)+p22);
764 p28 = __FADD( p25, p27 );
765
766 *yHead = p28;
767 *yTail = p25 - p28 + p27 + p27a;
768 }
769 // end of if ((M & 1) == 0)
770 else { // evaluation of cos series.
771 p1 = TrigPtr[index].cosHead *cos; // o(-16) exact
772 p2 = __FMSUB( TrigPtr[index].cosHead, cos, p1 ); // TrigPtr[index].cosHead *cos-p1; // o(-71) exact
773 p5 = sintab*arg; // o(-8) exact
774 p6 = __FMSUB( sintab, arg, p5 ); // sintab*arg - p5 ; // o(-62) exact
775 p7 = sintab*arg2; // o(-62)
776 p8 = sintab*sin; // o(-24) exact
777 p9 = __FMSUB( sintab, sin, p8 ); // sintab*sin - p8 ; // o(-78) exact
778 p12 = __FMSUB( TrigPtr[index].cosTail, cos, (sintab2*arg) ); // TrigPtr[index].cosTail *cos - (sintab2*arg); // o(-80)
779 p13 = __FNMSUB( sintab2, sin, p12 ); // p12 - (sintab2*sin); // o(-72)
780 p14 = p13 + p2; // o(-70)
781 p15 = __FMADD( TrigPtr[index].cosHead, cosl, p14 ); // p14 + TrigPtr[index].cosHead *cosl; // o(-69)
782 p16 = __FNMSUB( sintab, sinl, p15 - p9 );// p15 - p9 - sintab*sinl; // o(-67)
783 p17 = __FADD( TrigPtr[index].cosHead, p1 ); // o(-1) exact
784 p18 = TrigPtr[index].cosHead - p17 + p1; // o(-55) exact-16 bits
785 p19 = __FSUB( p17, p5 ); // o(-1) exact
786 p20 = p17 - p19 - p5; // o(-55) exact-8 bits
787 p21 = p18 + p20 ; // o(-54) exact-17 bits
788 p22 = __FSUB( p19, p8 ); // o(-1) exact
789 p23 = p19 - p22 - p8; // o(-53) exact-25 bits
790 p24 = TrigPtr[index].cosTail + p16; // o(-65)
791 p23a = __FADD( p23, p21 ); // o(-54)
792 p25 = __FSUB( p23a, p6 ); // o(-53) exact
793 p26 = p23a - p25 - p6; // o(-106)
794 p27 = __FSUB( p25, p7 ); // o(-53) exact
795 if (fabs(p23) > fabs(p21))
796 p23b = p23 - p23a + p21; // o(-106)
797 else
798 p23b = p21 - p23a + p23;
799 p28 = p25 - p27 - p7 + p26 + p23b; // o(-105)
800 p29 = __FADD( p27, p24 ); // o(-53) exact
801 p30 = p27 - p29 + p24 + p28; // o(-104)
802 p31 = __FADD( p22, p29 ); // o(-1) exact
803
804 *yHead = p31;
805 *yTail = p22 - p31 + p29 + p30;
806 } //end of else if ((M & 1) == 0)
807
808 } //end of else if (index < 6)
809
810
811 if((M & 2) != 0) { // change sign for below x axis.
812 *yHead = - *yHead;
813 *yTail = - *yTail;
814 }
815}
816
817long double sinl( long double x )
818{
819 uint32_t hibits;
820 double xHead, xTail, dleft, fenv, yHead, yTail;
821 DblDbl ddx;
822 int M;
823
824 FEGETENVD(fenv);
825 FESETENVD(0.0);
826
827 ddx.f = x;
828
829 xHead = ddx.d[0];
830 xTail = ddx.d[1];
831
832 hibits = ddx.i[0] & 0x7FFFFFFFu; // highest 32 bits as uint32_t
833
834 // NaNs and infinities
835
836 if (hibits >= 0x7ff00000u) { // special cases: NaN and infinity
837 if (xHead != xHead) { // x is a NaN
838 FESETENVD(fenv);
839 return x;
840 }
841 else {
842 ddx.d[0] = xHead - xHead; // inf - inf gives NaN
843 ddx.d[1] = ddx.d[0];
844 FESETENVD(fenv);
845 return (ddx.f);
846 }
847 }
848
849 // x = zero hence xHead is zero
850 if ((hibits | ddx.i[1]) == 0) {
851 FESETENVD(fenv);
852 return x;
853 }
854
855 // x is not infinity, NaN, or zero
856
857 //for tiny x, sin x = x
858 if (hibits < 0x29900000u) { // 2^(-358); 358 = (1022 + 52 )/3
859 FESETENVD(fenv);
860 return x;
861 }
862
863 M = 0; // M = 0 for Sine, M = 1 for Cosine
864 dleft = 0.0;
865 if (hibits >= 0x42d00000u) { // |x| > = 2 ^ 46
866 argReduce(xHead, xTail, &dleft, &M, &xHead, &xTail);
867 }
868
869 // Call core algorithm
870 SinCosCore( M, xHead, xTail, dleft, &yHead, &yTail );
871
872 ddx.d[0] = yHead + yTail; // Cannonize result
873 ddx.d[1] = (yHead - ddx.d[0]) + yTail;
874
875 FESETENVD(fenv);
876 return (ddx.f);
877}
878
879long double cosl( long double x )
880{
881 uint32_t hibits;
882 double xHead, xTail, dleft, fenv, yHead, yTail;
883 DblDbl ddx;
884 int M;
885
886 FEGETENVD(fenv);
887 FESETENVD(0.0);
888
889 ddx.f = x;
890
891 xHead = ddx.d[0];
892 xTail = ddx.d[1];
893
894 hibits = ddx.i[0] & 0x7FFFFFFFu; // highest 32 bits as uint32_t
895
896 // NaNs and infinities
897
898 if (hibits >= 0x7ff00000u) { // special cases: NaN and infinity
899 if (xHead != xHead) { // x is a NaN
900 FESETENVD(fenv);
901 return x;
902 }
903 else {
904 ddx.d[0] = xHead - xHead; // inf - inf gives NaN
905 ddx.d[1] = ddx.d[0];
906 FESETENVD(fenv);
907 return (ddx.f);
908 }
909 }
910
911 // x = zero hence xHead is zero
912 if ((hibits | ddx.i[1]) == 0) {
913 FESETENVD(fenv);
914 return (long double) 1.0;
915 }
916
917 // x is not infinity, NaN, or zero
918
919 //for tiny x, cos x = 1.0
920 if (hibits < 0x1e600000u) { // 2^(-537); 537 = (1022 + 52 )/2
921 FESETENVD(fenv);
922 return 1.0;
923 }
924
925 M = 1; // M = 0 for Sine, M = 1 for Cosine
926 dleft = 0.0;
927 if (hibits >= 0x42d00000u) { // |x| > = 2 ^ 46
928 argReduce(xHead, xTail, &dleft, &M, &xHead, &xTail);
929 }
930
931 // Call core algorithm
932 SinCosCore( M, xHead, xTail, dleft, &yHead, &yTail );
933
934 ddx.d[0] = __FADD( yHead, yTail ); // Cannonize result
935 ddx.d[1] = (yHead - ddx.d[0]) + yTail;
936
937 FESETENVD(fenv);
938 return (ddx.f);
939}
940
941long double tanl( long double x )
942{
943 double temp, temps;
944 double b, c, c1, c2, ca, cerr, d;
945 double frac, frac2, frac3, t1, t2;
946 double rec, res, reslow, resmid, resbot, restop, resin, resextra;
947 double ressup, restiny, resinf;
948 double rarg, rarg0, rarg2, rarg3, rarg2a, redarg, redarg1, redargold;
949 double p1, p2, p3, p4, p5, p6, p7, p8, p9, p10;
950 double p11, p12, p13, p14, p15, p16, p17, p18, p19, p20;
951 double p21, p22, p23, p24, p25, p26, p27, p28, p29, p30;
952 double p31, p32, p33, p34, p35, p36, p37, p38, p39, p40;
953 double p41, p42, p43, p44, p45, p46;
954 double arg, arg2, arg2a, arg2b, argsq, argsq2, argres, argmid, argnew;
955 double cot, cot2, cotl, cott, tan, tan2, tant, tanl;
956 double cotargtiny, tanargtiny;
957 double prods, prodsl, q;
958 double blow, carry, dleft;
959 int i, index, M;
960
961 uint32_t hibits;
962 Double DeN;
963 double fenv;
964
965 double xHead, xTail;
966 DblDbl ddx;
967
968 struct TanTblEntry *TanTablePtr = (struct TanTblEntry *)TanTable - 6;
969 struct TanCoeffEntry *TanPtr = (struct TanCoeffEntry *) TanCoeff - 1;
970 struct TanCoeff16Entry *Tan16Ptr = (struct TanCoeff16Entry *) TanCoeff16 - 1;
971 struct CotCoeff16Entry *Cot16Ptr = (struct CotCoeff16Entry *) CotCoeff16 - 1;
972
973 FEGETENVD(fenv);
974 FESETENVD(0.0);
975
976 ddx.f = x;
977
978 xHead = ddx.d[0];
979 xTail = ddx.d[1];
980
981 hibits = ddx.i[0] & 0x7FFFFFFFu; // highest 32 bits as uint32_t
982
983 if (hibits >= 0x7ff00000u){ // special cases: NaN and infinity
984 if (xHead != xHead) { // x is a NaN
985 FESETENVD(fenv);
986 return x;
987 }
988 else {
989 ddx.d[0] = xHead - xHead; // inf - inf gives NaN & raises invalid
990 ddx.d[1] = ddx.d[0];
991 FESETENVD(fenv);
992 return (ddx.f);
993 }
994 }
995
996 // x is not infinity nor a NaN
997
998 // x = zero
999 if ((hibits | ddx.i[1]) == 0) {
1000 FESETENVD(fenv);
1001 return x;
1002 }
1003
1004
1005 // finite x
1006
1007 //for tiny x, tan x = x
1008 if (hibits < 0x29900000u) { // 2^(-358); 358 = (1022 + 52 )/3
1009 FESETENVD(fenv);
1010 return x;
1011 }
1012
1013 // finite result but not tiny
1014
1015 M=0; // determine quadrant
1016 dleft=0.0;
1017
1018 if (hibits >= 0x42b00000u){ // |x| > = 2 ^ 44
1019 argReduce (xHead, xTail, &dleft, &M, &xHead, &xTail);
1020 }
1021
1022 // now finite result < 2 ^ 44
1023
1024 DeN.f = __FMADD( xHead, k2ByPi, k1_8To52 ); // xHead*k2ByPi + k1_8To52; // # of multiples of pi/2 to remove
1025 q = DeN.f - k1_8To52; // to remove integer
1026 redarg = __FNMSUB( q, kPiBy2, xHead ); // xHead - intquo*kPiBy2; // Exact result
1027
1028 M = M + DeN.i[1]; // M=quadrant in which to compute sin
1029 b = kPiBy2Tail1*q; // more argument reduction
1030 blow = __FMSUB( kPiBy2Tail1, q, b ); // kPiBy2Tail1*q - b;
1031 redarg1 = __FSUB( xTail, b ); // second word of reduced argument
1032 DeN.f = __FMADD( fabs(redarg), 64.0, k1_8To52 );// fabs(redarg)*64.0 + k1_8To52; // index into acctbl
1033 ca = kPiBy2Tail2*q;
1034
1035 if (fabs(b) > fabs(xTail))
1036 carry = xTail - (b+redarg1); // tail-b may be inexact
1037 else
1038 carry = (xTail-redarg1) - b;
1039
1040 redargold = redarg;
1041 index = DeN.i[1];
1042
1043 if (!(index < 6)) {
1044 if (redargold > 0.0)
1045 redarg = redarg - TanTablePtr[index].arg; // arg-q ln2-tbl(index) ->
1046 // mucho trailing zeroes in redarg
1047 else
1048 redarg = redarg + TanTablePtr[index].arg;
1049 }
1050
1051 arg = __FADD( redarg, redarg1 ); // get signif. bits from the right
1052 arg2a = (redarg - arg + redarg1);
1053 arg2b = __FSUB( arg2a, blow ); // fill in with bits from 2nd term
1054 c = __FSUB( carry, ca ); // 3rd word of reduced argument
1055
1056 argmid = __FADD( arg2b , c );
1057 argnew = __FADD( arg, argmid );
1058 argres = arg - argnew + argmid;
1059 arg = argnew;
1060 cerr = (carry-c) - ca;
1061
1062 if (fabs(arg2b) > fabs(c))
1063 d = arg2b - argmid + c + dleft;
1064 else
1065 d = c - argmid + arg2b + dleft;
1066
1067 // d = d + cerr + (ca-kPiBy2Tail2*q) - kPiBy2Tail3*q; // 4th word
1068 d = __FNMSUB( kPiBy2Tail3, q, d + cerr + __FNMSUB( kPiBy2Tail2, q, ca ) );
1069 arg2 = __FADD( argres, d ); // reduced arg is (arg, arg2)
1070 temp = 2.0*arg*arg2;
1071 argsq = __FMADD( arg, arg, temp ); // arg*arg + temp;
1072 argsq2 = __FMADD( arg2, arg2, __FMSUB( arg, arg, argsq ) + temp ); // arg*arg - argsq + temp + arg2*arg2;
1073
1074 if (index < 6) { // straight power series
1075 if ((M & 1) == 0) { // 2n multiple of pi/2 -> use TAN
1076 tan = Tan16Ptr[10].Head; // need extra terms when
1077 tanl = Tan16Ptr[10].Tail; // used without table lookup
1078 tanargtiny = argres - arg2 + d;
1079
1080 for(i = 9; i > 0; i--) {
1081 temps = __FMADD( tan, argsq2, tanl*argsq ); // tan*argsq2 + tanl*argsq;
1082 prods = __FMADD( tan, argsq, temps ); // tan*argsq + temps;
1083 prodsl = __FMSUB( tan, argsq, prods ) + temps; // tan*argsq - prods + temps;
1084 tan = __FADD( Tan16Ptr[i].Head, prods );
1085 tan2 = Tan16Ptr[i].Tail + prodsl;
1086 tanl = Tan16Ptr[i].Head - tan + prods + tan2;
1087 }
1088 temps = __FMADD( tan, argsq2, tanl*argsq ); // tan*argsq2 + tanl*argsq;
1089 tant = __FMADD( tan, argsq, temps ); // tan*argsq + temps;
1090 tanl = __FMSUB( tan, argsq, tant ) + temps; // tan*argsq - tant + temps;
1091 temps = __FMADD( tant, arg2, tanl*arg ); // tant*arg2 + tanl*arg;
1092 tan = __FMADD( tant, arg, temps ); // tant*arg + temps;
1093 tanl = __FMSUB( tant, arg, tan ) + temps; // tant*arg - tan + temps; // tan(arg)-arg complete
1094 res = __FADD( arg, tan ); // tan of a small angle
1095 reslow = __FADD( arg2, tanl ); // distill VERY CAREFULLY
1096 resmid = arg - res + tan;
1097 resbot = arg2 - reslow + tanl;
1098 restop = __FADD( res, (reslow+resmid) );
1099 resin = res - restop + (reslow+resmid);
1100 resextra = resmid - (reslow+resmid) + reslow;
1101 ressup = __FADD( restop, (resin+resbot) );
1102 restiny = resin - (resin+resbot) + resbot;
1103 resinf = restop - ressup + (resin+resbot) + ( resextra+restiny+tanargtiny);
1104
1105 //getting a cannonical result
1106
1107 ddx.d[0] = __FADD( ressup, resinf );
1108 ddx.d[1] = (ressup - ddx.d[0]) + resinf;
1109
1110 FESETENVD(fenv);
1111 return (ddx.f);
1112 }
1113 else { // near an 2k+1 multiple of pi/2
1114 // tan (x)=-cot(x+(2k+1) pi/2)
1115 cot = Cot16Ptr[9].Head; // need extra terms when
1116 cotl = Cot16Ptr[9].Tail; // used without table lookup
1117 cotargtiny = argres - arg2 + d;
1118 rarg0 = 1.0/arg; // We use the cotangent
1119
1120 for(i = 8; i > 0; i--){ // series, which is one term shorter
1121 temps = __FMADD( cot, argsq2, cotl*argsq ); // cot*argsq2 + cotl*argsq;
1122 prods = __FMADD( cot, argsq, temps ); // cot*argsq + temps;
1123 prodsl = __FMSUB( cot, argsq, prods ) + temps; // cot*argsq - prods + temps;
1124 cot = __FADD( Cot16Ptr[i].Head, prods );
1125 cot2 = Cot16Ptr[i].Tail + prodsl;
1126 cotl = Cot16Ptr[i].Head - cot + prods + cot2;
1127 }
1128 cott = cot;
1129 rarg = rarg0;
1130 temps = __FMADD( cott, arg2, cotl*arg ); // cott*arg2 + cotl*arg;
1131 cot = __FMADD( cott, arg, temps ); // cott*arg + temps;
1132 cotl = __FMSUB( cott, arg, cot ) + temps; // cott*arg - cot + temps; // 1/arg-cot(arg) completed
1133 p1 = __FNMSUB( arg, rarg, 1.0 ); // 1.0 - arg*rarg; // next complete hi precision
1134 p2 = arg2*rarg; // reciprocal of argument, to subtract
1135 p3 = __FSUB( p1, p2 ); // from series. The reciprocal will
1136 rarg2a = p3*rarg; // be the dominant term.
1137 rarg2 = rarg2a + __FNMSUB( arg, rarg2a, p3 )*rarg; // rarg2a + (p3-arg*rarg2a)*rarg
1138
1139 if (fabs(p1) > fabs(p2))
1140 p4 = p1 - p3 - p2;
1141 else
1142 p4 = p1 - (p3+p2);
1143
1144 p5 = p4 + __FNMSUB( arg2, rarg, p2 ); // p4 + (p2-arg2*rarg)
1145 p6 = __FNMSUB( cotargtiny, rarg, p5 ); // p5 - cotargtiny*rarg;
1146 p7 = __FNMSUB( rarg2, arg, p3 ); // p3 - rarg2*arg;
1147 rarg3 = __FNMSUB( arg2, rarg2, p6+p7 )*rarg; // ((p6+p7) - arg2*rarg2)*rarg; // low order recip. bits
1148
1149 //**********************************************************************
1150 // BEGIN the subtraction of 1/arg from the series *
1151 // |cot| much greater than |rarg| *
1152 //**********************************************************************
1153
1154 p8 = __FSUB( cot, rarg );
1155 p9 = cot - (rarg+p8);
1156 p10 = __FSUB( cotl, rarg2 );
1157 p11 = cotl - (rarg2+p10);
1158 p12 = p11 - rarg3;
1159 p13 = __FADD( p9, p10 );
1160
1161 if (fabs(p9) > fabs(p10))
1162 p14 = p9 - p13 + p10;
1163 else
1164 p14 = p10 - p13 + p9;
1165
1166 p15 = p12 + p14;
1167 p16 = __FADD( p8, p13 ); // HO term of result
1168 p17 = p8 - p16 + p13; // almost LO term
1169 p18 = p17 + p15; // assimilate any other LO bits
1170
1171 // getting a cannonical result
1172
1173 ddx.d[0] = __FADD( p16, p18 );
1174 ddx.d[1] = (p16 - ddx.d[0]) + p18;
1175
1176 FESETENVD(fenv);
1177 return (ddx.f);
1178 }
1179 }
1180 else { // Serious cases
1181 // Accurate tables will be used.
1182 tant = TanPtr[7].Head; // Need fewer terms when
1183 tant = __FMADD( tant, argsq, TanPtr[6].Head ); // TanPtr[6].Head + tant*argsq;
1184 tan = __FMADD( tant, argsq, TanPtr[5].Head ); // TanPtr[5].Head + tant*argsq;
1185 tanl = __FMADD( tant, argsq, TanPtr[5].Head - tan );// TanPtr[5].Head - tan + tant*argsq;
1186 tanargtiny = argres - arg2 + d;
1187
1188 for(i = 4; i > 0; i--) {
1189 temps = __FMADD( tan, argsq2, tanl*argsq ); // tan*argsq2 + tanl*argsq;
1190 prods = __FMADD( tan, argsq, temps ); // tan*argsq + temps;
1191 prodsl = __FMSUB( tan, argsq, prods ) + temps; // tan*argsq - prods + temps;
1192 tan = __FADD( TanPtr[i].Head, prods );
1193 tan2 = TanPtr[i].Tail + prodsl;
1194 tanl = TanPtr[i].Head -tan + prods + tan2;
1195 }
1196 tant = tan;
1197 temps = __FMADD( tant, arg2, tanl*arg ); // tant*arg2 + tanl*arg;
1198 tan = __FMADD( tant, arg, temps ); // tant*arg + temps;
1199 tanl = __FMSUB( tant, arg, tan ) + temps; // tant*arg - tan + temps; // tan(arg) complete
1200 if (redargold > 0) { // read out table values
1201 t1 = TanTablePtr[index].TanHead; // (t1,t2) are tan(x)
1202 t2 = TanTablePtr[index].TanTail;
1203 c1 = TanTablePtr[index].CotHead; // (c1,c2) are cot(x)
1204 c2 = TanTablePtr[index].CotTail;
1205 }
1206 else {
1207 t1 = -TanTablePtr[index].TanHead;
1208 t2 = -TanTablePtr[index].TanTail;
1209 c1 = -TanTablePtr[index].CotHead;
1210 c2 = -TanTablePtr[index].CotTail;
1211 }
1212
1213 p1 = __FADD( t1, c1 ); // Tan(x)+Cot(x) needed
1214 p2 = __FADD( t2, c2 ); // for all quadrants
1215 p3 = c1 - p1 + t1;
1216 p4 = c2 - p2 + t2;
1217 p5 = __FADD( p3, p2 );
1218 p6 = p3 - p5 + p2 + p4; // p1+p5+p6=Tan(x)+Cot(x) exactly
1219 // where x = table lookup argument
1220 // Next, this sum needs to be multiplied
1221 // by the tangent of the reduced argument
1222 p7 = p1*tan;
1223 p8 = __FMSUB( p1, tan, p7 ); // p1*tan - p7;
1224 p9 = p5*tan;
1225 p10 = __FMSUB( p5, tan, p9 ); // p5*tan - p9;
1226 p11 = p1*tanl;
1227 p12 = __FMSUB( p1, tanl, p11 ); // p1*tanl - p11;
1228 p13 = __FMADD( p5, tanl, p10 + p12 ); // p10 + p12 + p5*tanl;
1229 p14 = __FMADD( p6, tan, p13 ); // p13 + p6*tan;
1230 p15 = __FADD( p8, p9 );
1231 p17 = __FADD( p15, p11 );
1232
1233 if (fabs(p8) > fabs(p9))
1234 p16 = p8 - p15 + p9;
1235 else
1236 p16 = p9 - p15 + p8;
1237
1238 if (fabs(p11) > fabs(p15))
1239 p18 = p11 - p17 + p15 + (p14+p16);
1240 else
1241 p18 = p15 - p17 + p11 + (p14+p16);
1242
1243 p19 = __FADD( p7, p17 );
1244 p20 = p7 - p19 + p17;
1245 p21 = __FADD( p20, p18 );
1246 p22 = p20 - p21 + p18;
1247
1248 //***********************************************************************
1249 // At this point, different paths are taken for arguments *
1250 // "neighborhood" of even and odd multiples of pi/2. *
1251 //***********************************************************************
1252
1253 if ((M & 1) == 0) {
1254
1255 //***********************************************************************
1256 // Let the argument A be written as A=x+y, where *
1257 // x=TanTablePtr[index].arg. Then if T=tan(x) and C=cot(x) *
1258 // as read from tantbl, and t=tan(y), (computed above as *
1259 // (tan+tanl), then compute tan(A) as: *
1260 // *
1261 // t ( C+T) *
1262 // tan(A)= T + ----------- *
1263 // C-t *
1264 // *
1265 // numerator represented as (p19,p21,p22) in the code below. *
1266 // *
1267 //***********************************************************************
1268 p23 = __FSUB( c1, tan ); // High order part of C-t
1269 p24 = __FSUB( c2, tanl );
1270 rec = 1.0/p23;
1271 p25 = c1 - p23 - tan;
1272 frac = p19*rec; // start computation of fraction
1273
1274 if (fabs(c2) > fabs(tanl))
1275 p26 = c2 - p24 - tanl;
1276 else
1277 p26 = c2 - (p24+tanl);
1278
1279 p27 = __FADD( p25, p24 );
1280 p28 = p25 - p27 + p24 + p26; // C-t=(p23,p27,p28) done
1281 frac = frac + __FNMSUB( p23, frac, p19 )*rec; // frac + (p19-p23*frac)*rec; // refine lead fraction wd.
1282 p29 = __FNMSUB( p23, frac, p19 ); //p19 - p23*frac;
1283 p30 = frac*p27; // p30 a negative term
1284 p31 = __FNMSUB( frac, p27, p30 ) + p22; // p30 - frac*p27 + p22;
1285 p32 = __FNMSUB( frac, p28, p31 ); // p31 - frac*p28;
1286 p33 = __FSUB( p29, p30 );
1287 p35 = p33 + p21;
1288 frac2 = p35*rec;
1289 frac2 = frac2 + __FNMSUB( frac2, p23, p35)*rec; // frac2 + (p35-frac2*p23)*rec; // second fraction wd.
1290
1291 if (fabs(p29) > fabs(p30))
1292 p34 = p29 - p33 - p30 + p32;
1293 else
1294 p34 = p29 - (p33+p30) + p32;
1295
1296 p36 = __FNMSUB( frac2, p23, p35 ) + __FNMSUB( frac2, p27, p34 ); // p35 - frac2*p23 + p34 - frac2*p27;
1297
1298 if (fabs(p33) > fabs(p21))
1299 p37 = p33 - p35 + p21;
1300 else
1301 p37 = p21 - p35 + p33;
1302
1303 p38 = p37 + p36;
1304 frac3 = p38*rec; // Fraction of above formula done
1305
1306 //***********************************************************************
1307 // Finally, add (tan,tanl) and (frac,frac2,frac3) for result *
1308 //***********************************************************************
1309
1310 p39 = __FADD( t1, frac );
1311 p40 = t1 - p39 + frac;
1312 p41 = __FADD( t2, frac2 );
1313
1314 if (fabs(t2) > fabs(frac2))
1315 p42 = t2 - p41 + frac2 + frac3;
1316 else
1317 p42 = frac2 - p41 + t2 + frac3;
1318
1319 p43 = __FADD( p40, p41 );
1320 p44 = p40 - p43 + p41 + p42;
1321 p45 = __FADD( p39, p43 );
1322 p46 = p39 - p45 + p43 + p44;
1323
1324 // getting a cannonical result
1325
1326 ddx.d[0] = __FADD( p45, p46 );
1327 ddx.d[1] = (p45 - ddx.d[0]) + p46;
1328
1329 FESETENVD(fenv);
1330 return (ddx.f);
1331 }
1332 else {
1333
1334 //***********************************************************************
1335 // The remaining case is for argments in the neighborhood of *
1336 // even multiples of pi/2. *
1337 // Use the identity *
1338 // tan(A+(2k+1)pi/2)=-cot(A) *
1339 // *
1340 // Let the reduced argument A be written as *
1341 // A=x+y, where x= TanTablePtr[index].arg *
1342 // Then if T=tan(x) and C=cot(x) as read from tantbl, and t=tan(y),*
1343 // (computed above as (tan+tanl), compute cot(A) as: *
1344 // *
1345 // t(C+T) *
1346 // cot(A)=C - ---------- *
1347 // T+t *
1348 //***********************************************************************
1349 p23 = __FADD( t1, tan ); // High order part of T+t
1350 p24 = __FADD( t2, tanl );
1351 rec = 1.0/p23;
1352 p25 = t1 - p23 + tan;
1353 frac = p19*rec; // start computation of fraction
1354
1355 if (fabs(t2) > fabs(tanl))
1356 p26 = t2 - p24 + tanl;
1357 else
1358 p26 = tanl - p24 + t2;
1359
1360 p27 = __FADD( p25, p24 );
1361 p28 = p25 - p27 + p24 + p26; // T+t= (p23,p27,p28) done
1362 frac = frac + __FNMSUB( p23, frac, p19 )*rec; // frac + (p19-p23*frac)*rec; // refine lead fraction wd.
1363 p29 = __FNMSUB( p23, frac, p19 ); // p19 - p23*frac;
1364 p30 = frac*p27; // p30 a negative term
1365 p31 = __FNMSUB( frac, p27, p30 ) + p22; // p30 - frac*p27 + p22;
1366 p32 = __FNMSUB( frac, p28, p31 ); // p31 - frac*p28;
1367 p33 = __FSUB( p29, p30 );
1368 p35 = __FADD( p33, p21 );
1369 frac2 = p35*rec;
1370 frac2 = frac2, __FNMSUB( frac2, p23, p35 )*rec; // frac2 + (p35-frac2*p23)*rec; // second fraction wd.
1371
1372 if (fabs(p29) > fabs(p30))
1373 p34 = p29 - p33 - p30 + p32;
1374 else
1375 p34 = p29 - (p33+p30) + p32;
1376
1377 p36 = __FNMSUB( frac2, p23, p35 ) + __FNMSUB( frac2, p27, p34 ); // p35 - frac2*p23 + p34 - frac2*p27;
1378
1379 if (fabs(p33) > fabs(p21))
1380 p37 = p33 - p35 + p21;
1381 else
1382 p37 = p21 - p35 + p33;
1383
1384 p38 = p37 + p36;
1385 frac3 = p38*rec; // Fraction done
1386
1387 //*************************************************************************
1388 // At last, subtract (cot,cotl) from (frac,frac2,frac3) for the result *
1389 //*************************************************************************
1390 p39 = __FSUB( frac, c1 );
1391 p40 = frac - (c1+p39);
1392 p41 = __FSUB( frac2, c2 );
1393 if (fabs(c2) > fabs(frac2))
1394 p42 = frac2 - (c2+p41) + frac3;
1395 else
1396 p42 = frac2 - p41 - c2 + frac3;
1397
1398 p43 = __FADD( p40, p41 );
1399 p44 = p40 - p43 + p41 + p42;
1400 p45 = __FADD( p39, p43 );
1401 p46 = p39 - p45 + p43 + p44;
1402
1403 //getting a cannonical result
1404
1405 ddx.d[0] = __FADD( p45, p46 );
1406 ddx.d[1] = (p45 - ddx.d[0]) + p46;
1407
1408 FESETENVD(fenv);
1409 return (ddx.f);
1410 }
1411 }
1412}
1413#endif