source: mainline/uspace/lib/softfloat/div.c@ 5ef16903

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since 5ef16903 was a35b458, checked in by Jiří Zárevúcky <zarevucky.jiri@…>, 8 years ago

style: Remove trailing whitespace on _all_ lines, including empty ones, for particular file types.

Command used: tools/srepl '\s\+$' '' -- *.c *.h *.py *.sh *.s *.S *.ag

Currently, whitespace on empty lines is very inconsistent.
There are two basic choices: Either remove the whitespace, or keep empty lines
indented to the level of surrounding code. The former is AFAICT more common,
and also much easier to do automatically.

Alternatively, we could write script for automatic indentation, and use that
instead. However, if such a script exists, it's possible to use the indented
style locally, by having the editor apply relevant conversions on load/save,
without affecting remote repository. IMO, it makes more sense to adopt
the simpler rule.
  • Property mode set to 100644
File size: 13.4 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.