VirtualBox

source: vbox/trunk/src/libs/softfloat-3e/source/extF80_sincos.c

Last change on this file was 98103, checked in by vboxsync, 17 months ago

Copyright year updates by scm.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.2 KB
Line 
1/* $Id: extF80_sincos.c 98103 2023-01-17 14:15:46Z vboxsync $ */
2/** @file
3 * SoftFloat - VBox Extension - extF80_sin, extF80_cos, extF80_sincos, extF80_atan2.
4 */
5
6/*
7 * Copyright (C) 2022-2023 Oracle and/or its affiliates.
8 *
9 * This file is part of VirtualBox base platform packages, as
10 * available from https://www.virtualbox.org.
11 *
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation, in version 3 of the
15 * License.
16 *
17 * This program is distributed in the hope that it will be useful, but
18 * WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with this program; if not, see <https://www.gnu.org/licenses>.
24 *
25 * SPDX-License-Identifier: GPL-3.0-only
26 */
27
28
29/*********************************************************************************************************************************
30* Header Files *
31*********************************************************************************************************************************/
32#include <stdbool.h>
33#include <stdint.h>
34#include "platform.h"
35#include "internals.h"
36#include "specialize.h"
37#include "softfloat.h"
38#include <iprt/types.h>
39#include <iprt/x86.h>
40
41#include "extF80_sincos.h"
42
43
44static void cordic_sincos( float128_t z, float128_t *pv1, float128_t *pv2 SOFTFLOAT_STATE_DECL_COMMA )
45{
46 float128_t v1 = { { 0, 0 } }; /* MSC thinks it can be used uninitialized */
47 float128_t v2 = { { 0, 0 } }; /* MSC thinks it can be used uninitialized */
48 /** @todo TBD: CORDIC kernel should be easily implemented in assembly * */
49
50 float128_t x1 = ui32_to_f128(1, pState);
51 float128_t x2 = ui32_to_f128(0, pState);
52 float128_t zz = ui32_to_f128(0, pState);
53
54 float128_t p2m = ui32_to_f128(1, pState);
55 float128_t two = ui32_to_f128(2, pState);
56
57 for (unsigned k = 0; k < RT_ELEMENTS(g_ar128FsincosCORDICConsts); k++)
58 {
59 float128_t atg = *(float128_t *)&g_ar128FsincosCORDICConsts[k];
60 float128_t scale = *(float128_t *)&g_ar128FsincosCORDICConsts2[k];
61
62 float128_t px1 = f128_mul(x1, p2m, pState);
63 float128_t px2 = f128_mul(x2, p2m, pState);
64
65 if (f128_le(zz, z, pState))
66 {
67 x1 = f128_sub(x1, px2, pState);
68 x2 = f128_add(x2, px1, pState);
69 zz = f128_add(zz, atg, pState);
70 }
71 else
72 {
73 x1 = f128_add(x1, px2, pState);
74 x2 = f128_sub(x2, px1, pState);
75 zz = f128_sub(zz, atg, pState);
76 }
77
78 p2m = f128_div(p2m, two, pState);
79
80 v1 = f128_mul(x1, scale, pState);
81 v2 = f128_mul(x2, scale, pState);
82 }
83
84 *pv1 = v1;
85 *pv2 = v2;
86}
87
88static float128_t cordic_atan2( float128_t y, float128_t x SOFTFLOAT_STATE_DECL_COMMA )
89{
90 float128_t v1 = { { 0, 0 } }; /* MSC thinks it can be used uninitialized */
91 float128_t v2 = { { 0, 0 } }; /* MSC thinks it can be used uninitialized */
92 /** @todo TBD: CORDIC kernel should be easily implemented in assembly * */
93
94 float128_t x1 = x, x2 = y;
95 float128_t z = ui32_to_f128(0, pState);
96 float128_t zero = ui32_to_f128(0, pState);
97 float128_t p2m = ui32_to_f128(1, pState);
98 float128_t two = ui32_to_f128(2, pState);
99
100 for (unsigned k = 0; k < RT_ELEMENTS(g_ar128FsincosCORDICConsts); k++)
101 {
102 float128_t atg = *(float128_t *)&g_ar128FsincosCORDICConsts[k];
103 float128_t scale = *(float128_t *)&g_ar128FsincosCORDICConsts2[k];
104
105 float128_t px1 = f128_mul(x1, p2m, pState);
106 float128_t px2 = f128_mul(x2, p2m, pState);
107
108 if (f128_le(x2, zero, pState))
109 {
110 x1 = f128_sub(x1, px2, pState);
111 x2 = f128_add(x2, px1, pState);
112 z = f128_sub(z, atg, pState);
113 }
114 else
115 {
116 x1 = f128_add(x1, px2, pState);
117 x2 = f128_sub(x2, px1, pState);
118 z = f128_add(z, atg, pState);
119 }
120
121 p2m = f128_div(p2m, two, pState);
122
123 v1 = f128_mul(x1, scale, pState);
124 v2 = f128_mul(x2, scale, pState);
125 }
126
127 return z;
128}
129
130extFloat80_t extF80_sin( extFloat80_t x SOFTFLOAT_STATE_DECL_COMMA )
131{
132 int32_t fSign = 0;
133 extFloat80_t f80zero = ui32_to_extF80(0, pState);
134 if (extF80_le(x, f80zero, pState))
135 {
136 x = extF80_sub(f80zero, x, pState);
137 fSign = 1;
138 }
139
140 extFloat80_t f80pi2 = f128_to_extF80(*(float128_t *)&g_r128pi2, pState);
141
142 /** @todo TBD: Partial remainder should be calculated using float128 value of pi2 to increase precision **/
143 uint16_t fCxFlags = 0;
144 extFloat80_t rem = extF80_partialRem(x, f80pi2, pState->roundingMode, &fCxFlags, pState);
145 int32_t const quo = X86_FSW_CX_TO_QUOTIENT(fCxFlags);
146
147 float128_t z = extF80_to_f128(rem, pState);
148 float128_t f128zero = ui32_to_f128(0, pState);
149
150 float128_t v1, v2;
151 cordic_sincos(z, &v1, &v2, pState);
152
153 float128_t v;
154 switch(quo % 4)
155 {
156#ifdef _MSC_VER /* stupid MSC thinks v might be used uninitialized otherwise: */
157 default:
158#endif
159 case 0:
160 v = v2;
161 break;
162
163 case 1:
164 v = v1;
165 break;
166
167 case 2:
168 v = f128_sub(f128zero, v2, pState);
169 break;
170
171 case 3:
172 v = f128_sub(f128zero, v1, pState);
173 break;
174 }
175
176 if (fSign)
177 v = f128_sub(f128zero, v, pState);
178
179 return f128_to_extF80(v, pState);
180}
181
182extFloat80_t extF80_cos( extFloat80_t x SOFTFLOAT_STATE_DECL_COMMA )
183{
184 extFloat80_t f80zero = ui32_to_extF80(0, pState);
185 if (extF80_le(x, f80zero, pState))
186 x = extF80_sub(f80zero, x, pState);
187
188 extFloat80_t f80pi2 = f128_to_extF80(*(float128_t *)&g_r128pi2, pState);
189
190 /** TBD: Partial remainder should be calculated using float128 value of pi2 to increase precision **/
191 uint16_t fCxFlags = 0;
192 extFloat80_t rem = extF80_partialRem(x, f80pi2, pState->roundingMode, &fCxFlags, pState);
193 int32_t const quo = X86_FSW_CX_TO_QUOTIENT(fCxFlags);
194
195 float128_t z = extF80_to_f128(rem, pState);
196 float128_t f128zero = ui32_to_f128(0, pState);
197
198 float128_t v1, v2;
199 cordic_sincos(z, &v1, &v2, pState);
200
201 float128_t v;
202 switch(quo % 4)
203 {
204#ifdef _MSC_VER /* stupid MSC thinks v might be used uninitialized otherwise: */
205 default:
206#endif
207 case 0:
208 v = v1;
209 break;
210
211 case 1:
212 v = f128_sub(f128zero, v2, pState);;
213 break;
214
215 case 2:
216 v = f128_sub(f128zero, v1, pState);
217 break;
218
219 case 3:
220 v = v2;
221 break;
222 }
223
224 return f128_to_extF80(v, pState);
225}
226
227void extF80_sincos( extFloat80_t x, extFloat80_t* pSin, extFloat80_t* pCos SOFTFLOAT_STATE_DECL_COMMA )
228{
229 uint16_t fCxFlags = 0;
230 int32_t quo;
231 extFloat80_t rem, f80pi2, f80zero;
232 int32_t fSign = 0;
233
234 f80zero = ui32_to_extF80(0, pState);
235 if (extF80_le(x, f80zero, pState))
236 {
237 x = extF80_sub(f80zero, x, pState);
238 fSign = 1;
239 }
240
241 f80pi2 = f128_to_extF80(*(float128_t const *)&g_r128pi2, pState);
242
243 /** @todo TBD: Partial remainder should be calculated using float128 value of pi2 to increase precision **/
244 rem = extF80_partialRem(x, f80pi2, pState->roundingMode, &fCxFlags, pState);
245 quo = X86_FSW_CX_TO_QUOTIENT(fCxFlags);
246
247 float128_t z = extF80_to_f128(rem, pState);
248 float128_t f128zero = ui32_to_f128(0, pState);
249
250 float128_t v1, v2;
251 cordic_sincos(z, &v1, &v2, pState);
252
253 float128_t vCos, vSin;
254 switch(quo % 4)
255 {
256#ifdef _MSC_VER /* stupid MSC thinks vCos & cSin might be used uninitialized otherwise: */
257 default:
258#endif
259 case 0:
260 vCos = v1;
261 vSin = v2;
262 break;
263
264 case 1:
265 vCos = f128_sub(f128zero, v2, pState);
266 vSin = v1;
267 break;
268
269 case 2:
270 vCos = f128_sub(f128zero, v1, pState);
271 vSin = f128_sub(f128zero, v2, pState);
272 break;
273
274 case 3:
275 vCos = v2;
276 vSin = f128_sub(f128zero, v1, pState);
277 break;
278 }
279
280 if (fSign)
281 vSin = f128_sub(f128zero, vSin, pState);
282
283 *pCos = f128_to_extF80(vCos, pState);
284 *pSin = f128_to_extF80(vSin, pState);
285}
286
287extFloat80_t extF80_atan2( extFloat80_t f80y, extFloat80_t f80x SOFTFLOAT_STATE_DECL_COMMA )
288{
289 float128_t v;
290 int32_t fSignX = 0, fSignY = 0;
291 float128_t f128zero = ui32_to_f128(0, pState);
292 float128_t y = extF80_to_f128(f80y, pState);
293 float128_t x = extF80_to_f128(f80x, pState);
294
295 if (f128_le(x, f128zero, pState))
296 {
297 x = f128_sub(f128zero, x, pState);
298 fSignX = 1;
299 }
300
301 if (f128_le(y, f128zero, pState))
302 {
303 y = f128_sub(f128zero, y, pState);
304 fSignY = 1;
305 }
306
307 v = cordic_atan2(y, x, pState);
308
309 if (fSignX)
310 {
311 if (fSignY)
312 v = f128_sub(v, *(float128_t const *)&g_r128pi, pState);
313 else
314 v = f128_sub(*(float128_t const *)&g_r128pi, v, pState);
315 }
316 else
317 {
318 if (fSignY)
319 v = f128_sub(f128zero, v, pState);
320 }
321
322 return f128_to_extF80(v, pState);
323}
Note: See TracBrowser for help on using the repository browser.

© 2023 Oracle
Contact – Privacy policy – Terms of Use