source: mainline/uspace/lib/softfloat/div.c@ c0c38c7c

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since c0c38c7c was c0c38c7c, checked in by Martin Decky <martin@…>, 10 years ago

software floating point overhaul
use proper type mapping
fix cosine calculation
  • Property mode set to 100644
File size: 13.5 KB
Line 
1/*
2 * Copyright (c) 2005 Josef Cejka
3 * Copyright (c) 2011 Petr Koupy
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * - Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * - Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * - The name of the author may not be used to endorse or promote products
16 * derived from this software without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
19 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
20 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
22 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
23 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
25 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
26 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
27 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30/** @addtogroup softfloat
31 * @{
32 */
33/** @file Division functions.
34 */
35
36#include "add.h"
37#include "div.h"
38#include "comparison.h"
39#include "mul.h"
40#include "common.h"
41
42/** Divide two single-precision floats.
43 *
44 * @param a Nominator.
45 * @param b Denominator.
46 *
47 * @return Result of division.
48 *
49 */
50float32 div_float32(float32 a, float32 b)
51{
52 float32 result;
53 int32_t aexp, bexp, cexp;
54 uint64_t afrac, bfrac, cfrac;
55
56 result.parts.sign = a.parts.sign ^ b.parts.sign;
57
58 if (is_float32_nan(a)) {
59 if (is_float32_signan(a)) {
60 // FIXME: SigNaN
61 }
62 /* NaN */
63 return a;
64 }
65
66 if (is_float32_nan(b)) {
67 if (is_float32_signan(b)) {
68 // FIXME: SigNaN
69 }
70 /* NaN */
71 return b;
72 }
73
74 if (is_float32_infinity(a)) {
75 if (is_float32_infinity(b)) {
76 /*FIXME: inf / inf */
77 result.bin = FLOAT32_NAN;
78 return result;
79 }
80 /* inf / num */
81 result.parts.exp = a.parts.exp;
82 result.parts.fraction = a.parts.fraction;
83 return result;
84 }
85
86 if (is_float32_infinity(b)) {
87 if (is_float32_zero(a)) {
88 /* FIXME 0 / inf */
89 result.parts.exp = 0;
90 result.parts.fraction = 0;
91 return result;
92 }
93 /* FIXME: num / inf*/
94 result.parts.exp = 0;
95 result.parts.fraction = 0;
96 return result;
97 }
98
99 if (is_float32_zero(b)) {
100 if (is_float32_zero(a)) {
101 /*FIXME: 0 / 0*/
102 result.bin = FLOAT32_NAN;
103 return result;
104 }
105 /* FIXME: division by zero */
106 result.parts.exp = 0;
107 result.parts.fraction = 0;
108 return result;
109 }
110
111 afrac = a.parts.fraction;
112 aexp = a.parts.exp;
113 bfrac = b.parts.fraction;
114 bexp = b.parts.exp;
115
116 /* denormalized numbers */
117 if (aexp == 0) {
118 if (afrac == 0) {
119 result.parts.exp = 0;
120 result.parts.fraction = 0;
121 return result;
122 }
123
124 /* normalize it*/
125 afrac <<= 1;
126 /* afrac is nonzero => it must stop */
127 while (!(afrac & FLOAT32_HIDDEN_BIT_MASK)) {
128 afrac <<= 1;
129 aexp--;
130 }
131 }
132
133 if (bexp == 0) {
134 bfrac <<= 1;
135 /* bfrac is nonzero => it must stop */
136 while (!(bfrac & FLOAT32_HIDDEN_BIT_MASK)) {
137 bfrac <<= 1;
138 bexp--;
139 }
140 }
141
142 afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE - 1);
143 bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE);
144
145 if (bfrac <= (afrac << 1)) {
146 afrac >>= 1;
147 aexp++;
148 }
149
150 cexp = aexp - bexp + FLOAT32_BIAS - 2;
151
152 cfrac = (afrac << 32) / bfrac;
153 if ((cfrac & 0x3F) == 0) {
154 cfrac |= (bfrac * cfrac != afrac << 32);
155 }
156
157 /* pack and round */
158
159 /* find first nonzero digit and shift result and detect possibly underflow */
160 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)))) {
161 cexp--;
162 cfrac <<= 1;
163 /* TODO: fix underflow */
164 }
165
166 cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
167
168 if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
169 ++cexp;
170 cfrac >>= 1;
171 }
172
173 /* check overflow */
174 if (cexp >= FLOAT32_MAX_EXPONENT) {
175 /* FIXME: overflow, return infinity */
176 result.parts.exp = FLOAT32_MAX_EXPONENT;
177 result.parts.fraction = 0;
178 return result;
179 }
180
181 if (cexp < 0) {
182 /* FIXME: underflow */
183 result.parts.exp = 0;
184 if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
185 result.parts.fraction = 0;
186 return result;
187 }
188 cfrac >>= 1;
189 while (cexp < 0) {
190 cexp++;
191 cfrac >>= 1;
192 }
193 } else {
194 result.parts.exp = (uint32_t) cexp;
195 }
196
197 result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
198
199 return result;
200}
201
202/** Divide two double-precision floats.
203 *
204 * @param a Nominator.
205 * @param b Denominator.
206 *
207 * @return Result of division.
208 *
209 */
210float64 div_float64(float64 a, float64 b)
211{
212 float64 result;
213 int64_t aexp, bexp, cexp;
214 uint64_t afrac, bfrac, cfrac;
215 uint64_t remlo, remhi;
216 uint64_t tmplo, tmphi;
217
218 result.parts.sign = a.parts.sign ^ b.parts.sign;
219
220 if (is_float64_nan(a)) {
221 if (is_float64_signan(b)) {
222 // FIXME: SigNaN
223 return b;
224 }
225
226 if (is_float64_signan(a)) {
227 // FIXME: SigNaN
228 }
229 /* NaN */
230 return a;
231 }
232
233 if (is_float64_nan(b)) {
234 if (is_float64_signan(b)) {
235 // FIXME: SigNaN
236 }
237 /* NaN */
238 return b;
239 }
240
241 if (is_float64_infinity(a)) {
242 if (is_float64_infinity(b) || is_float64_zero(b)) {
243 // FIXME: inf / inf
244 result.bin = FLOAT64_NAN;
245 return result;
246 }
247 /* inf / num */
248 result.parts.exp = a.parts.exp;
249 result.parts.fraction = a.parts.fraction;
250 return result;
251 }
252
253 if (is_float64_infinity(b)) {
254 if (is_float64_zero(a)) {
255 /* FIXME 0 / inf */
256 result.parts.exp = 0;
257 result.parts.fraction = 0;
258 return result;
259 }
260 /* FIXME: num / inf*/
261 result.parts.exp = 0;
262 result.parts.fraction = 0;
263 return result;
264 }
265
266 if (is_float64_zero(b)) {
267 if (is_float64_zero(a)) {
268 /*FIXME: 0 / 0*/
269 result.bin = FLOAT64_NAN;
270 return result;
271 }
272 /* FIXME: division by zero */
273 result.parts.exp = 0;
274 result.parts.fraction = 0;
275 return result;
276 }
277
278 afrac = a.parts.fraction;
279 aexp = a.parts.exp;
280 bfrac = b.parts.fraction;
281 bexp = b.parts.exp;
282
283 /* denormalized numbers */
284 if (aexp == 0) {
285 if (afrac == 0) {
286 result.parts.exp = 0;
287 result.parts.fraction = 0;
288 return result;
289 }
290
291 /* normalize it*/
292 aexp++;
293 /* afrac is nonzero => it must stop */
294 while (!(afrac & FLOAT64_HIDDEN_BIT_MASK)) {
295 afrac <<= 1;
296 aexp--;
297 }
298 }
299
300 if (bexp == 0) {
301 bexp++;
302 /* bfrac is nonzero => it must stop */
303 while (!(bfrac & FLOAT64_HIDDEN_BIT_MASK)) {
304 bfrac <<= 1;
305 bexp--;
306 }
307 }
308
309 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 2);
310 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 1);
311
312 if (bfrac <= (afrac << 1)) {
313 afrac >>= 1;
314 aexp++;
315 }
316
317 cexp = aexp - bexp + FLOAT64_BIAS - 2;
318
319 cfrac = div128est(afrac, 0x0ll, bfrac);
320
321 if ((cfrac & 0x1FF) <= 2) {
322 mul64(bfrac, cfrac, &tmphi, &tmplo);
323 sub128(afrac, 0x0ll, tmphi, tmplo, &remhi, &remlo);
324
325 while ((int64_t) remhi < 0) {
326 cfrac--;
327 add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo);
328 }
329 cfrac |= (remlo != 0);
330 }
331
332 /* round and shift */
333 result = finish_float64(cexp, cfrac, result.parts.sign);
334 return result;
335}
336
337/** Divide two quadruple-precision floats.
338 *
339 * @param a Nominator.
340 * @param b Denominator.
341 *
342 * @return Result of division.
343 *
344 */
345float128 div_float128(float128 a, float128 b)
346{
347 float128 result;
348 int64_t aexp, bexp, cexp;
349 uint64_t afrac_hi, afrac_lo, bfrac_hi, bfrac_lo, cfrac_hi, cfrac_lo;
350 uint64_t shift_out;
351 uint64_t rem_hihi, rem_hilo, rem_lohi, rem_lolo;
352 uint64_t tmp_hihi, tmp_hilo, tmp_lohi, tmp_lolo;
353
354 result.parts.sign = a.parts.sign ^ b.parts.sign;
355
356 if (is_float128_nan(a)) {
357 if (is_float128_signan(b)) {
358 // FIXME: SigNaN
359 return b;
360 }
361
362 if (is_float128_signan(a)) {
363 // FIXME: SigNaN
364 }
365 /* NaN */
366 return a;
367 }
368
369 if (is_float128_nan(b)) {
370 if (is_float128_signan(b)) {
371 // FIXME: SigNaN
372 }
373 /* NaN */
374 return b;
375 }
376
377 if (is_float128_infinity(a)) {
378 if (is_float128_infinity(b) || is_float128_zero(b)) {
379 // FIXME: inf / inf
380 result.bin.hi = FLOAT128_NAN_HI;
381 result.bin.lo = FLOAT128_NAN_LO;
382 return result;
383 }
384 /* inf / num */
385 result.parts.exp = a.parts.exp;
386 result.parts.frac_hi = a.parts.frac_hi;
387 result.parts.frac_lo = a.parts.frac_lo;
388 return result;
389 }
390
391 if (is_float128_infinity(b)) {
392 if (is_float128_zero(a)) {
393 // FIXME 0 / inf
394 result.parts.exp = 0;
395 result.parts.frac_hi = 0;
396 result.parts.frac_lo = 0;
397 return result;
398 }
399 // FIXME: num / inf
400 result.parts.exp = 0;
401 result.parts.frac_hi = 0;
402 result.parts.frac_lo = 0;
403 return result;
404 }
405
406 if (is_float128_zero(b)) {
407 if (is_float128_zero(a)) {
408 // FIXME: 0 / 0
409 result.bin.hi = FLOAT128_NAN_HI;
410 result.bin.lo = FLOAT128_NAN_LO;
411 return result;
412 }
413 // FIXME: division by zero
414 result.parts.exp = 0;
415 result.parts.frac_hi = 0;
416 result.parts.frac_lo = 0;
417 return result;
418 }
419
420 afrac_hi = a.parts.frac_hi;
421 afrac_lo = a.parts.frac_lo;
422 aexp = a.parts.exp;
423 bfrac_hi = b.parts.frac_hi;
424 bfrac_lo = b.parts.frac_lo;
425 bexp = b.parts.exp;
426
427 /* denormalized numbers */
428 if (aexp == 0) {
429 if (eq128(afrac_hi, afrac_lo, 0x0ll, 0x0ll)) {
430 result.parts.exp = 0;
431 result.parts.frac_hi = 0;
432 result.parts.frac_lo = 0;
433 return result;
434 }
435
436 /* normalize it*/
437 aexp++;
438 /* afrac is nonzero => it must stop */
439 and128(afrac_hi, afrac_lo,
440 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
441 &tmp_hihi, &tmp_lolo);
442 while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
443 lshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
444 aexp--;
445 }
446 }
447
448 if (bexp == 0) {
449 bexp++;
450 /* bfrac is nonzero => it must stop */
451 and128(bfrac_hi, bfrac_lo,
452 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
453 &tmp_hihi, &tmp_lolo);
454 while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
455 lshift128(bfrac_hi, bfrac_lo, 1, &bfrac_hi, &bfrac_lo);
456 bexp--;
457 }
458 }
459
460 or128(afrac_hi, afrac_lo,
461 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
462 &afrac_hi, &afrac_lo);
463 lshift128(afrac_hi, afrac_lo,
464 (128 - FLOAT128_FRACTION_SIZE - 1), &afrac_hi, &afrac_lo);
465 or128(bfrac_hi, bfrac_lo,
466 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
467 &bfrac_hi, &bfrac_lo);
468 lshift128(bfrac_hi, bfrac_lo,
469 (128 - FLOAT128_FRACTION_SIZE - 1), &bfrac_hi, &bfrac_lo);
470
471 if (le128(bfrac_hi, bfrac_lo, afrac_hi, afrac_lo)) {
472 rshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
473 aexp++;
474 }
475
476 cexp = aexp - bexp + FLOAT128_BIAS - 2;
477
478 cfrac_hi = div128est(afrac_hi, afrac_lo, bfrac_hi);
479
480 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_hi,
481 &tmp_lolo /* dummy */, &tmp_hihi, &tmp_hilo, &tmp_lohi);
482
483 /* sub192(afrac_hi, afrac_lo, 0,
484 * tmp_hihi, tmp_hilo, tmp_lohi
485 * &rem_hihi, &rem_hilo, &rem_lohi); */
486 sub128(afrac_hi, afrac_lo, tmp_hihi, tmp_hilo, &rem_hihi, &rem_hilo);
487 if (tmp_lohi > 0) {
488 sub128(rem_hihi, rem_hilo, 0x0ll, 0x1ll, &rem_hihi, &rem_hilo);
489 }
490 rem_lohi = -tmp_lohi;
491
492 while ((int64_t) rem_hihi < 0) {
493 --cfrac_hi;
494 /* add192(rem_hihi, rem_hilo, rem_lohi,
495 * 0, bfrac_hi, bfrac_lo,
496 * &rem_hihi, &rem_hilo, &rem_lohi); */
497 add128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo, &rem_hilo, &rem_lohi);
498 if (lt128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo)) {
499 ++rem_hihi;
500 }
501 }
502
503 cfrac_lo = div128est(rem_hilo, rem_lohi, bfrac_lo);
504
505 if ((cfrac_lo & 0x3FFF) <= 4) {
506 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_lo,
507 &tmp_hihi /* dummy */, &tmp_hilo, &tmp_lohi, &tmp_lolo);
508
509 /* sub192(rem_hilo, rem_lohi, 0,
510 * tmp_hilo, tmp_lohi, tmp_lolo,
511 * &rem_hilo, &rem_lohi, &rem_lolo); */
512 sub128(rem_hilo, rem_lohi, tmp_hilo, tmp_lohi, &rem_hilo, &rem_lohi);
513 if (tmp_lolo > 0) {
514 sub128(rem_hilo, rem_lohi, 0x0ll, 0x1ll, &rem_hilo, &rem_lohi);
515 }
516 rem_lolo = -tmp_lolo;
517
518 while ((int64_t) rem_hilo < 0) {
519 --cfrac_lo;
520 /* add192(rem_hilo, rem_lohi, rem_lolo,
521 * 0, bfrac_hi, bfrac_lo,
522 * &rem_hilo, &rem_lohi, &rem_lolo); */
523 add128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo, &rem_lohi, &rem_lolo);
524 if (lt128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo)) {
525 ++rem_hilo;
526 }
527 }
528
529 cfrac_lo |= ((rem_hilo | rem_lohi | rem_lolo) != 0 );
530 }
531
532 shift_out = cfrac_lo << (64 - (128 - FLOAT128_FRACTION_SIZE - 1));
533 rshift128(cfrac_hi, cfrac_lo, (128 - FLOAT128_FRACTION_SIZE - 1),
534 &cfrac_hi, &cfrac_lo);
535
536 result = finish_float128(cexp, cfrac_hi, cfrac_lo, result.parts.sign, shift_out);
537 return result;
538}
539
540#ifdef float32_t
541
542float32_t __divsf3(float32_t a, float32_t b)
543{
544 float32_u ua;
545 ua.val = a;
546
547 float32_u ub;
548 ub.val = b;
549
550 float32_u res;
551 res.data = div_float32(ua.data, ub.data);
552
553 return res.val;
554}
555
556float32_t __aeabi_fdiv(float32_t a, float32_t b)
557{
558 float32_u ua;
559 ua.val = a;
560
561 float32_u ub;
562 ub.val = b;
563
564 float32_u res;
565 res.data = div_float32(ua.data, ub.data);
566
567 return res.val;
568}
569
570#endif
571
572#ifdef float64_t
573
574float64_t __divdf3(float64_t a, float64_t b)
575{
576 float64_u ua;
577 ua.val = a;
578
579 float64_u ub;
580 ub.val = b;
581
582 float64_u res;
583 res.data = div_float64(ua.data, ub.data);
584
585 return res.val;
586}
587
588float64_t __aeabi_ddiv(float64_t a, float64_t b)
589{
590 float64_u ua;
591 ua.val = a;
592
593 float64_u ub;
594 ub.val = b;
595
596 float64_u res;
597 res.data = div_float64(ua.data, ub.data);
598
599 return res.val;
600}
601
602#endif
603
604#ifdef float128_t
605
606float128_t __divtf3(float128_t a, float128_t b)
607{
608 float128_u ua;
609 ua.val = a;
610
611 float128_u ub;
612 ub.val = b;
613
614 float128_u res;
615 res.data = div_float128(ua.data, ub.data);
616
617 return res.val;
618}
619
620void _Qp_div(float128_t *c, float128_t *a, float128_t *b)
621{
622 *c = __divtf3(*a, *b);
623}
624
625#endif
626
627/** @}
628 */
Note: See TracBrowser for help on using the repository browser.