| blob | 9a00dcea25ed32ca91f2253f76d6f878aa711d34 |
1 /* ieee754-df.S double-precision floating point support for ARM
3 Copyright (C) 2003 Free Software Foundation, Inc.
4 Contributed by Nicolas Pitre (nico@cam.org)
6 This file is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version.
11 In addition to the permissions in the GNU General Public License, the
12 Free Software Foundation gives you unlimited permission to link the
13 compiled version of this file into combinations with other programs,
14 and to distribute those combinations without any restriction coming
15 from the use of this file. (The General Public License restrictions
16 do apply in other respects; for example, they cover modification of
17 the file, and distribution when not linked into a combine
18 executable.)
20 This file is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of
22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23 General Public License for more details.
25 You should have received a copy of the GNU General Public License
26 along with this program; see the file COPYING. If not, write to
27 the Free Software Foundation, 59 Temple Place - Suite 330,
28 Boston, MA 02111-1307, USA. */
30 /*
31 * Notes:
32 *
33 * The goal of this code is to be as fast as possible. This is
34 * not meant to be easy to understand for the casual reader.
35 * For slightly simpler code please see the single precision version
36 * of this file.
37 *
38 * Only the default rounding mode is intended for best performances.
39 * Exceptions aren't supported yet, but that can be added quite easily
40 * if necessary without impacting performances.
41 */
43 @ This selects the minimum architecture level required.
44 #undef __ARM_ARCH__
45 #define __ARM_ARCH__ 3
47 #if defined(__ARM_ARCH_3M__) || defined(__ARM_ARCH_4__) \
48 || defined(__ARM_ARCH_4T__)
49 #undef __ARM_ARCH__
50 /* We use __ARM_ARCH__ set to 4 here, but in reality it's any processor with
51 long multiply instructions. That includes v3M. */
52 #define __ARM_ARCH__ 4
53 #endif
55 #if defined(__ARM_ARCH_5__) || defined(__ARM_ARCH_5T__) \
56 || defined(__ARM_ARCH_5TE__)
57 #undef __ARM_ARCH__
58 #define __ARM_ARCH__ 5
59 #endif
61 #if (__ARM_ARCH__ > 4) || defined(__ARM_ARCH_4T__)
62 #undef RET
63 #undef RETc
64 #define RET bx lr
65 #define RETc(x) bx##x lr
66 #if (__ARM_ARCH__ == 4) && (defined(__thumb__) || defined(__THUMB_INTERWORK__))
67 #define __FP_INTERWORKING__
68 #endif
69 #endif
71 @ For FPA, float words are always big-endian.
72 @ For VFP, floats words follow the memory system mode.
73 #if defined(__VFP_FP__) && !defined(__ARMEB__)
74 #define xl r0
75 #define xh r1
76 #define yl r2
77 #define yh r3
78 #else
79 #define xh r0
80 #define xl r1
81 #define yh r2
82 #define yl r3
83 #endif
86 #if defined(__thumb__) && !defined(__THUMB_INTERWORK__)
87 .macro ARM_FUNC_START name
88 FUNC_START \name
89 bx pc
90 nop
91 .arm
92 .endm
93 #else
94 .macro ARM_FUNC_START name
95 FUNC_START \name
96 .endm
97 #endif
99 ARM_FUNC_START negdf2
100 @ flip sign bit
101 eor xh, xh, #0x80000000
102 RET
104 ARM_FUNC_START subdf3
105 @ flip sign bit of second arg
106 eor yh, yh, #0x80000000
107 #if defined(__thumb__) && !defined(__THUMB_INTERWORK__)
108 b 1f @ Skip Thumb-code prologue
109 #endif
111 ARM_FUNC_START adddf3
113 1: @ Compare both args, return zero if equal but the sign.
114 teq xl, yl
115 eoreq ip, xh, yh
116 teqeq ip, #0x80000000
117 beq LSYM(Lad_z)
119 @ If first arg is 0 or -0, return second arg.
120 @ If second arg is 0 or -0, return first arg.
121 orrs ip, xl, xh, lsl #1
122 moveq xl, yl
123 moveq xh, yh
124 orrnes ip, yl, yh, lsl #1
125 RETc(eq)
127 stmfd sp!, {r4, r5, lr}
129 @ Mask out exponents.
130 mov ip, #0x7f000000
131 orr ip, ip, #0x00f00000
132 and r4, xh, ip
133 and r5, yh, ip
135 @ If either of them is 0x7ff, result will be INF or NAN
136 teq r4, ip
137 teqne r5, ip
138 beq LSYM(Lad_i)
140 @ Compute exponent difference. Make largest exponent in r4,
141 @ corresponding arg in xh-xl, and positive exponent difference in r5.
142 subs r5, r5, r4
143 rsblt r5, r5, #0
144 ble 1f
145 add r4, r4, r5
146 eor yl, xl, yl
147 eor yh, xh, yh
148 eor xl, yl, xl
149 eor xh, yh, xh
150 eor yl, xl, yl
151 eor yh, xh, yh
152 1:
154 @ If exponent difference is too large, return largest argument
155 @ already in xh-xl. We need up to 54 bit to handle proper rounding
156 @ of 0x1p54 - 1.1.
157 cmp r5, #(54 << 20)
158 #ifdef __FP_INTERWORKING__
159 ldmhifd sp!, {r4, r5, lr}
160 bxhi lr
161 #else
162 ldmhifd sp!, {r4, r5, pc}RETCOND
163 #endif
165 @ Convert mantissa to signed integer.
166 tst xh, #0x80000000
167 bic xh, xh, ip, lsl #1
168 orr xh, xh, #0x00100000
169 beq 1f
170 rsbs xl, xl, #0
171 rsc xh, xh, #0
172 1:
173 tst yh, #0x80000000
174 bic yh, yh, ip, lsl #1
175 orr yh, yh, #0x00100000
176 beq 1f
177 rsbs yl, yl, #0
178 rsc yh, yh, #0
179 1:
180 @ If exponent == difference, one or both args were denormalized.
181 @ Since this is not common case, rescale them off line.
182 teq r4, r5
183 beq LSYM(Lad_d)
184 LSYM(Lad_x):
185 @ Scale down second arg with exponent difference.
186 @ Apply shift one bit left to first arg and the rest to second arg
187 @ to simplify things later, but only if exponent does not become 0.
188 mov ip, #0
189 movs r5, r5, lsr #20
190 beq 3f
191 teq r4, #(1 << 20)
192 beq 1f
193 movs xl, xl, lsl #1
194 adc xh, ip, xh, lsl #1
195 sub r4, r4, #(1 << 20)
196 subs r5, r5, #1
197 beq 3f
199 @ Shift yh-yl right per r5, keep leftover bits into ip.
200 1: rsbs lr, r5, #32
201 blt 2f
202 mov ip, yl, lsl lr
203 mov yl, yl, lsr r5
204 orr yl, yl, yh, lsl lr
205 mov yh, yh, asr r5
206 b 3f
207 2: sub r5, r5, #32
208 add lr, lr, #32
209 cmp yl, #1
210 adc ip, ip, yh, lsl lr
211 mov yl, yh, asr r5
212 mov yh, yh, asr #32
213 3:
214 @ the actual addition
215 adds xl, xl, yl
216 adc xh, xh, yh
218 @ We now have a result in xh-xl-ip.
219 @ Keep absolute value in xh-xl-ip, sign in r5.
220 ands r5, xh, #0x80000000
221 bpl LSYM(Lad_p)
222 rsbs ip, ip, #0
223 rscs xl, xl, #0
224 rsc xh, xh, #0
226 @ Determine how to normalize the result.
227 LSYM(Lad_p):
228 cmp xh, #0x00100000
229 bcc LSYM(Lad_l)
230 cmp r0, #0x00200000
231 bcc LSYM(Lad_r0)
232 cmp r0, #0x00400000
233 bcc LSYM(Lad_r1)
235 @ Result needs to be shifted right.
236 movs xh, xh, lsr #1
237 movs xl, xl, rrx
238 movs ip, ip, rrx
239 orrcs ip, ip, #1
240 add r4, r4, #(1 << 20)
241 LSYM(Lad_r1):
242 movs xh, xh, lsr #1
243 movs xl, xl, rrx
244 movs ip, ip, rrx
245 orrcs ip, ip, #1
246 add r4, r4, #(1 << 20)
248 @ Our result is now properly aligned into xh-xl, remaining bits in ip.
249 @ Round with MSB of ip. If halfway between two numbers, round towards
250 @ LSB of xl = 0.
251 LSYM(Lad_r0):
252 adds xl, xl, ip, lsr #31
253 adc xh, xh, #0
254 teq ip, #0x80000000
255 biceq xl, xl, #1
257 @ One extreme rounding case may add a new MSB. Adjust exponent.
258 @ That MSB will be cleared when exponent is merged below.
259 tst xh, #0x00200000
260 addne r4, r4, #(1 << 20)
262 @ Make sure we did not bust our exponent.
263 adds ip, r4, #(1 << 20)
264 bmi LSYM(Lad_o)
266 @ Pack final result together.
267 LSYM(Lad_e):
268 bic xh, xh, #0x00300000
269 orr xh, xh, r4
270 orr xh, xh, r5
271 #ifdef __FP_INTERWORKING__
272 ldmfd sp!, {r4, r5, lr}
273 bx lr
274 #else
275 ldmfd sp!, {r4, r5, pc}RETCOND
276 #endif
278 LSYM(Lad_l): @ Result must be shifted left and exponent adjusted.
279 @ No rounding necessary since ip will always be 0.
280 #if __ARM_ARCH__ < 5
282 teq xh, #0
283 movne r3, #-11
284 moveq r3, #21
285 moveq xh, xl
286 moveq xl, #0
287 mov r2, xh
288 movs ip, xh, lsr #16
289 moveq r2, r2, lsl #16
290 addeq r3, r3, #16
291 tst r2, #0xff000000
292 moveq r2, r2, lsl #8
293 addeq r3, r3, #8
294 tst r2, #0xf0000000
295 moveq r2, r2, lsl #4
296 addeq r3, r3, #4
297 tst r2, #0xc0000000
298 moveq r2, r2, lsl #2
299 addeq r3, r3, #2
300 tst r2, #0x80000000
301 addeq r3, r3, #1
303 #else
305 teq xh, #0
306 moveq xh, xl
307 moveq xl, #0
308 clz r3, xh
309 addeq r3, r3, #32
310 sub r3, r3, #11
312 #endif
314 @ determine how to shift the value.
315 subs r2, r3, #32
316 bge 2f
317 adds r2, r2, #12
318 ble 1f
320 @ shift value left 21 to 31 bits, or actually right 11 to 1 bits
321 @ since a register switch happened above.
322 add ip, r2, #20
323 rsb r2, r2, #12
324 mov xl, xh, lsl ip
325 mov xh, xh, lsr r2
326 b 3f
328 @ actually shift value left 1 to 20 bits, which might also represent
329 @ 32 to 52 bits if counting the register switch that happened earlier.
330 1: add r2, r2, #20
331 2: rsble ip, r2, #32
332 mov xh, xh, lsl r2
333 orrle xh, xh, xl, lsr ip
334 movle xl, xl, lsl r2
336 @ adjust exponent accordingly.
337 3: subs r4, r4, r3, lsl #20
338 bgt LSYM(Lad_e)
340 @ Exponent too small, denormalize result.
341 @ Find out proper shift value.
342 mvn r4, r4, asr #20
343 subs r4, r4, #30
344 bge 2f
345 adds r4, r4, #12
346 bgt 1f
348 @ shift result right of 1 to 20 bits, sign is in r5.
349 add r4, r4, #20
350 rsb r2, r4, #32
351 mov xl, xl, lsr r4
352 orr xl, xl, xh, lsl r2
353 orr xh, r5, xh, lsr r4
354 #ifdef __FP_INTERWORKING
355 ldmfd sp!, {r4, r5, lr}
356 bx lr
357 #else
358 ldmfd sp!, {r4, r5, pc}RETCOND
359 #endif
361 @ shift result right of 21 to 31 bits, or left 11 to 1 bits after
362 @ a register switch from xh to xl.
363 1: rsb r4, r4, #12
364 rsb r2, r4, #32
365 mov xl, xl, lsr r2
366 orr xl, xl, xh, lsl r4
367 mov xh, r5
368 #ifdef __FP_INTERWORKING__
369 ldmfd sp!, {r4, r5, lr}
370 bx lr
371 #else
372 ldmfd sp!, {r4, r5, pc}RETCOND
373 #endif
375 @ Shift value right of 32 to 64 bits, or 0 to 32 bits after a switch
376 @ from xh to xl.
377 2: mov xl, xh, lsr r4
378 mov xh, r5
379 #ifdef __FP_INTERWORKING__
380 ldmfd sp!, {r4, r5, lr}
381 bx lr
382 #else
383 ldmfd sp!, {r4, r5, pc}RETCOND
384 #endif
386 @ Adjust exponents for denormalized arguments.
387 LSYM(Lad_d):
388 teq r4, #0
389 eoreq xh, xh, #0x00100000
390 addeq r4, r4, #(1 << 20)
391 eor yh, yh, #0x00100000
392 subne r5, r5, #(1 << 20)
393 b LSYM(Lad_x)
395 @ Result is x - x = 0, unless x = INF or NAN.
396 LSYM(Lad_z):
397 sub ip, ip, #0x00100000 @ ip becomes 0x7ff00000
398 and r2, xh, ip
399 teq r2, ip
400 orreq xh, ip, #0x00080000
401 movne xh, #0
402 mov xl, #0
403 RET
405 @ Overflow: return INF.
406 LSYM(Lad_o):
407 orr xh, r5, #0x7f000000
408 orr xh, xh, #0x00f00000
409 mov xl, #0
410 #ifdef __FP_INTERWORKING__
411 ldmfd sp!, {r4, r5, lr}
412 bx lr
413 #else
414 ldmfd sp!, {r4, r5, pc}RETCOND
415 #endif
417 @ At least one of x or y is INF/NAN.
418 @ if xh-xl != INF/NAN: return yh-yl (which is INF/NAN)
419 @ if yh-yl != INF/NAN: return xh-xl (which is INF/NAN)
420 @ if either is NAN: return NAN
421 @ if opposite sign: return NAN
422 @ return xh-xl (which is INF or -INF)
423 LSYM(Lad_i):
424 teq r4, ip
425 movne xh, yh
426 movne xl, yl
427 teqeq r5, ip
428 #ifdef __FP_INTERWORKING__
429 ldmnefd sp!, {r4, r5, lr}
430 bxne lr
431 #else
432 ldmnefd sp!, {r4, r5, pc}RETCOND
433 #endif
434 orrs r4, xl, xh, lsl #12
435 orreqs r4, yl, yh, lsl #12
436 teqeq xh, yh
437 orrne xh, r5, #0x00080000
438 movne xl, #0
439 #ifdef __FP_INTERWORKING__
440 ldmfd sp!, {r4, r5, lr}
441 bx lr
442 #else
443 ldmfd sp!, {r4, r5, pc}RETCOND
444 #endif
447 ARM_FUNC_START floatunsidf
448 teq r0, #0
449 moveq r1, #0
450 RETc(eq)
451 stmfd sp!, {r4, r5, lr}
452 mov r4, #(0x400 << 20) @ initial exponent
453 add r4, r4, #((52-1) << 20)
454 mov r5, #0 @ sign bit is 0
455 mov xl, r0
456 mov xh, #0
457 b LSYM(Lad_l)
460 ARM_FUNC_START floatsidf
461 teq r0, #0
462 moveq r1, #0
463 RETc(eq)
464 stmfd sp!, {r4, r5, lr}
465 mov r4, #(0x400 << 20) @ initial exponent
466 add r4, r4, #((52-1) << 20)
467 ands r5, r0, #0x80000000 @ sign bit in r5
468 rsbmi r0, r0, #0 @ absolute value
469 mov xl, r0
470 mov xh, #0
471 b LSYM(Lad_l)
474 ARM_FUNC_START extendsfdf2
475 movs r2, r0, lsl #1
476 beq 1f @ value is 0.0 or -0.0
477 mov xh, r2, asr #3 @ stretch exponent
478 mov xh, xh, rrx @ retrieve sign bit
479 mov xl, r2, lsl #28 @ retrieve remaining bits
480 ands r2, r2, #0xff000000 @ isolate exponent
481 beq 2f @ exponent was 0 but not mantissa
482 teq r2, #0xff000000 @ check if INF or NAN
483 eorne xh, xh, #0x38000000 @ fixup exponent otherwise.
484 RET
486 1: mov xh, r0
487 mov xl, #0
488 RET
490 2: @ value was denormalized. We can normalize it now.
491 stmfd sp!, {r4, r5, lr}
492 mov r4, #(0x380 << 20) @ setup corresponding exponent
493 add r4, r4, #(1 << 20)
494 and r5, xh, #0x80000000 @ move sign bit in r5
495 bic xh, xh, #0x80000000
496 b LSYM(Lad_l)
499 ARM_FUNC_START muldf3
501 stmfd sp!, {r4, r5, r6, lr}
503 @ Mask out exponents.
504 mov ip, #0x7f000000
505 orr ip, ip, #0x00f00000
506 and r4, xh, ip
507 and r5, yh, ip
509 @ Trap any INF/NAN.
510 teq r4, ip
511 teqne r5, ip
512 beq LSYM(Lml_s)
514 @ Trap any multiplication by 0.
515 orrs r6, xl, xh, lsl #1
516 orrnes r6, yl, yh, lsl #1
517 beq LSYM(Lml_z)
519 @ Shift exponents right one bit to make room for overflow bit.
520 @ If either of them is 0, scale denormalized arguments off line.
521 @ Then add both exponents together.
522 movs r4, r4, lsr #1
523 teqne r5, #0
524 beq LSYM(Lml_d)
525 LSYM(Lml_x):
526 add r4, r4, r5, asr #1
528 @ Preserve final sign in r4 along with exponent for now.
529 teq xh, yh
530 orrmi r4, r4, #0x8000
532 @ Convert mantissa to unsigned integer.
533 bic xh, xh, ip, lsl #1
534 bic yh, yh, ip, lsl #1
535 orr xh, xh, #0x00100000
536 orr yh, yh, #0x00100000
538 #if __ARM_ARCH__ < 4
540 @ Well, no way to make it shorter without the umull instruction.
541 @ We must perform that 53 x 53 bit multiplication by hand.
542 stmfd sp!, {r7, r8, r9, sl, fp}
543 mov r7, xl, lsr #16
544 mov r8, yl, lsr #16
545 mov r9, xh, lsr #16
546 mov sl, yh, lsr #16
547 bic xl, xl, r7, lsl #16
548 bic yl, yl, r8, lsl #16
549 bic xh, xh, r9, lsl #16
550 bic yh, yh, sl, lsl #16
551 mul ip, xl, yl
552 mul fp, xl, r8
553 mov lr, #0
554 adds ip, ip, fp, lsl #16
555 adc lr, lr, fp, lsr #16
556 mul fp, r7, yl
557 adds ip, ip, fp, lsl #16
558 adc lr, lr, fp, lsr #16
559 mul fp, xl, sl
560 mov r5, #0
561 adds lr, lr, fp, lsl #16
562 adc r5, r5, fp, lsr #16
563 mul fp, r7, yh
564 adds lr, lr, fp, lsl #16
565 adc r5, r5, fp, lsr #16
566 mul fp, xh, r8
567 adds lr, lr, fp, lsl #16
568 adc r5, r5, fp, lsr #16
569 mul fp, r9, yl
570 adds lr, lr, fp, lsl #16
571 adc r5, r5, fp, lsr #16
572 mul fp, xh, sl
573 mul r6, r9, sl
574 adds r5, r5, fp, lsl #16
575 adc r6, r6, fp, lsr #16
576 mul fp, r9, yh
577 adds r5, r5, fp, lsl #16
578 adc r6, r6, fp, lsr #16
579 mul fp, xl, yh
580 adds lr, lr, fp
581 mul fp, r7, sl
582 adcs r5, r5, fp
583 mul fp, xh, yl
584 adc r6, r6, #0
585 adds lr, lr, fp
586 mul fp, r9, r8
587 adcs r5, r5, fp
588 mul fp, r7, r8
589 adc r6, r6, #0
590 adds lr, lr, fp
591 mul fp, xh, yh
592 adcs r5, r5, fp
593 adc r6, r6, #0
594 ldmfd sp!, {r7, r8, r9, sl, fp}
596 #else
598 @ Here is the actual multiplication: 53 bits * 53 bits -> 106 bits.
599 umull ip, lr, xl, yl
600 mov r5, #0
601 umlal lr, r5, xl, yh
602 umlal lr, r5, xh, yl
603 mov r6, #0
604 umlal r5, r6, xh, yh
606 #endif
608 @ The LSBs in ip are only significant for the final rounding.
609 @ Fold them into one bit of lr.
610 teq ip, #0
611 orrne lr, lr, #1
613 @ Put final sign in xh.
614 mov xh, r4, lsl #16
615 bic r4, r4, #0x8000
617 @ Adjust result if one extra MSB appeared (one of four times).
618 tst r6, #(1 << 9)
619 beq 1f
620 add r4, r4, #(1 << 19)
621 movs r6, r6, lsr #1
622 movs r5, r5, rrx
623 movs lr, lr, rrx
624 orrcs lr, lr, #1
625 1:
626 @ Scale back to 53 bits.
627 @ xh contains sign bit already.
628 orr xh, xh, r6, lsl #12
629 orr xh, xh, r5, lsr #20
630 mov xl, r5, lsl #12
631 orr xl, xl, lr, lsr #20
633 @ Apply exponent bias, check range for underflow.
634 sub r4, r4, #0x00f80000
635 subs r4, r4, #0x1f000000
636 ble LSYM(Lml_u)
638 @ Round the result.
639 movs lr, lr, lsl #12
640 bpl 1f
641 adds xl, xl, #1
642 adc xh, xh, #0
643 teq lr, #0x80000000
644 biceq xl, xl, #1
646 @ Rounding may have produced an extra MSB here.
647 @ The extra bit is cleared before merging the exponent below.
648 tst xh, #0x00200000
649 addne r4, r4, #(1 << 19)
650 1:
651 @ Check exponent for overflow.
652 adds ip, r4, #(1 << 19)
653 tst ip, #(1 << 30)
654 bne LSYM(Lml_o)
656 @ Add final exponent.
657 bic xh, xh, #0x00300000
658 orr xh, xh, r4, lsl #1
659 #ifdef __FP_INTERWORKING__
660 ldmfd sp!, {r4, r5, r6, lr}
661 bx lr
662 #else
663 ldmfd sp!, {r4, r5, r6, pc}RETCOND
664 #endif
666 @ Result is 0, but determine sign anyway.
667 LSYM(Lml_z):
668 eor xh, xh, yh
669 LSYM(Ldv_z):
670 bic xh, xh, #0x7fffffff
671 mov xl, #0
672 #ifdef __FP_INTERWORKING__
673 ldmfd sp!, {r4, r5, r6, lr}
674 bx lr
675 #else
676 ldmfd sp!, {r4, r5, r6, pc}RETCOND
677 #endif
679 @ Check if denormalized result is possible, otherwise return signed 0.
680 LSYM(Lml_u):
681 cmn r4, #(53 << 19)
682 movle xl, #0
683 #ifdef __FP_INTERWORKING__
684 ldmlefd sp!, {r4, r5, r6, lr}
685 bxle lr
686 #else
687 ldmlefd sp!, {r4, r5, r6, pc}RETCOND
688 #endif
690 @ Find out proper shift value.
691 LSYM(Lml_r):
692 mvn r4, r4, asr #19
693 subs r4, r4, #30
694 bge 2f
695 adds r4, r4, #12
696 bgt 1f
698 @ shift result right of 1 to 20 bits, preserve sign bit, round, etc.
699 add r4, r4, #20
700 rsb r5, r4, #32
701 mov r3, xl, lsl r5
702 mov xl, xl, lsr r4
703 orr xl, xl, xh, lsl r5
704 movs xh, xh, lsl #1
705 mov xh, xh, lsr r4
706 mov xh, xh, rrx
707 adds xl, xl, r3, lsr #31
708 adc xh, xh, #0
709 teq lr, #0
710 teqeq r3, #0x80000000
711 biceq xl, xl, #1
712 #ifdef __FP_INTERWORKING__
713 ldmfd sp!, {r4, r5, r6, lr}
714 bx lr
715 #else
716 ldmfd sp!, {r4, r5, r6, pc}RETCOND
717 #endif
719 @ shift result right of 21 to 31 bits, or left 11 to 1 bits after
720 @ a register switch from xh to xl. Then round.
721 1: rsb r4, r4, #12
722 rsb r5, r4, #32
723 mov r3, xl, lsl r4
724 mov xl, xl, lsr r5
725 orr xl, xl, xh, lsl r4
726 bic xh, xh, #0x7fffffff
727 adds xl, xl, r3, lsr #31
728 adc xh, xh, #0
729 teq lr, #0
730 teqeq r3, #0x80000000
731 biceq xl, xl, #1
732 #ifdef __FP_INTERWORKING__
733 ldmfd sp!, {r4, r5, r6, lr}
734 bx lr
735 #else
736 ldmfd sp!, {r4, r5, r6, pc}RETCOND
737 #endif
739 @ Shift value right of 32 to 64 bits, or 0 to 32 bits after a switch
740 @ from xh to xl. Leftover bits are in r3-r6-lr for rounding.
741 2: rsb r5, r4, #32
742 mov r6, xl, lsl r5
743 mov r3, xl, lsr r4
744 orr r3, r3, xh, lsl r5
745 mov xl, xh, lsr r4
746 bic xh, xh, #0x7fffffff
747 adds xl, xl, r3, lsr #31
748 adc xh, xh, #0
749 orrs r6, r6, lr
750 teqeq r3, #0x80000000
751 biceq xl, xl, #1
752 #ifdef __FP_INTERWORKING__
753 ldmfd sp!, {r4, r5, r6, lr}
754 bx lr
755 #else
756 ldmfd sp!, {r4, r5, r6, pc}RETCOND
757 #endif
759 @ One or both arguments are denormalized.
760 @ Scale them leftwards and preserve sign bit.
761 LSYM(Lml_d):
762 mov lr, #0
763 teq r4, #0
764 bne 2f
765 and r6, xh, #0x80000000
766 1: movs xl, xl, lsl #1
767 adc xh, lr, xh, lsl #1
768 tst xh, #0x00100000
769 subeq r4, r4, #(1 << 19)
770 beq 1b
771 orr xh, xh, r6
772 teq r5, #0
773 bne LSYM(Lml_x)
774 2: and r6, yh, #0x80000000
775 3: movs yl, yl, lsl #1
776 adc yh, lr, yh, lsl #1
777 tst yh, #0x00100000
778 subeq r5, r5, #(1 << 20)
779 beq 3b
780 orr yh, yh, r6
781 b LSYM(Lml_x)
783 @ One or both args are INF or NAN.
784 LSYM(Lml_s):
785 orrs r6, xl, xh, lsl #1
786 orrnes r6, yl, yh, lsl #1
787 beq LSYM(Lml_n) @ 0 * INF or INF * 0 -> NAN
788 teq r4, ip
789 bne 1f
790 orrs r6, xl, xh, lsl #12
791 bne LSYM(Lml_n) @ NAN * <anything> -> NAN
792 1: teq r5, ip
793 bne LSYM(Lml_i)
794 orrs r6, yl, yh, lsl #12
795 bne LSYM(Lml_n) @ <anything> * NAN -> NAN
797 @ Result is INF, but we need to determine its sign.
798 LSYM(Lml_i):
799 eor xh, xh, yh
801 @ Overflow: return INF (sign already in xh).
802 LSYM(Lml_o):
803 and xh, xh, #0x80000000
804 orr xh, xh, #0x7f000000
805 orr xh, xh, #0x00f00000
806 mov xl, #0
807 #ifdef __FP_INTERWORKING__
808 ldmfd sp!, {r4, r5, r6, lr}
809 bx lr
810 #else
811 ldmfd sp!, {r4, r5, r6, pc}RETCOND
812 #endif
814 @ Return NAN.
815 LSYM(Lml_n):
816 mov xh, #0x7f000000
817 orr xh, xh, #0x00f80000
818 #ifdef __FP_INTERWORKING__
819 ldmfd sp!, {r4, r5, r6, lr}
820 bx lr
821 #else
822 ldmfd sp!, {r4, r5, r6, pc}RETCOND
823 #endif
826 ARM_FUNC_START divdf3
828 stmfd sp!, {r4, r5, r6, lr}
830 @ Mask out exponents.
831 mov ip, #0x7f000000
832 orr ip, ip, #0x00f00000
833 and r4, xh, ip
834 and r5, yh, ip
836 @ Trap any INF/NAN or zeroes.
837 teq r4, ip
838 teqne r5, ip
839 orrnes r6, xl, xh, lsl #1
840 orrnes r6, yl, yh, lsl #1
841 beq LSYM(Ldv_s)
843 @ Shift exponents right one bit to make room for overflow bit.
844 @ If either of them is 0, scale denormalized arguments off line.
845 @ Then substract divisor exponent from dividend''s.
846 movs r4, r4, lsr #1
847 teqne r5, #0
848 beq LSYM(Ldv_d)
849 LSYM(Ldv_x):
850 sub r4, r4, r5, asr #1
852 @ Preserve final sign into lr.
853 eor lr, xh, yh
855 @ Convert mantissa to unsigned integer.
856 @ Dividend -> r5-r6, divisor -> yh-yl.
857 mov r5, #0x10000000
858 mov yh, yh, lsl #12
859 orr yh, r5, yh, lsr #4
860 orr yh, yh, yl, lsr #24
861 movs yl, yl, lsl #8
862 mov xh, xh, lsl #12
863 teqeq yh, r5
864 beq LSYM(Ldv_1)
865 orr r5, r5, xh, lsr #4
866 orr r5, r5, xl, lsr #24
867 mov r6, xl, lsl #8
869 @ Initialize xh with final sign bit.
870 and xh, lr, #0x80000000
872 @ Ensure result will land to known bit position.
873 cmp r5, yh
874 cmpeq r6, yl
875 bcs 1f
876 sub r4, r4, #(1 << 19)
877 movs yh, yh, lsr #1
878 mov yl, yl, rrx
879 1:
880 @ Apply exponent bias, check range for over/underflow.
881 add r4, r4, #0x1f000000
882 add r4, r4, #0x00f80000
883 cmn r4, #(53 << 19)
884 ble LSYM(Ldv_z)
885 cmp r4, ip, lsr #1
886 bge LSYM(Lml_o)
888 @ Perform first substraction to align result to a nibble.
889 subs r6, r6, yl
890 sbc r5, r5, yh
891 movs yh, yh, lsr #1
892 mov yl, yl, rrx
893 mov xl, #0x00100000
894 mov ip, #0x00080000
896 @ The actual division loop.
897 1: subs lr, r6, yl
898 sbcs lr, r5, yh
899 subcs r6, r6, yl
900 movcs r5, lr
901 orrcs xl, xl, ip
902 movs yh, yh, lsr #1
903 mov yl, yl, rrx
904 subs lr, r6, yl
905 sbcs lr, r5, yh
906 subcs r6, r6, yl
907 movcs r5, lr
908 orrcs xl, xl, ip, lsr #1
909 movs yh, yh, lsr #1
910 mov yl, yl, rrx
911 subs lr, r6, yl
912 sbcs lr, r5, yh
913 subcs r6, r6, yl
914 movcs r5, lr
915 orrcs xl, xl, ip, lsr #2
916 movs yh, yh, lsr #1
917 mov yl, yl, rrx
918 subs lr, r6, yl
919 sbcs lr, r5, yh
920 subcs r6, r6, yl
921 movcs r5, lr
922 orrcs xl, xl, ip, lsr #3
924 orrs lr, r5, r6
925 beq 2f
926 mov r5, r5, lsl #4
927 orr r5, r5, r6, lsr #28
928 mov r6, r6, lsl #4
929 mov yh, yh, lsl #3
930 orr yh, yh, yl, lsr #29
931 mov yl, yl, lsl #3
932 movs ip, ip, lsr #4
933 bne 1b
935 @ We are done with a word of the result.
936 @ Loop again for the low word if this pass was for the high word.
937 tst xh, #0x00100000
938 bne 3f
939 orr xh, xh, xl
940 mov xl, #0
941 mov ip, #0x80000000
942 b 1b
943 2:
944 @ Be sure result starts in the high word.
945 tst xh, #0x00100000
946 orreq xh, xh, xl
947 moveq xl, #0
948 3:
949 @ Check if denormalized result is needed.
950 cmp r4, #0
951 ble LSYM(Ldv_u)
953 @ Apply proper rounding.
954 subs ip, r5, yh
955 subeqs ip, r6, yl
956 adcs xl, xl, #0
957 adc xh, xh, #0
958 teq ip, #0
959 biceq xl, xl, #1
961 @ Add exponent to result.
962 bic xh, xh, #0x00100000
963 orr xh, xh, r4, lsl #1
964 #ifdef __FP_INTERWORKING__
965 ldmfd sp!, {r4, r5, r6, lr}
966 bx lr
967 #else
968 ldmfd sp!, {r4, r5, r6, pc}RETCOND
969 #endif
971 @ Division by 0x1p*: shortcut a lot of code.
972 LSYM(Ldv_1):
973 and lr, lr, #0x80000000
974 orr xh, lr, xh, lsr #12
975 add r4, r4, #0x1f000000
976 add r4, r4, #0x00f80000
977 cmp r4, ip, lsr #1
978 bge LSYM(Lml_o)
979 cmp r4, #0
980 orrgt xh, xh, r4, lsl #1
981 #ifdef __FP_INTERWORKING__
982 ldmgtfd sp!, {r4, r5, r6, lr}
983 bxgt lr
984 #else
985 ldmgtfd sp!, {r4, r5, r6, pc}RETCOND
986 #endif
987 cmn r4, #(53 << 19)
988 ble LSYM(Ldv_z)
989 orr xh, xh, #0x00100000
990 mov lr, #0
991 b LSYM(Lml_r)
993 @ Result must be denormalized: put remainder in lr for
994 @ rounding considerations.
995 LSYM(Ldv_u):
996 orr lr, r5, r6
997 b LSYM(Lml_r)
999 @ One or both arguments are denormalized.
1000 @ Scale them leftwards and preserve sign bit.
1001 LSYM(Ldv_d):
1002 mov lr, #0
1003 teq r4, #0
1004 bne 2f
1005 and r6, xh, #0x80000000
1006 1: movs xl, xl, lsl #1
1007 adc xh, lr, xh, lsl #1
1008 tst xh, #0x00100000
1009 subeq r4, r4, #(1 << 19)
1010 beq 1b
1011 orr xh, xh, r6
1012 teq r5, #0
1013 bne LSYM(Ldv_x)
1014 2: and r6, yh, #0x80000000
1015 3: movs yl, yl, lsl #1
1016 adc yh, lr, yh, lsl #1
1017 tst yh, #0x00100000
1018 subeq r5, r5, #(1 << 20)
1019 beq 3b
1020 orr yh, yh, r6
1021 b LSYM(Ldv_x)
1023 @ One or both arguments is either INF, NAN or zero.
1024 LSYM(Ldv_s):
1025 teq r4, ip
1026 teqeq r5, ip
1027 beq LSYM(Lml_n) @ INF/NAN / INF/NAN -> NAN
1028 teq r4, ip
1029 bne 1f
1030 orrs r4, xl, xh, lsl #12
1031 bne LSYM(Lml_n) @ NAN / <anything> -> NAN
1032 b LSYM(Lml_i) @ INF / <anything> -> INF
1033 1: teq r5, ip
1034 bne 2f
1035 orrs r5, yl, yh, lsl #12
1036 bne LSYM(Lml_n) @ <anything> / NAN -> NAN
1037 b LSYM(Lml_z) @ <anything> / INF -> 0
1038 2: @ One or both arguments are 0.
1039 orrs r4, xl, xh, lsl #1
1040 bne LSYM(Lml_i) @ <non_zero> / 0 -> INF
1041 orrs r5, yl, yh, lsl #1
1042 bne LSYM(Lml_z) @ 0 / <non_zero> -> 0
1043 b LSYM(Lml_n) @ 0 / 0 -> NAN
1046 FUNC_START gedf2
1047 ARM_FUNC_START gtdf2
1048 mov ip, #-1
1049 b 1f
1051 FUNC_START ledf2
1052 ARM_FUNC_START ltdf2
1053 mov ip, #1
1054 b 1f
1056 FUNC_START nedf2
1057 FUNC_START eqdf2
1058 ARM_FUNC_START cmpdf2
1059 mov ip, #1 @ how should we specify unordered here?
1061 1: stmfd sp!, {r4, r5, lr}
1063 @ Trap any INF/NAN first.
1064 mov lr, #0x7f000000
1065 orr lr, lr, #0x00f00000
1066 and r4, xh, lr
1067 and r5, yh, lr
1068 teq r4, lr
1069 teqne r5, lr
1070 beq 3f
1072 @ Test for equality.
1073 @ Note that 0.0 is equal to -0.0.
1074 2: orrs ip, xl, xh, lsl #1 @ if x == 0.0 or -0.0
1075 orreqs ip, yl, yh, lsl #1 @ and y == 0.0 or -0.0
1076 teqne xh, yh @ or xh == yh
1077 teqeq xl, yl @ and xl == yl
1078 moveq r0, #0 @ then equal.
1079 #ifdef __FP_INTERWORKING__
1080 ldmeqfd sp!, {r4, r5, lr}
1081 bxeq lr
1082 #else
1083 ldmeqfd sp!, {r4, r5, pc}RETCOND
1084 #endif
1086 @ Check for sign difference.
1087 teq xh, yh
1088 movmi r0, xh, asr #31
1089 orrmi r0, r0, #1
1090 #ifdef __FP_INTERWORKING__
1091 ldmmifd sp!, {r4, r5, lr}
1092 bxmi lr
1093 #else
1094 ldmmifd sp!, {r4, r5, pc}RETCOND
1095 #endif
1097 @ Compare exponents.
1098 cmp r4, r5
1100 @ Compare mantissa if exponents are equal.
1101 moveq xh, xh, lsl #12
1102 cmpeq xh, yh, lsl #12
1103 cmpeq xl, yl
1104 movcs r0, yh, asr #31
1105 mvncc r0, yh, asr #31
1106 orr r0, r0, #1
1107 #ifdef __FP_INTERWORKING__
1108 ldmfd sp!, {r4, r5, lr}
1109 bx lr
1110 #else
1111 ldmfd sp!, {r4, r5, pc}RETCOND
1112 #endif
1114 @ Look for a NAN.
1115 3: teq r4, lr
1116 bne 4f
1117 orrs xl, xl, xh, lsl #12
1118 bne 5f @ x is NAN
1119 4: teq r5, lr
1120 bne 2b
1121 orrs yl, yl, yh, lsl #12
1122 beq 2b @ y is not NAN
1123 5: mov r0, ip @ return unordered code from ip
1124 #ifdef __FP_INTERWORKING__
1125 ldmfd sp!, {r4, r5, lr}
1126 bx lr
1127 #else
1128 ldmfd sp!, {r4, r5, pc}RETCOND
1129 #endif
1132 ARM_FUNC_START unorddf2
1133 str lr, [sp, #-4]!
1134 mov ip, #0x7f000000
1135 orr ip, ip, #0x00f00000
1136 and lr, xh, ip
1137 teq lr, ip
1138 bne 1f
1139 orrs xl, xl, xh, lsl #12
1140 bne 3f @ x is NAN
1141 1: and lr, yh, ip
1142 teq lr, ip
1143 bne 2f
1144 orrs yl, yl, yh, lsl #12
1145 bne 3f @ y is NAN
1146 2: mov r0, #0 @ arguments are ordered.
1147 #ifdef __FP_INTERWORKING__
1148 ldr lr, [sp], #4
1149 bx lr
1150 #elif defined (__APCS_26__)
1151 ldmia sp!, {pc}^
1152 #else
1153 ldr pc, [sp], #4
1154 #endif
1155 3: mov r0, #1 @ arguments are unordered.
1156 #ifdef __FP_INTERWORKING__
1157 ldr lr, [sp], #4
1158 bx lr
1159 #elif defined (__APCS_26__)
1160 ldmia sp!, {pc}^
1161 #else
1162 ldr pc, [sp], #4
1163 #endif
1166 ARM_FUNC_START fixdfsi
1167 orrs ip, xl, xh, lsl #1
1168 beq 1f @ value is 0.
1170 @ preserve C flag (the actual sign)
1171 #ifdef __APCS_26__
1172 mov r3, pc
1173 #else
1174 mrs r3, cpsr
1175 #endif
1177 @ check exponent range.
1178 mov ip, #0x7f000000
1179 orr ip, ip, #0x00f00000
1180 and r2, xh, ip
1181 teq r2, ip
1182 beq 2f @ value is INF or NAN
1183 bic ip, ip, #0x40000000
1184 cmp r2, ip
1185 bcc 1f @ value is too small
1186 add ip, ip, #(31 << 20)
1187 cmp r2, ip
1188 bcs 3f @ value is too large
1190 rsb r2, r2, ip
1191 mov ip, xh, lsl #11
1192 orr ip, ip, #0x80000000
1193 orr ip, ip, xl, lsr #21
1194 mov r2, r2, lsr #20
1195 mov r0, ip, lsr r2
1196 tst r3, #0x20000000 @ the sign bit
1197 rsbne r0, r0, #0
1198 RET
1200 1: mov r0, #0
1201 RET
1203 2: orrs xl, xl, xh, lsl #12
1204 bne 4f @ r0 is NAN.
1205 3: tst r3, #0x20000000 @ the sign bit
1206 moveq r0, #0x7fffffff @ maximum signed positive si
1207 movne r0, #0x80000000 @ maximum signed negative si
1208 RET
1210 4: mov r0, #0 @ How should we convert NAN?
1211 RET
1213 ARM_FUNC_START fixunsdfsi
1214 orrs ip, xl, xh, lsl #1
1215 beq 1b @ value is 0
1216 bcs 1b @ value is negative
1218 @ check exponent range.
1219 mov ip, #0x7f000000
1220 orr ip, ip, #0x00f00000
1221 and r2, xh, ip
1222 teq r2, ip
1223 beq 1f @ value is INF or NAN
1224 bic ip, ip, #0x40000000
1225 cmp r2, ip
1226 bcc 1b @ value is too small
1227 add ip, ip, #(31 << 20)
1228 cmp r2, ip
1229 bhi 2f @ value is too large
1231 rsb r2, r2, ip
1232 mov ip, xh, lsl #11
1233 orr ip, ip, #0x80000000
1234 orr ip, ip, xl, lsr #21
1235 mov r2, r2, lsr #20
1236 mov r0, ip, lsr r2
1237 RET
1239 1: orrs xl, xl, xh, lsl #12
1240 bne 4b @ value is NAN.
1241 2: mov r0, #0xffffffff @ maximum unsigned si
1242 RET
1245 ARM_FUNC_START truncdfsf2
1246 orrs r2, xl, xh, lsl #1
1247 moveq r0, r2, rrx
1248 RETc(eq) @ value is 0.0 or -0.0
1250 @ check exponent range.
1251 mov ip, #0x7f000000
1252 orr ip, ip, #0x00f00000
1253 and r2, ip, xh
1254 teq r2, ip
1255 beq 2f @ value is INF or NAN
1256 bic xh, xh, ip
1257 cmp r2, #(0x380 << 20)
1258 bls 4f @ value is too small
1260 @ shift and round mantissa
1261 1: movs r3, xl, lsr #29
1262 adc r3, r3, xh, lsl #3
1264 @ if halfway between two numbers, round towards LSB = 0.
1265 mov xl, xl, lsl #3
1266 teq xl, #0x80000000
1267 biceq r3, r3, #1
1269 @ rounding might have created an extra MSB. If so adjust exponent.
1270 tst r3, #0x00800000
1271 addne r2, r2, #(1 << 20)
1272 bicne r3, r3, #0x00800000
1274 @ check exponent for overflow
1275 mov ip, #(0x400 << 20)
1276 orr ip, ip, #(0x07f << 20)
1277 cmp r2, ip
1278 bcs 3f @ overflow
1280 @ adjust exponent, merge with sign bit and mantissa.
1281 movs xh, xh, lsl #1
1282 mov r2, r2, lsl #4
1283 orr r0, r3, r2, rrx
1284 eor r0, r0, #0x40000000
1285 RET
1287 2: @ chech for NAN
1288 orrs xl, xl, xh, lsl #12
1289 movne r0, #0x7f000000
1290 orrne r0, r0, #0x00c00000
1291 RETc(ne) @ return NAN
1293 3: @ return INF with sign
1294 and r0, xh, #0x80000000
1295 orr r0, r0, #0x7f000000
1296 orr r0, r0, #0x00800000
1297 RET
1299 4: @ check if denormalized value is possible
1300 subs r2, r2, #((0x380 - 24) << 20)
1301 andle r0, xh, #0x80000000 @ too small, return signed 0.
1302 RETc(le)
1304 @ denormalize value so we can resume with the code above afterwards.
1305 orr xh, xh, #0x00100000
1306 mov r2, r2, lsr #20
1307 rsb r2, r2, #25
1308 cmp r2, #20
1309 bgt 6f
1311 rsb ip, r2, #32
1312 mov r3, xl, lsl ip
1313 mov xl, xl, lsr r2
1314 orr xl, xl, xh, lsl ip
1315 movs xh, xh, lsl #1
1316 mov xh, xh, lsr r2
1317 mov xh, xh, rrx
1318 5: teq r3, #0 @ fold r3 bits into the LSB
1319 orrne xl, xl, #1 @ for rounding considerations.
1320 mov r2, #(0x380 << 20) @ equivalent to the 0 float exponent
1321 b 1b
1323 6: rsb r2, r2, #(12 + 20)
1324 rsb ip, r2, #32
1325 mov r3, xl, lsl r2
1326 mov xl, xl, lsr ip
1327 orr xl, xl, xh, lsl r2
1328 and xh, xh, #0x80000000
1329 b 5b