1 /* libgcc routines for
68000 w
/o floating
-point hardware.
2 Copyright
(C
) 1994, 1996, 1997, 1998, 2008, 2009 Free Software Foundation
, Inc.
4 This file is part of GCC.
6 GCC 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 3, or (at your option) any
11 This file is distributed
in the hope that it will be useful
, but
12 WITHOUT ANY WARRANTY
; without even the implied warranty of
13 MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 General
Public License for more details.
16 Under
Section 7 of GPL version
3, you are granted additional
17 permissions described
in the GCC Runtime Library Exception
, version
18 3.1, as published by the Free Software Foundation.
20 You should have received a copy of the GNU General
Public License
and
21 a copy of the GCC Runtime Library Exception along with
this program
;
22 see the files COPYING3
and COPYING.RUNTIME respectively. If
not, see
23 <http://www.gnu.
org/licenses
/>.
*/
25 /* Use
this one for any
680x0
; assumes no floating point hardware.
26 The trailing
" '" appearing on some lines is for ANSI preprocessors. Yuk.
27 Some of
this code comes from MINIX
, via the folks at ericsson.
28 D. V. Henkel
-Wallace
(gumby
@cygnus.com
) Fete Bastille
, 1992
31 /* These are predefined by new versions of GNU cpp.
*/
33 #ifndef __USER_LABEL_PREFIX__
34 #define __USER_LABEL_PREFIX__ _
37 #ifndef __REGISTER_PREFIX__
38 #define __REGISTER_PREFIX__
41 #ifndef __IMMEDIATE_PREFIX__
42 #define __IMMEDIATE_PREFIX__ #
45 /* ANSI concatenation macros.
*/
47 #define CONCAT1
(a
, b
) CONCAT2
(a
, b
)
48 #define CONCAT2
(a
, b
) a ## b
50 /* Use the right prefix for
global labels.
*/
52 #define SYM
(x
) CONCAT1
(__USER_LABEL_PREFIX__
, x
)
54 /* Note that X is a function.
*/
57 #define FUNC
(x
) .
type SYM
(x
),function
59 /* The .
proc pseudo
-op is accepted
, but ignored
, by GAS. We could just
60 define
this to the empty string for non
-ELF systems
, but defining it
61 to .
proc means that the information is available to the assembler if
66 /* Use the right prefix for registers.
*/
68 #define REG
(x
) CONCAT1
(__REGISTER_PREFIX__
, x
)
70 /* Use the right prefix for immediate values.
*/
72 #define IMM
(x
) CONCAT1
(__IMMEDIATE_PREFIX__
, x
)
93 /* Provide a few macros to allow for PIC code support.
94 * With PIC
, data is stored A5 relative so we
've got to take a bit of special
95 * care to ensure that all loads of global data is via A5. PIC also requires
96 * jumps and subroutine calls to be PC relative rather than absolute. We cheat
97 * a little on this and in the PIC case, we use short offset branches and
98 * hope that the final object code is within range (which it should be).
102 /* Non PIC (absolute/relocatable) versions */
112 .macro PICLEA sym, reg
116 .macro PICPEA sym, areg
122 # if defined (__uClinux__)
124 /* Versions for uClinux */
126 # if defined(__ID_SHARED_LIBRARY__)
128 /* -mid-shared-library versions */
130 .macro PICLEA sym, reg
131 movel a5@(_current_shared_library_a5_offset_), \reg
132 movel \sym@GOT(\reg), \reg
135 .macro PICPEA sym, areg
136 movel a5@(_current_shared_library_a5_offset_), \areg
137 movel \sym@GOT(\areg), sp@-
150 # else /* !__ID_SHARED_LIBRARY__ */
152 /* Versions for -msep-data */
154 .macro PICLEA sym, reg
155 movel \sym@GOT(a5), \reg
158 .macro PICPEA sym, areg
159 movel \sym@GOT(a5), sp@-
163 #if defined (__mcoldfire__) && !defined (__mcfisab__) && !defined (__mcfisac__)
172 /* ISA C has no bra.l instruction, and since this assembly file
173 gets assembled into multiple object files, we avoid the
174 bra instruction entirely. */
175 #if defined (__mcoldfire__) && !defined (__mcfisab__)
185 # else /* !__uClinux__ */
187 /* Versions for Linux */
189 .macro PICLEA sym, reg
190 movel #_GLOBAL_OFFSET_TABLE_@GOTPC, \reg
191 lea (-6, pc, \reg), \reg
192 movel \sym@GOT(\reg), \reg
195 .macro PICPEA sym, areg
196 movel #_GLOBAL_OFFSET_TABLE_@GOTPC, \areg
197 lea (-6, pc, \areg), \areg
198 movel \sym@GOT(\areg), sp@-
202 #if defined (__mcoldfire__) && !defined (__mcfisab__) && !defined (__mcfisac__)
211 /* ISA C has no bra.l instruction, and since this assembly file
212 gets assembled into multiple object files, we avoid the
213 bra instruction entirely. */
214 #if defined (__mcoldfire__) && !defined (__mcfisab__)
228 | This is an attempt at a decent floating point (single, double and
229 | extended double) code for the GNU C compiler. It should be easy to
230 | adapt to other compilers (but beware of the local labels!).
232 | Starting date: 21 October, 1990
234 | It is convenient to introduce the notation (s,e,f) for a floating point
235 | number, where s=sign, e=exponent, f=fraction. We will call a floating
236 | point number fpn to abbreviate, independently of the precision.
237 | Let MAX_EXP be in each case the maximum exponent (255 for floats, 1023
238 | for doubles and 16383 for long doubles). We then have the following
240 | 1. Normalized fpns have 0 < e < MAX_EXP. They correspond to
241 | (-1)^s x 1.f x 2^(e-bias-1).
242 | 2. Denormalized fpns have e=0. They correspond to numbers of the form
243 | (-1)^s x 0.f x 2^(-bias).
244 | 3. +/-INFINITY have e=MAX_EXP, f=0.
245 | 4. Quiet NaN (Not a Number) have all bits set.
246 | 5. Signaling NaN (Not a Number) have s=0, e=MAX_EXP, f=1.
248 |=============================================================================
250 |=============================================================================
252 | This is the floating point condition code register (_fpCCR):
255 | short _exception_bits;
256 | short _trap_enable_bits;
257 | short _sticky_bits;
258 | short _rounding_mode;
260 | short _last_operation;
284 .word ROUND_TO_NEAREST
297 EBITS = __exception_bits - SYM (_fpCCR)
298 TRAPE = __trap_enable_bits - SYM (_fpCCR)
299 STICK = __sticky_bits - SYM (_fpCCR)
300 ROUND = __rounding_mode - SYM (_fpCCR)
301 FORMT = __format - SYM (_fpCCR)
302 LASTO = __last_operation - SYM (_fpCCR)
303 OPER1 = __operand1 - SYM (_fpCCR)
304 OPER2 = __operand2 - SYM (_fpCCR)
306 | The following exception types are supported:
307 INEXACT_RESULT = 0x0001
310 DIVIDE_BY_ZERO = 0x0008
311 INVALID_OPERATION = 0x0010
313 | The allowed rounding modes are:
315 ROUND_TO_NEAREST = 0 | round result to nearest representable value
316 ROUND_TO_ZERO = 1 | round result towards zero
317 ROUND_TO_PLUS = 2 | round result towards plus infinity
318 ROUND_TO_MINUS = 3 | round result towards minus infinity
320 | The allowed values of format are:
326 | The allowed values for the last operation are:
336 |=============================================================================
337 | __clear_sticky_bits
338 |=============================================================================
340 | The sticky bits are normally not cleared (thus the name), whereas the
341 | exception type and exception value reflect the last computation.
342 | This routine is provided to clear them (you can also write to _fpCCR,
343 | since it is globally visible).
345 .globl SYM (__clear_sticky_bit)
350 | void __clear_sticky_bits(void);
351 SYM (__clear_sticky_bit):
352 PICLEA SYM (_fpCCR),a0
353 #ifndef __mcoldfire__
354 movew IMM (0),a0@(STICK)
360 |=============================================================================
361 | $_exception_handler
362 |=============================================================================
364 .globl $_exception_handler
369 | This is the common exit point if an exception occurs.
370 | NOTE: it is NOT callable from C!
371 | It expects the exception type in d7, the format (SINGLE_FLOAT,
372 | DOUBLE_FLOAT or LONG_FLOAT) in d6, and the last operation code in d5.
373 | It sets the corresponding exception and sticky bits, and the format.
374 | Depending on the format if fills the corresponding slots for the
375 | operands which produced the exception (all this information is provided
376 | so if you write your own exception handlers you have enough information
377 | to deal with the problem).
378 | Then checks to see if the corresponding exception is trap-enabled,
379 | in which case it pushes the address of _fpCCR and traps through
380 | trap FPTRAP (15 for the moment).
385 PICLEA SYM (_fpCCR),a0
386 movew d7,a0@(EBITS) | set __exception_bits
387 #ifndef __mcoldfire__
388 orw d7,a0@(STICK) | and __sticky_bits
394 movew d6,a0@(FORMT) | and __format
395 movew d5,a0@(LASTO) | and __last_operation
397 | Now put the operands in place:
398 #ifndef __mcoldfire__
399 cmpw IMM (SINGLE_FLOAT),d6
401 cmpl IMM (SINGLE_FLOAT),d6
404 movel a6@(8),a0@(OPER1)
405 movel a6@(12),a0@(OPER1+4)
406 movel a6@(16),a0@(OPER2)
407 movel a6@(20),a0@(OPER2+4)
409 1: movel a6@(8),a0@(OPER1)
410 movel a6@(12),a0@(OPER2)
412 | And check whether the exception is trap-enabled:
413 #ifndef __mcoldfire__
414 andw a0@(TRAPE),d7 | is exception trap-enabled?
421 PICPEA SYM (_fpCCR),a1 | yes, push address of _fpCCR
422 trap IMM (FPTRAP) | and trap
423 #ifndef __mcoldfire__
424 1: moveml sp@+,d2-d7 | restore data registers
427 | XXX if frame pointer is ever removed, stack pointer must
432 #endif /* L_floatex */
437 .globl SYM (__mulsi3)
439 movew sp@(4), d0 /* x0 -> d0 */
440 muluw sp@(10), d0 /* x0*y1 */
441 movew sp@(6), d1 /* x1 -> d1 */
442 muluw sp@(8), d1 /* x1*y0 */
443 #ifndef __mcoldfire__
450 movew sp@(6), d1 /* x1 -> d1 */
451 muluw sp@(10), d1 /* x1*y1 */
455 #endif /* L_mulsi3 */
460 .globl SYM (__udivsi3)
462 #ifndef __mcoldfire__
464 movel sp@(12), d1 /* d1 = divisor */
465 movel sp@(8), d0 /* d0 = dividend */
467 cmpl IMM (0x10000), d1 /* divisor >= 2 ^ 16 ? */
468 jcc L3 /* then try next algorithm */
472 divu d1, d2 /* high quotient in lower word */
473 movew d2, d0 /* save high quotient */
475 movew sp@(10), d2 /* get low dividend + high rest */
476 divu d1, d2 /* low quotient */
480 L3: movel d1, d2 /* use d2 as divisor backup */
481 L4: lsrl IMM (1), d1 /* shift divisor */
482 lsrl IMM (1), d0 /* shift dividend */
483 cmpl IMM (0x10000), d1 /* still divisor >= 2 ^ 16 ? */
485 divu d1, d0 /* now we have 16-bit divisor */
486 andl IMM (0xffff), d0 /* mask out divisor, ignore remainder */
488 /* Multiply the 16-bit tentative quotient with the 32-bit divisor. Because of
489 the operand ranges, this might give a 33-bit product. If this product is
490 greater than the dividend, the tentative quotient was too large. */
492 mulu d0, d1 /* low part, 32 bits */
494 mulu d0, d2 /* high part, at most 17 bits */
495 swap d2 /* align high part with low part */
496 tstw d2 /* high part 17 bits? */
497 jne L5 /* if 17 bits, quotient was too large */
498 addl d2, d1 /* add parts */
499 jcs L5 /* if sum is 33 bits, quotient was too large */
500 cmpl sp@(8), d1 /* compare the sum with the dividend */
501 jls L6 /* if sum > dividend, quotient was too large */
502 L5: subql IMM (1), d0 /* adjust quotient */
507 #else /* __mcoldfire__ */
509 /* ColdFire implementation of non-restoring division algorithm from
510 Hennessy & Patterson, Appendix A. */
517 L1: addl d0,d0 | shift reg pair (p,a) one bit left
519 movl d2,d3 | subtract b from p, store in tmp.
521 jcs L2 | if no carry,
522 bset IMM (0),d0 | set the low order bit of a to 1,
523 movl d3,d2 | and store tmp in p.
526 moveml sp@,d2-d4 | restore data registers
529 #endif /* __mcoldfire__ */
531 #endif /* L_udivsi3 */
536 .globl SYM (__divsi3)
540 moveq IMM (1), d2 /* sign of result stored in d2 (=1 or =-1) */
541 movel sp@(12), d1 /* d1 = divisor */
544 #ifndef __mcoldfire__
545 negb d2 /* change sign because divisor <0 */
547 negl d2 /* change sign because divisor <0 */
549 L1: movel sp@(8), d0 /* d0 = dividend */
552 #ifndef __mcoldfire__
560 PICCALL SYM (__udivsi3) /* divide abs(dividend) by abs(divisor) */
569 #endif /* L_divsi3 */
574 .globl SYM (__umodsi3)
576 movel sp@(8), d1 /* d1 = divisor */
577 movel sp@(4), d0 /* d0 = dividend */
580 PICCALL SYM (__udivsi3)
582 movel sp@(8), d1 /* d1 = divisor */
583 #ifndef __mcoldfire__
586 PICCALL SYM (__mulsi3) /* d0 = (a/b)*b */
591 movel sp@(4), d1 /* d1 = dividend */
592 subl d0, d1 /* d1 = a - (a/b)*b */
595 #endif /* L_umodsi3 */
600 .globl SYM (__modsi3)
602 movel sp@(8), d1 /* d1 = divisor */
603 movel sp@(4), d0 /* d0 = dividend */
606 PICCALL SYM (__divsi3)
608 movel sp@(8), d1 /* d1 = divisor */
609 #ifndef __mcoldfire__
612 PICCALL SYM (__mulsi3) /* d0 = (a/b)*b */
617 movel sp@(4), d1 /* d1 = dividend */
618 subl d0, d1 /* d1 = a - (a/b)*b */
621 #endif /* L_modsi3 */
627 .globl $_exception_handler
629 QUIET_NaN = 0xffffffff
633 DBL_MAX_EXP = D_MAX_EXP - D_BIAS
634 DBL_MIN_EXP = 1 - D_BIAS
637 INEXACT_RESULT = 0x0001
640 DIVIDE_BY_ZERO = 0x0008
641 INVALID_OPERATION = 0x0010
655 ROUND_TO_NEAREST = 0 | round result to nearest representable value
656 ROUND_TO_ZERO = 1 | round result towards zero
657 ROUND_TO_PLUS = 2 | round result towards plus infinity
658 ROUND_TO_MINUS = 3 | round result towards minus infinity
662 .globl SYM (__adddf3)
663 .globl SYM (__subdf3)
664 .globl SYM (__muldf3)
665 .globl SYM (__divdf3)
666 .globl SYM (__negdf2)
667 .globl SYM (__cmpdf2)
668 .globl SYM (__cmpdf2_internal)
669 .hidden SYM (__cmpdf2_internal)
674 | These are common routines to return and signal exceptions.
677 | Return and signal a denormalized number
679 movew IMM (INEXACT_RESULT+UNDERFLOW),d7
680 moveq IMM (DOUBLE_FLOAT),d6
681 PICJUMP $_exception_handler
685 | Return a properly signed INFINITY and set the exception flags
686 movel IMM (0x7ff00000),d0
689 movew IMM (INEXACT_RESULT+OVERFLOW),d7
690 moveq IMM (DOUBLE_FLOAT),d6
691 PICJUMP $_exception_handler
694 | Return 0 and set the exception flags
697 movew IMM (INEXACT_RESULT+UNDERFLOW),d7
698 moveq IMM (DOUBLE_FLOAT),d6
699 PICJUMP $_exception_handler
702 | Return a quiet NaN and set the exception flags
703 movel IMM (QUIET_NaN),d0
705 movew IMM (INEXACT_RESULT+INVALID_OPERATION),d7
706 moveq IMM (DOUBLE_FLOAT),d6
707 PICJUMP $_exception_handler
710 | Return a properly signed INFINITY and set the exception flags
711 movel IMM (0x7ff00000),d0
714 movew IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
715 moveq IMM (DOUBLE_FLOAT),d6
716 PICJUMP $_exception_handler
718 |=============================================================================
719 |=============================================================================
720 | double precision routines
721 |=============================================================================
722 |=============================================================================
724 | A double precision floating point number (double) has the format:
727 | unsigned int sign : 1; /* sign bit */
728 | unsigned int exponent : 11; /* exponent, shifted by 126 */
729 | unsigned int fraction : 52; /* fraction */
732 | Thus sizeof(double) = 8 (64 bits).
734 | All the routines are callable from C programs, and return the result
735 | in the register pair d0-d1. They also preserve all registers except
738 |=============================================================================
740 |=============================================================================
742 | double __subdf3(double, double);
745 bchg IMM (31),sp@(12) | change sign of second operand
746 | and fall through, so we always add
747 |=============================================================================
749 |=============================================================================
751 | double __adddf3(double, double);
754 #ifndef __mcoldfire__
755 link a6,IMM (0) | everything will be done in registers
756 moveml d2-d7,sp@- | save all data registers and a2 (but d0-d1)
761 movel a6@(8),d0 | get first operand
763 movel a6@(16),d2 | get second operand
766 movel d0,d7 | get d0's sign bit
in d7
'
767 addl d1,d1 | check and clear sign bit of a, and gain one
768 addxl d0,d0 | bit of extra precision
769 beq Ladddf$b | if zero return second operand
771 movel d2,d6 | save sign in d6
772 addl d3,d3 | get rid of sign bit and gain one bit of
773 addxl d2,d2 | extra precision
774 beq Ladddf$a | if zero return first operand
776 andl IMM (0x80000000),d7 | isolate a's sign bit
'
777 swap d6 | and also b's sign bit
'
778 #ifndef __mcoldfire__
779 andw IMM (0x8000),d6 |
780 orw d6,d7 | and combine them into d7, so that a's sign
'
781 | bit is in the high word and b's is
in the
'
782 | low word, so d6 is free to be used
787 movel d7,a0 | now save d7 into a0, so d7 is free to
790 | Get the exponents and check for denormalized and/or infinity.
792 movel IMM (0x001fffff),d6 | mask for the fraction
793 movel IMM (0x00200000),d7 | mask to put hidden bit back
796 andl d6,d0 | get fraction in d0
797 notl d6 | make d6 into mask for the exponent
798 andl d6,d4 | get exponent in d4
799 beq Ladddf$a$den | branch if a is denormalized
800 cmpl d6,d4 | check for INFINITY or NaN
802 orl d7,d0 | and put hidden bit back
804 swap d4 | shift right exponent so that it starts
805 #ifndef __mcoldfire__
806 lsrw IMM (5),d4 | in bit 0 and not bit 20
808 lsrl IMM (5),d4 | in bit 0 and not bit 20
810 | Now we have a's exponent
in d4
and fraction
in d0
-d1
'
811 movel d2,d5 | save b to get exponent
812 andl d6,d5 | get exponent in d5
813 beq Ladddf$b$den | branch if b is denormalized
814 cmpl d6,d5 | check for INFINITY or NaN
816 notl d6 | make d6 into mask for the fraction again
817 andl d6,d2 | and get fraction in d2
818 orl d7,d2 | and put hidden bit back
820 swap d5 | shift right exponent so that it starts
821 #ifndef __mcoldfire__
822 lsrw IMM (5),d5 | in bit 0 and not bit 20
824 lsrl IMM (5),d5 | in bit 0 and not bit 20
827 | Now we have b's exponent
in d5
and fraction
in d2
-d3.
'
829 | The situation now is as follows: the signs are combined in a0, the
830 | numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
831 | and d5 (b). To do the rounding correctly we need to keep all the
832 | bits until the end, so we need to use d0-d1-d2-d3 for the first number
833 | and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
834 | exponents in a2-a3.
836 #ifndef __mcoldfire__
837 moveml a2-a3,sp@- | save the address registers
844 movel d4,a2 | save the exponents
847 movel IMM (0),d7 | and move the numbers around
854 | Here we shift the numbers until the exponents are the same, and put
855 | the largest exponent in a2.
856 #ifndef __mcoldfire__
857 exg d4,a2 | get exponents back
859 cmpw d4,d5 | compare the exponents
861 movel d4,a4 | get exponents back
867 cmpl d4,d5 | compare the exponents
869 beq Ladddf$3 | if equal don't shift
'
870 bhi 9f | branch if second exponent is higher
872 | Here we have a's exponent larger than b
's, so we have to shift b. We do
873 | this by using as counter d2:
874 1: movew d4,d2 | move largest exponent to d2
875 #ifndef __mcoldfire__
876 subw d5,d2 | and subtract second exponent
877 exg d4,a2 | get back the longs we saved
880 subl d5,d2 | and subtract second exponent
881 movel d4,a4 | get back the longs we saved
888 | if difference is too large we don't shift
(actually
, we can just exit
) '
889 #ifndef __mcoldfire__
890 cmpw IMM (DBL_MANT_DIG+2),d2
892 cmpl IMM (DBL_MANT_DIG+2),d2
895 #ifndef __mcoldfire__
896 cmpw IMM (32),d2 | if difference >= 32, shift by longs
898 cmpl IMM (32),d2 | if difference >= 32, shift by longs
902 #ifndef __mcoldfire__
903 cmpw IMM (16),d2 | if difference >= 16, shift by words
905 cmpl IMM (16),d2 | if difference >= 16, shift by words
908 bra 3f | enter dbra loop
911 #ifndef __mcoldfire__
932 #ifndef __mcoldfire__
946 #ifndef __mcoldfire__
961 #ifndef __mcoldfire__
969 #ifndef __mcoldfire__
972 subw d5,d6 | keep d5 (largest exponent) in d4
987 | if difference is too large we don't shift
(actually
, we can just exit
) '
988 #ifndef __mcoldfire__
989 cmpw IMM (DBL_MANT_DIG+2),d6
991 cmpl IMM (DBL_MANT_DIG+2),d6
994 #ifndef __mcoldfire__
995 cmpw IMM (32),d6 | if difference >= 32, shift by longs
997 cmpl IMM (32),d6 | if difference >= 32, shift by longs
1001 #ifndef __mcoldfire__
1002 cmpw IMM (16),d6 | if difference >= 16, shift by words
1004 cmpl IMM (16),d6 | if difference >= 16, shift by words
1007 bra 3f | enter dbra loop
1010 #ifndef __mcoldfire__
1031 #ifndef __mcoldfire__
1045 #ifndef __mcoldfire__
1060 #ifndef __mcoldfire__
1067 #ifndef __mcoldfire__
1079 | Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
1082 | Here we have to decide whether to add or subtract the numbers:
1083 #ifndef __mcoldfire__
1084 exg d7,a0 | get the signs
1085 exg d6,a3 | a3 is free to be used
1095 movew IMM (0),d7 | get a's sign
in d7
'
1097 movew IMM (0),d6 | and b's sign
in d6
'
1098 eorl d7,d6 | compare the signs
1099 bmi Lsubdf$0 | if the signs are different we have
1101 #ifndef __mcoldfire__
1102 exg d7,a0 | else we add the numbers
1117 movel a2,d4 | return exponent to d4
1119 andl IMM (0x80000000),d7 | d7 now has the sign
1121 #ifndef __mcoldfire__
1129 | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
1130 | the case of denormalized numbers in the rounding routine itself).
1131 | As in the addition (not in the subtraction!) we could have set
1132 | one more bit we check this:
1133 btst IMM (DBL_MANT_DIG+1),d0
1135 #ifndef __mcoldfire__
1158 lea pc@(Ladddf$5),a0 | to return from rounding routine
1159 PICLEA SYM (_fpCCR),a1 | check the rounding mode
1160 #ifdef __mcoldfire__
1163 movew a1@(6),d6 | rounding mode in d6
1164 beq Lround$to$nearest
1165 #ifndef __mcoldfire__
1166 cmpw IMM (ROUND_TO_PLUS),d6
1168 cmpl IMM (ROUND_TO_PLUS),d6
1174 | Put back the exponent and check for overflow
1175 #ifndef __mcoldfire__
1176 cmpw IMM (0x7ff),d4 | is the exponent big?
1178 cmpl IMM (0x7ff),d4 | is the exponent big?
1181 bclr IMM (DBL_MANT_DIG-1),d0
1182 #ifndef __mcoldfire__
1183 lslw IMM (4),d4 | put exponent back into position
1185 lsll IMM (4),d4 | put exponent back into position
1188 #ifndef __mcoldfire__
1200 | Here we do the subtraction.
1201 #ifndef __mcoldfire__
1202 exg d7,a0 | put sign back in a0
1216 beq Ladddf$ret$1 | if zero just exit
1217 bpl 1f | if positive skip the following
1219 bchg IMM (31),d7 | change sign bit in d7
1223 negxl d1 | and negate result
1226 movel a2,d4 | return exponent to d4
1228 andl IMM (0x80000000),d7 | isolate sign bit
1229 #ifndef __mcoldfire__
1237 | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
1238 | the case of denormalized numbers in the rounding routine itself).
1239 | As in the addition (not in the subtraction!) we could have set
1240 | one more bit we check this:
1241 btst IMM (DBL_MANT_DIG+1),d0
1243 #ifndef __mcoldfire__
1266 lea pc@(Lsubdf$1),a0 | to return from rounding routine
1267 PICLEA SYM (_fpCCR),a1 | check the rounding mode
1268 #ifdef __mcoldfire__
1271 movew a1@(6),d6 | rounding mode in d6
1272 beq Lround$to$nearest
1273 #ifndef __mcoldfire__
1274 cmpw IMM (ROUND_TO_PLUS),d6
1276 cmpl IMM (ROUND_TO_PLUS),d6
1282 | Put back the exponent and sign (we don't have overflow
).
'
1283 bclr IMM (DBL_MANT_DIG-1),d0
1284 #ifndef __mcoldfire__
1285 lslw IMM (4),d4 | put exponent back into position
1287 lsll IMM (4),d4 | put exponent back into position
1290 #ifndef __mcoldfire__
1298 | If one of the numbers was too small (difference of exponents >=
1299 | DBL_MANT_DIG+1) we return the other (and now we don't have to
'
1300 | check for finiteness or zero).
1302 #ifndef __mcoldfire__
1311 PICLEA SYM (_fpCCR),a0
1313 #ifndef __mcoldfire__
1314 moveml sp@+,d2-d7 | restore data registers
1317 | XXX if frame pointer is ever removed, stack pointer must
1320 unlk a6 | and return
1324 #ifndef __mcoldfire__
1333 PICLEA SYM (_fpCCR),a0
1335 #ifndef __mcoldfire__
1336 moveml sp@+,d2-d7 | restore data registers
1339 | XXX if frame pointer is ever removed, stack pointer must
1342 unlk a6 | and return
1346 movel d7,d4 | d7 contains 0x00200000
1350 movel d7,d5 | d7 contains 0x00200000
1355 | Return b (if a is zero)
1358 bne 1f | Check if b is -0
1359 cmpl IMM (0x80000000),d0
1361 andl IMM (0x80000000),d7 | Use the sign of a
1369 | Check for NaN and +/-INFINITY.
1371 andl IMM (0x80000000),d7 |
1373 cmpl IMM (0x7ff00000),d0 |
1375 movel d0,d0 | check for zero, since we don't
'
1376 bne Ladddf$ret | want to return -0 by mistake
1380 andl IMM (0x000fffff),d0 | check for NaN (nonzero fraction)
1386 #ifndef __mcoldfire__
1387 moveml sp@+,a2-a3 | restore regs and exit
1396 PICLEA SYM (_fpCCR),a0
1398 orl d7,d0 | put sign bit back
1399 #ifndef __mcoldfire__
1403 | XXX if frame pointer is ever removed, stack pointer must
1410 | Return a denormalized number.
1411 #ifndef __mcoldfire__
1412 lsrl IMM (1),d0 | shift right once more
1425 | This could be faster but it is not worth the effort, since it is not
1426 | executed very often. We sacrifice speed for clarity here.
1427 movel a6@(8),d0 | get the numbers back (remember that we
1428 movel a6@(12),d1 | did some processing already)
1431 movel IMM (0x7ff00000),d4 | useful constant (INFINITY)
1432 movel d0,d7 | save sign bits
1434 bclr IMM (31),d0 | clear sign bits
1436 | We know that one of them is either NaN of +/-INFINITY
1437 | Check for NaN (if either one is NaN return NaN)
1438 cmpl d4,d0 | check first a (d0)
1439 bhi Ld$inop | if d0 > 0x7ff00000 or equal and
1441 tstl d1 | d1 > 0, a is NaN
1443 2: cmpl d4,d2 | check now b (d1)
1449 | Now comes the check for +/-INFINITY. We know that both are (maybe not
1450 | finite) numbers, but we have to check if both are infinite whether we
1451 | are adding or subtracting them.
1452 eorl d7,d6 | to check sign bits
1454 andl IMM (0x80000000),d7 | get (common) sign bit
1457 | We know one (or both) are infinite, so we test for equality between the
1458 | two numbers (if they are equal they have to be infinite both, so we
1460 cmpl d2,d0 | are both infinite?
1461 bne 1f | if d0 <> d2 they are not equal
1462 cmpl d3,d1 | if d0 == d2 test d3 and d1
1463 beq Ld$inop | if equal return NaN
1465 andl IMM (0x80000000),d7 | get a's sign bit
'
1466 cmpl d4,d0 | test now for infinity
1467 beq Ld$infty | if a is INFINITY return with this sign
1468 bchg IMM (31),d7 | else we know b is INFINITY and has
1469 bra Ld$infty | the opposite sign
1471 |=============================================================================
1473 |=============================================================================
1475 | double __muldf3(double, double);
1478 #ifndef __mcoldfire__
1485 movel a6@(8),d0 | get a into d0-d1
1487 movel a6@(16),d2 | and b into d2-d3
1489 movel d0,d7 | d7 will hold the sign of the product
1491 andl IMM (0x80000000),d7 |
1492 movel d7,a0 | save sign bit into a0
1493 movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
1494 movel d7,d6 | another (mask for fraction)
1496 bclr IMM (31),d0 | get rid of a's sign bit
'
1499 beq Lmuldf$a$0 | branch if a is zero
1501 bclr IMM (31),d2 | get rid of b's sign bit
'
1504 beq Lmuldf$b$0 | branch if b is zero
1506 cmpl d7,d0 | is a big?
1507 bhi Lmuldf$inop | if a is NaN return NaN
1508 beq Lmuldf$a$nf | we still have to check d1 and b ...
1509 cmpl d7,d2 | now compare b with INFINITY
1510 bhi Lmuldf$inop | is b NaN?
1511 beq Lmuldf$b$nf | we still have to check d3 ...
1512 | Here we have both numbers finite and nonzero (and with no sign bit).
1513 | Now we get the exponents into d4 and d5.
1514 andl d7,d4 | isolate exponent in d4
1515 beq Lmuldf$a$den | if exponent zero, have denormalized
1516 andl d6,d0 | isolate fraction
1517 orl IMM (0x00100000),d0 | and put hidden bit back
1518 swap d4 | I like exponents in the first byte
1519 #ifndef __mcoldfire__
1528 orl IMM (0x00100000),d2 | and put hidden bit back
1530 #ifndef __mcoldfire__
1536 #ifndef __mcoldfire__
1537 addw d5,d4 | add exponents
1538 subw IMM (D_BIAS+1),d4 | and subtract bias (plus one)
1540 addl d5,d4 | add exponents
1541 subl IMM (D_BIAS+1),d4 | and subtract bias (plus one)
1544 | We are now ready to do the multiplication. The situation is as follows:
1545 | both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were
1546 | denormalized to start with!), which means that in the product bit 104
1547 | (which will correspond to bit 8 of the fourth long) is set.
1549 | Here we have to do the product.
1550 | To do it we have to juggle the registers back and forth, as there are not
1551 | enough to keep everything in them. So we use the address registers to keep
1552 | some intermediate data.
1554 #ifndef __mcoldfire__
1555 moveml a2-a3,sp@- | save a2 and a3 for temporary use
1561 movel IMM (0),a2 | a2 is a null register
1562 movel d4,a3 | and a3 will preserve the exponent
1564 | First, shift d2-d3 so bit 20 becomes bit 31:
1565 #ifndef __mcoldfire__
1566 rorl IMM (5),d2 | rotate d2 5 places right
1567 swap d2 | and swap it
1568 rorl IMM (5),d3 | do the same thing with d3
1570 movew d3,d6 | get the rightmost 11 bits of d3
1571 andw IMM (0x07ff),d6 |
1572 orw d6,d2 | and put them into d2
1573 andw IMM (0xf800),d3 | clear those bits in d3
1575 moveq IMM (11),d7 | left shift d2 11 bits
1577 movel d3,d6 | get a copy of d3
1578 lsll d7,d3 | left shift d3 11 bits
1579 andl IMM (0xffe00000),d6 | get the top 11 bits of d3
1580 moveq IMM (21),d7 | right shift them 21 bits
1582 orl d6,d2 | stick them at the end of d2
1585 movel d2,d6 | move b into d6-d7
1586 movel d3,d7 | move a into d4-d5
1587 movel d0,d4 | and clear d0-d1-d2-d3 (to put result)
1594 | We use a1 as counter:
1595 movel IMM (DBL_MANT_DIG-1),a1
1596 #ifndef __mcoldfire__
1605 #ifndef __mcoldfire__
1606 exg d7,a1 | put counter back in a1
1612 addl d3,d3 | shift sum once left
1618 bcc 2f | if bit clear skip the following
1619 #ifndef __mcoldfire__
1626 addl d5,d3 | else add a to the sum
1630 #ifndef __mcoldfire__
1638 #ifndef __mcoldfire__
1639 exg d7,a1 | put counter in d7
1640 dbf d7,1b | decrement and branch
1649 movel a3,d4 | restore exponent
1650 #ifndef __mcoldfire__
1658 | Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The
1659 | first thing to do now is to normalize it so bit 8 becomes bit
1660 | DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
1669 #ifndef __mcoldfire__
1699 | Now round, check for over- and underflow, and exit.
1700 movel a0,d7 | get sign bit back into d7
1701 moveq IMM (MULTIPLY),d5
1703 btst IMM (DBL_MANT_DIG+1-32),d0
1705 #ifndef __mcoldfire__
1720 moveq IMM (MULTIPLY),d5
1724 moveq IMM (MULTIPLY),d5
1725 movel a0,d7 | get sign bit back into d7
1726 tstl d3 | we know d2 == 0x7ff00000, so check d3
1727 bne Ld$inop | if d3 <> 0 b is NaN
1728 bra Ld$overflow | else we have overflow (since a is finite)
1731 moveq IMM (MULTIPLY),d5
1732 movel a0,d7 | get sign bit back into d7
1733 tstl d1 | we know d0 == 0x7ff00000, so check d1
1734 bne Ld$inop | if d1 <> 0 a is NaN
1735 bra Ld$overflow | else signal overflow
1737 | If either number is zero return zero, unless the other is +/-INFINITY or
1738 | NaN, in which case we return NaN.
1740 moveq IMM (MULTIPLY),d5
1741 #ifndef __mcoldfire__
1742 exg d2,d0 | put b (==0) into d0-d1
1743 exg d3,d1 | and a (with sign bit cleared) into d2-d3
1744 movel a0,d0 | set result sign
1746 movel d0,d2 | put a into d2-d3
1748 movel a0,d0 | put result zero into d0-d1
1753 movel a0,d0 | set result sign
1754 movel a6@(16),d2 | put b into d2-d3 again
1756 bclr IMM (31),d2 | clear sign bit
1757 1: cmpl IMM (0x7ff00000),d2 | check for non-finiteness
1758 bge Ld$inop | in case NaN or +/-INFINITY return NaN
1759 PICLEA SYM (_fpCCR),a0
1761 #ifndef __mcoldfire__
1765 | XXX if frame pointer is ever removed, stack pointer must
1771 | If a number is denormalized we put an exponent of 1 but do not put the
1772 | hidden bit back into the fraction; instead we shift left until bit 21
1773 | (the hidden bit) is set, adjusting the exponent accordingly. We do this
1774 | to ensure that the product of the fractions is close to 1.
1778 1: addl d1,d1 | shift a left until bit 20 is set
1780 #ifndef __mcoldfire__
1781 subw IMM (1),d4 | and adjust exponent
1783 subl IMM (1),d4 | and adjust exponent
1792 1: addl d3,d3 | shift b left until bit 20 is set
1794 #ifndef __mcoldfire__
1795 subw IMM (1),d5 | and adjust exponent
1797 subql IMM (1),d5 | and adjust exponent
1804 |=============================================================================
1806 |=============================================================================
1808 | double __divdf3(double, double);
1811 #ifndef __mcoldfire__
1818 movel a6@(8),d0 | get a into d0-d1
1820 movel a6@(16),d2 | and b into d2-d3
1822 movel d0,d7 | d7 will hold the sign of the result
1824 andl IMM (0x80000000),d7
1825 movel d7,a0 | save sign into a0
1826 movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
1827 movel d7,d6 | another (mask for fraction)
1829 bclr IMM (31),d0 | get rid of a's sign bit
'
1832 beq Ldivdf$a$0 | branch if a is zero
1834 bclr IMM (31),d2 | get rid of b's sign bit
'
1837 beq Ldivdf$b$0 | branch if b is zero
1839 cmpl d7,d0 | is a big?
1840 bhi Ldivdf$inop | if a is NaN return NaN
1841 beq Ldivdf$a$nf | if d0 == 0x7ff00000 we check d1
1842 cmpl d7,d2 | now compare b with INFINITY
1843 bhi Ldivdf$inop | if b is NaN return NaN
1844 beq Ldivdf$b$nf | if d2 == 0x7ff00000 we check d3
1845 | Here we have both numbers finite and nonzero (and with no sign bit).
1846 | Now we get the exponents into d4 and d5 and normalize the numbers to
1847 | ensure that the ratio of the fractions is around 1. We do this by
1848 | making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
1849 | set, even if they were denormalized to start with.
1850 | Thus, the result will satisfy: 2 > result > 1/2.
1851 andl d7,d4 | and isolate exponent in d4
1852 beq Ldivdf$a$den | if exponent is zero we have a denormalized
1853 andl d6,d0 | and isolate fraction
1854 orl IMM (0x00100000),d0 | and put hidden bit back
1855 swap d4 | I like exponents in the first byte
1856 #ifndef __mcoldfire__
1865 orl IMM (0x00100000),d2
1867 #ifndef __mcoldfire__
1873 #ifndef __mcoldfire__
1874 subw d5,d4 | subtract exponents
1875 addw IMM (D_BIAS),d4 | and add bias
1877 subl d5,d4 | subtract exponents
1878 addl IMM (D_BIAS),d4 | and add bias
1881 | We are now ready to do the division. We have prepared things in such a way
1882 | that the ratio of the fractions will be less than 2 but greater than 1/2.
1883 | At this point the registers in use are:
1884 | d0-d1 hold a (first operand, bit DBL_MANT_DIG-32=0, bit
1885 | DBL_MANT_DIG-1-32=1)
1886 | d2-d3 hold b (second operand, bit DBL_MANT_DIG-32=1)
1887 | d4 holds the difference of the exponents, corrected by the bias
1888 | a0 holds the sign of the ratio
1890 | To do the rounding correctly we need to keep information about the
1891 | nonsignificant bits. One way to do this would be to do the division
1892 | using four registers; another is to use two registers (as originally
1893 | I did), but use a sticky bit to preserve information about the
1894 | fractional part. Note that we can keep that info in a1, which is not
1896 movel IMM (0),d6 | d6-d7 will hold the result
1898 movel IMM (0),a1 | and a1 will hold the sticky bit
1900 movel IMM (DBL_MANT_DIG-32+1),d5
1902 1: cmpl d0,d2 | is a < b?
1903 bhi 3f | if b > a skip the following
1904 beq 4f | if d0==d2 check d1 and d3
1906 subxl d2,d0 | a <-- a - b
1907 bset d5,d6 | set the corresponding bit in d6
1908 3: addl d1,d1 | shift a by 1
1910 #ifndef __mcoldfire__
1911 dbra d5,1b | and branch back
1917 4: cmpl d1,d3 | here d0==d2, so check d1 and d3
1918 bhi 3b | if d1 > d2 skip the subtraction
1919 bra 2b | else go do it
1921 | Here we have to start setting the bits in the second long.
1922 movel IMM (31),d5 | again d5 is counter
1924 1: cmpl d0,d2 | is a < b?
1925 bhi 3f | if b > a skip the following
1926 beq 4f | if d0==d2 check d1 and d3
1928 subxl d2,d0 | a <-- a - b
1929 bset d5,d7 | set the corresponding bit in d7
1930 3: addl d1,d1 | shift a by 1
1932 #ifndef __mcoldfire__
1933 dbra d5,1b | and branch back
1939 4: cmpl d1,d3 | here d0==d2, so check d1 and d3
1940 bhi 3b | if d1 > d2 skip the subtraction
1941 bra 2b | else go do it
1943 | Now go ahead checking until we hit a one, which we store in d2.
1944 movel IMM (DBL_MANT_DIG),d5
1945 1: cmpl d2,d0 | is a < b?
1946 bhi 4f | if b < a, exit
1947 beq 3f | if d0==d2 check d1 and d3
1948 2: addl d1,d1 | shift a by 1
1950 #ifndef __mcoldfire__
1951 dbra d5,1b | and branch back
1956 movel IMM (0),d2 | here no sticky bit was found
1959 3: cmpl d1,d3 | here d0==d2, so check d1 and d3
1960 bhi 2b | if d1 > d2 go back
1962 | Here put the sticky bit in d2-d3 (in the position which actually corresponds
1963 | to it; if you don't do
this the algorithm loses
in some cases
).
'
1966 #ifndef __mcoldfire__
1967 subw IMM (DBL_MANT_DIG),d5
1971 subl IMM (DBL_MANT_DIG),d5
1978 #ifndef __mcoldfire__
1985 | Finally we are finished! Move the longs in the address registers to
1986 | their final destination:
1991 | Here we have finished the division, with the result in d0-d1-d2-d3, with
1992 | 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
1993 | If it is not, then definitely bit 21 is set. Normalize so bit 22 is
1995 btst IMM (DBL_MANT_DIG-32+1),d0
1997 #ifndef __mcoldfire__
2020 | Now round, check for over- and underflow, and exit.
2021 movel a0,d7 | restore sign bit to d7
2022 moveq IMM (DIVIDE),d5
2026 moveq IMM (DIVIDE),d5
2030 | If a is zero check to see whether b is zero also. In that case return
2031 | NaN; then check if b is NaN, and return NaN also in that case. Else
2032 | return a properly signed zero.
2033 moveq IMM (DIVIDE),d5
2037 beq Ld$inop | if b is also zero return NaN
2038 cmpl IMM (0x7ff00000),d2 | check for NaN
2043 1: movel a0,d0 | else return signed zero
2045 PICLEA SYM (_fpCCR),a0 | clear exception flags
2047 #ifndef __mcoldfire__
2051 | XXX if frame pointer is ever removed, stack pointer must
2058 moveq IMM (DIVIDE),d5
2059 | If we got here a is not zero. Check if a is NaN; in that case return NaN,
2060 | else return +/-INFINITY. Remember that a is in d0 with the sign bit
2062 movel a0,d7 | put a's sign bit back
in d7
'
2063 cmpl IMM (0x7ff00000),d0 | compare d0 with INFINITY
2064 bhi Ld$inop | if larger it is NaN
2067 bra Ld$div$0 | else signal DIVIDE_BY_ZERO
2070 moveq IMM (DIVIDE),d5
2071 | If d2 == 0x7ff00000 we have to check d3.
2073 bne Ld$inop | if d3 <> 0, b is NaN
2074 bra Ld$underflow | else b is +/-INFINITY, so signal underflow
2077 moveq IMM (DIVIDE),d5
2078 | If d0 == 0x7ff00000 we have to check d1.
2080 bne Ld$inop | if d1 <> 0, a is NaN
2081 | If a is INFINITY we have to check b
2082 cmpl d7,d2 | compare b with INFINITY
2083 bge Ld$inop | if b is NaN or INFINITY return NaN
2086 bra Ld$overflow | else return overflow
2088 | If a number is denormalized we put an exponent of 1 but do not put the
2089 | bit back into the fraction.
2093 1: addl d1,d1 | shift a left until bit 20 is set
2095 #ifndef __mcoldfire__
2096 subw IMM (1),d4 | and adjust exponent
2098 subl IMM (1),d4 | and adjust exponent
2100 btst IMM (DBL_MANT_DIG-32-1),d0
2107 1: addl d3,d3 | shift b left until bit 20 is set
2109 #ifndef __mcoldfire__
2110 subw IMM (1),d5 | and adjust exponent
2112 subql IMM (1),d5 | and adjust exponent
2114 btst IMM (DBL_MANT_DIG-32-1),d2
2119 | This is a common exit point for __muldf3 and __divdf3. When they enter
2120 | this point the sign of the result is in d7, the result in d0-d1, normalized
2121 | so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.
2123 | First check for underlow in the exponent:
2124 #ifndef __mcoldfire__
2125 cmpw IMM (-DBL_MANT_DIG-1),d4
2127 cmpl IMM (-DBL_MANT_DIG-1),d4
2130 | It could happen that the exponent is less than 1, in which case the
2131 | number is denormalized. In this case we shift right and adjust the
2132 | exponent until it becomes 1 or the fraction is zero (in the latter case
2133 | we signal underflow and return zero).
2135 movel IMM (0),d6 | use d6-d7 to collect bits flushed right
2136 movel d6,d7 | use d6-d7 to collect bits flushed right
2137 #ifndef __mcoldfire__
2138 cmpw IMM (1),d4 | if the exponent is less than 1 we
2140 cmpl IMM (1),d4 | if the exponent is less than 1 we
2142 bge 2f | have to shift right (denormalize)
2144 #ifndef __mcoldfire__
2145 addw IMM (1),d4 | adjust the exponent
2146 lsrl IMM (1),d0 | shift right once
2152 cmpw IMM (1),d4 | is the exponent 1 already?
2154 addl IMM (1),d4 | adjust the exponent
2176 cmpl IMM (1),d4 | is the exponent 1 already?
2178 beq 2f | if not loop back
2180 bra Ld$underflow | safety check, shouldn't execute
'
2181 2: orl d6,d2 | this is a trick so we don't lose
'
2182 orl d7,d3 | the bits which were flushed right
2183 movel a0,d7 | get back sign bit into d7
2184 | Now call the rounding routine (which takes care of denormalized numbers):
2185 lea pc@(Lround$0),a0 | to return from rounding routine
2186 PICLEA SYM (_fpCCR),a1 | check the rounding mode
2187 #ifdef __mcoldfire__
2190 movew a1@(6),d6 | rounding mode in d6
2191 beq Lround$to$nearest
2192 #ifndef __mcoldfire__
2193 cmpw IMM (ROUND_TO_PLUS),d6
2195 cmpl IMM (ROUND_TO_PLUS),d6
2201 | Here we have a correctly rounded result (either normalized or denormalized).
2203 | Here we should have either a normalized number or a denormalized one, and
2204 | the exponent is necessarily larger or equal to 1 (so we don't have to
'
2205 | check again for underflow!). We have to check for overflow or for a
2206 | denormalized number (which also signals underflow).
2207 | Check for overflow (i.e., exponent >= 0x7ff).
2208 #ifndef __mcoldfire__
2209 cmpw IMM (0x07ff),d4
2211 cmpl IMM (0x07ff),d4
2214 | Now check for a denormalized number (exponent==0):
2218 | Put back the exponents and sign and return.
2219 #ifndef __mcoldfire__
2220 lslw IMM (4),d4 | exponent back to fourth byte
2222 lsll IMM (4),d4 | exponent back to fourth byte
2224 bclr IMM (DBL_MANT_DIG-32-1),d0
2225 swap d0 | and put back exponent
2226 #ifndef __mcoldfire__
2232 orl d7,d0 | and sign also
2234 PICLEA SYM (_fpCCR),a0
2236 #ifndef __mcoldfire__
2240 | XXX if frame pointer is ever removed, stack pointer must
2246 |=============================================================================
2248 |=============================================================================
2250 | double __negdf2(double, double);
2253 #ifndef __mcoldfire__
2260 moveq IMM (NEGATE),d5
2261 movel a6@(8),d0 | get number to negate in d0-d1
2263 bchg IMM (31),d0 | negate
2264 movel d0,d2 | make a positive copy (for the tests)
2266 movel d2,d4 | check for zero
2268 beq 2f | if zero (either sign) return +zero
2269 cmpl IMM (0x7ff00000),d2 | compare to +INFINITY
2270 blt 1f | if finite, return
2271 bhi Ld$inop | if larger (fraction not zero) is NaN
2272 tstl d1 | if d2 == 0x7ff00000 check d1
2274 movel d0,d7 | else get sign and return INFINITY
2275 andl IMM (0x80000000),d7
2277 1: PICLEA SYM (_fpCCR),a0
2279 #ifndef __mcoldfire__
2283 | XXX if frame pointer is ever removed, stack pointer must
2291 |=============================================================================
2293 |=============================================================================
2299 | int __cmpdf2_internal(double, double, int);
2300 SYM (__cmpdf2_internal):
2301 #ifndef __mcoldfire__
2303 moveml d2-d7,sp@- | save registers
2308 moveq IMM (COMPARE),d5
2309 movel a6@(8),d0 | get first operand
2311 movel a6@(16),d2 | get second operand
2313 | First check if a and/or b are (+/-) zero and in that case clear
2315 movel d0,d6 | copy signs into d6 (a) and d7(b)
2316 bclr IMM (31),d0 | and clear signs in d0 and d2
2319 cmpl IMM (0x7ff00000),d0 | check for a == NaN
2320 bhi Lcmpd$inop | if d0 > 0x7ff00000, a is NaN
2321 beq Lcmpdf$a$nf | if equal can be INFINITY, so check d1
2322 movel d0,d4 | copy into d4 to test for zero
2326 cmpl IMM (0x7ff00000),d2 | check for b == NaN
2327 bhi Lcmpd$inop | if d2 > 0x7ff00000, b is NaN
2328 beq Lcmpdf$b$nf | if equal can be INFINITY, so check d3
2336 | If the signs are not equal check if a >= 0
2338 bpl Lcmpdf$a$gt$b | if (a >= 0 && b < 0) => a > b
2339 bmi Lcmpdf$b$gt$a | if (a < 0 && b >= 0) => a < b
2341 | If the signs are equal check for < 0
2344 | If both are negative exchange them
2345 #ifndef __mcoldfire__
2357 | Now that they are positive we just compare them as longs (does this also
2358 | work for denormalized numbers?).
2360 bhi Lcmpdf$b$gt$a | |b| > |a|
2361 bne Lcmpdf$a$gt$b | |b| < |a|
2362 | If we got here d0 == d2, so we compare d1 and d3.
2364 bhi Lcmpdf$b$gt$a | |b| > |a|
2365 bne Lcmpdf$a$gt$b | |b| < |a|
2366 | If we got here a == b.
2367 movel IMM (EQUAL),d0
2368 #ifndef __mcoldfire__
2369 moveml sp@+,d2-d7 | put back the registers
2372 | XXX if frame pointer is ever removed, stack pointer must
2378 movel IMM (GREATER),d0
2379 #ifndef __mcoldfire__
2380 moveml sp@+,d2-d7 | put back the registers
2383 | XXX if frame pointer is ever removed, stack pointer must
2390 #ifndef __mcoldfire__
2391 moveml sp@+,d2-d7 | put back the registers
2394 | XXX if frame pointer is ever removed, stack pointer must
2419 moveq IMM (INEXACT_RESULT+INVALID_OPERATION),d7
2420 moveq IMM (DOUBLE_FLOAT),d6
2421 PICJUMP $_exception_handler
2423 | int __cmpdf2(double, double);
2432 PICCALL SYM (__cmpdf2_internal)
2436 |=============================================================================
2438 |=============================================================================
2440 | The rounding routines expect the number to be normalized in registers
2441 | d0-d1-d2-d3, with the exponent in register d4. They assume that the
2442 | exponent is larger or equal to 1. They return a properly normalized number
2443 | if possible, and a denormalized number otherwise. The exponent is returned
2447 | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
2448 | Here we assume that the exponent is not too small (this should be checked
2449 | before entering the rounding routine), but the number could be denormalized.
2451 | Check for denormalized numbers:
2452 1: btst IMM (DBL_MANT_DIG-32),d0
2453 bne 2f | if set the number is normalized
2454 | Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent
2455 | is one (remember that a denormalized number corresponds to an
2456 | exponent of -D_BIAS+1).
2457 #ifndef __mcoldfire__
2458 cmpw IMM (1),d4 | remember that the exponent is at least one
2460 cmpl IMM (1),d4 | remember that the exponent is at least one
2462 beq 2f | an exponent of one means denormalized
2463 addl d3,d3 | else shift and adjust the exponent
2467 #ifndef __mcoldfire__
2474 | Now round: we do it as follows: after the shifting we can write the
2475 | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
2476 | If delta < 1, do nothing. If delta > 1, add 1 to f.
2477 | If delta == 1, we make sure the rounded number will be even (odd?)
2479 btst IMM (0),d1 | is delta < 1?
2480 beq 2f | if so, do not do anything
2481 orl d2,d3 | is delta == 1?
2482 bne 1f | if so round to even
2484 andl IMM (2),d3 | bit 1 is the last significant bit
2489 1: movel IMM (1),d3 | else add 1
2493 | Shift right once (because we used bit #DBL_MANT_DIG-32!).
2495 #ifndef __mcoldfire__
2506 | Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
2507 | 'fraction overflow
' ...).
2508 btst IMM (DBL_MANT_DIG-32),d0
2510 #ifndef __mcoldfire__
2523 | If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we
2524 | have to put the exponent to zero and return a denormalized number.
2525 btst IMM (DBL_MANT_DIG-32-1),d0
2535 #endif /* L_double */
2540 .globl $_exception_handler
2542 QUIET_NaN = 0xffffffff
2543 SIGNL_NaN = 0x7f800001
2544 INFINITY = 0x7f800000
2548 FLT_MAX_EXP = F_MAX_EXP - F_BIAS
2549 FLT_MIN_EXP = 1 - F_BIAS
2552 INEXACT_RESULT = 0x0001
2555 DIVIDE_BY_ZERO = 0x0008
2556 INVALID_OPERATION = 0x0010
2570 ROUND_TO_NEAREST = 0 | round result to nearest representable value
2571 ROUND_TO_ZERO = 1 | round result towards zero
2572 ROUND_TO_PLUS = 2 | round result towards plus infinity
2573 ROUND_TO_MINUS = 3 | round result towards minus infinity
2577 .globl SYM (__addsf3)
2578 .globl SYM (__subsf3)
2579 .globl SYM (__mulsf3)
2580 .globl SYM (__divsf3)
2581 .globl SYM (__negsf2)
2582 .globl SYM (__cmpsf2)
2583 .globl SYM (__cmpsf2_internal)
2584 .hidden SYM (__cmpsf2_internal)
2586 | These are common routines to return and signal exceptions.
2592 | Return and signal a denormalized number
2594 moveq IMM (INEXACT_RESULT+UNDERFLOW),d7
2595 moveq IMM (SINGLE_FLOAT),d6
2596 PICJUMP $_exception_handler
2600 | Return a properly signed INFINITY and set the exception flags
2601 movel IMM (INFINITY),d0
2603 moveq IMM (INEXACT_RESULT+OVERFLOW),d7
2604 moveq IMM (SINGLE_FLOAT),d6
2605 PICJUMP $_exception_handler
2608 | Return 0 and set the exception flags
2610 moveq IMM (INEXACT_RESULT+UNDERFLOW),d7
2611 moveq IMM (SINGLE_FLOAT),d6
2612 PICJUMP $_exception_handler
2615 | Return a quiet NaN and set the exception flags
2616 movel IMM (QUIET_NaN),d0
2617 moveq IMM (INEXACT_RESULT+INVALID_OPERATION),d7
2618 moveq IMM (SINGLE_FLOAT),d6
2619 PICJUMP $_exception_handler
2622 | Return a properly signed INFINITY and set the exception flags
2623 movel IMM (INFINITY),d0
2625 moveq IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
2626 moveq IMM (SINGLE_FLOAT),d6
2627 PICJUMP $_exception_handler
2629 |=============================================================================
2630 |=============================================================================
2631 | single precision routines
2632 |=============================================================================
2633 |=============================================================================
2635 | A single precision floating point number (float) has the format:
2638 | unsigned int sign : 1; /* sign bit */
2639 | unsigned int exponent : 8; /* exponent, shifted by 126 */
2640 | unsigned int fraction : 23; /* fraction */
2643 | Thus sizeof(float) = 4 (32 bits).
2645 | All the routines are callable from C programs, and return the result
2646 | in the single register d0. They also preserve all registers except
2649 |=============================================================================
2651 |=============================================================================
2653 | float __subsf3(float, float);
2656 bchg IMM (31),sp@(8) | change sign of second operand
2658 |=============================================================================
2660 |=============================================================================
2662 | float __addsf3(float, float);
2665 #ifndef __mcoldfire__
2666 link a6,IMM (0) | everything will be done in registers
2667 moveml d2-d7,sp@- | save all data registers but d0-d1
2672 movel a6@(8),d0 | get first operand
2673 movel a6@(12),d1 | get second operand
2674 movel d0,a0 | get d0's sign bit
'
2675 addl d0,d0 | check and clear sign bit of a
2676 beq Laddsf$b | if zero return second operand
2677 movel d1,a1 | save b's sign bit
'
2678 addl d1,d1 | get rid of sign bit
2679 beq Laddsf$a | if zero return first operand
2681 | Get the exponents and check for denormalized and/or infinity.
2683 movel IMM (0x00ffffff),d4 | mask to get fraction
2684 movel IMM (0x01000000),d5 | mask to put hidden bit back
2686 movel d0,d6 | save a to get exponent
2687 andl d4,d0 | get fraction in d0
2688 notl d4 | make d4 into a mask for the exponent
2689 andl d4,d6 | get exponent in d6
2690 beq Laddsf$a$den | branch if a is denormalized
2691 cmpl d4,d6 | check for INFINITY or NaN
2693 swap d6 | put exponent into first word
2694 orl d5,d0 | and put hidden bit back
2696 | Now we have a's exponent
in d6
(second
byte) and the mantissa
in d0.
'
2697 movel d1,d7 | get exponent in d7
2699 beq Laddsf$b$den | branch if b is denormalized
2700 cmpl d4,d7 | check for INFINITY or NaN
2702 swap d7 | put exponent into first word
2703 notl d4 | make d4 into a mask for the fraction
2704 andl d4,d1 | get fraction in d1
2705 orl d5,d1 | and put hidden bit back
2707 | Now we have b's exponent
in d7
(second
byte) and the mantissa
in d1.
'
2709 | Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we
2710 | shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
2713 movel d1,d2 | move b to d2, since we want to use
2714 | two registers to do the sum
2715 movel IMM (0),d1 | and clear the new ones
2718 | Here we shift the numbers in registers d0 and d1 so the exponents are the
2719 | same, and put the largest exponent in d6. Note that we are using two
2720 | registers for each number (see the discussion by D. Knuth in "Seminumerical
2722 #ifndef __mcoldfire__
2723 cmpw d6,d7 | compare exponents
2725 cmpl d6,d7 | compare exponents
2727 beq Laddsf$3 | if equal don't shift
'
2728 bhi 5f | branch if second exponent largest
2730 subl d6,d7 | keep the largest exponent
2732 #ifndef __mcoldfire__
2733 lsrw IMM (8),d7 | put difference in lower byte
2735 lsrl IMM (8),d7 | put difference in lower byte
2737 | if difference is too large we don't shift
(actually
, we can just exit
) '
2738 #ifndef __mcoldfire__
2739 cmpw IMM (FLT_MANT_DIG+2),d7
2741 cmpl IMM (FLT_MANT_DIG+2),d7
2744 #ifndef __mcoldfire__
2745 cmpw IMM (16),d7 | if difference >= 16 swap
2747 cmpl IMM (16),d7 | if difference >= 16 swap
2751 #ifndef __mcoldfire__
2757 #ifndef __mcoldfire__
2758 lsrl IMM (1),d2 | shift right second operand
2776 #ifndef __mcoldfire__
2781 bne 2b | if still more bits, go back to normal case
2784 #ifndef __mcoldfire__
2785 exg d6,d7 | exchange the exponents
2791 subl d6,d7 | keep the largest exponent
2793 #ifndef __mcoldfire__
2794 lsrw IMM (8),d7 | put difference in lower byte
2796 lsrl IMM (8),d7 | put difference in lower byte
2798 | if difference is too large we don't shift
(and exit
!) '
2799 #ifndef __mcoldfire__
2800 cmpw IMM (FLT_MANT_DIG+2),d7
2802 cmpl IMM (FLT_MANT_DIG+2),d7
2805 #ifndef __mcoldfire__
2806 cmpw IMM (16),d7 | if difference >= 16 swap
2808 cmpl IMM (16),d7 | if difference >= 16 swap
2812 #ifndef __mcoldfire__
2818 #ifndef __mcoldfire__
2819 lsrl IMM (1),d0 | shift right first operand
2837 #ifndef __mcoldfire__
2842 bne 6b | if still more bits, go back to normal case
2843 | otherwise we fall through
2845 | Now we have a in d0-d1, b in d2-d3, and the largest exponent in d6 (the
2846 | signs are stored in a0 and a1).
2849 | Here we have to decide whether to add or subtract the numbers
2850 #ifndef __mcoldfire__
2851 exg d6,a0 | get signs back
2852 exg d7,a1 | and save the exponents
2861 eorl d6,d7 | combine sign bits
2862 bmi Lsubsf$0 | if negative a and b have opposite
2863 | sign so we actually subtract the
2866 | Here we have both positive or both negative
2867 #ifndef __mcoldfire__
2868 exg d6,a0 | now we have the exponent in d6
2874 movel a0,d7 | and sign in d7
2875 andl IMM (0x80000000),d7
2876 | Here we do the addition.
2879 | Note: now we have d2, d3, d4 and d5 to play with!
2881 | Put the exponent, in the first byte, in d2, to use the "standard" rounding
2884 #ifndef __mcoldfire__
2890 | Before rounding normalize so bit #FLT_MANT_DIG is set (we will consider
2891 | the case of denormalized numbers in the rounding routine itself).
2892 | As in the addition (not in the subtraction!) we could have set
2893 | one more bit we check this:
2894 btst IMM (FLT_MANT_DIG+1),d0
2896 #ifndef __mcoldfire__
2908 lea pc@(Laddsf$4),a0 | to return from rounding routine
2909 PICLEA SYM (_fpCCR),a1 | check the rounding mode
2910 #ifdef __mcoldfire__
2913 movew a1@(6),d6 | rounding mode in d6
2914 beq Lround$to$nearest
2915 #ifndef __mcoldfire__
2916 cmpw IMM (ROUND_TO_PLUS),d6
2918 cmpl IMM (ROUND_TO_PLUS),d6
2924 | Put back the exponent, but check for overflow.
2925 #ifndef __mcoldfire__
2931 bclr IMM (FLT_MANT_DIG-1),d0
2932 #ifndef __mcoldfire__
2945 | We are here if a > 0 and b < 0 (sign bits cleared).
2946 | Here we do the subtraction.
2947 movel d6,d7 | put sign in d7
2948 andl IMM (0x80000000),d7
2950 subl d3,d1 | result in d0-d1
2952 beq Laddsf$ret | if zero just exit
2953 bpl 1f | if positive skip the following
2954 bchg IMM (31),d7 | change sign bit in d7
2958 #ifndef __mcoldfire__
2959 exg d2,a0 | now we have the exponent in d2
2960 lsrw IMM (8),d2 | put it in the first byte
2965 lsrl IMM (8),d2 | put it in the first byte
2968 | Now d0-d1 is positive and the sign bit is in d7.
2970 | Note that we do not have to normalize, since in the subtraction bit
2971 | #FLT_MANT_DIG+1 is never set, and denormalized numbers are handled by
2972 | the rounding routines themselves.
2973 lea pc@(Lsubsf$1),a0 | to return from rounding routine
2974 PICLEA SYM (_fpCCR),a1 | check the rounding mode
2975 #ifdef __mcoldfire__
2978 movew a1@(6),d6 | rounding mode in d6
2979 beq Lround$to$nearest
2980 #ifndef __mcoldfire__
2981 cmpw IMM (ROUND_TO_PLUS),d6
2983 cmpl IMM (ROUND_TO_PLUS),d6
2989 | Put back the exponent (we can't have overflow
!).
'
2990 bclr IMM (FLT_MANT_DIG-1),d0
2991 #ifndef __mcoldfire__
3000 | If one of the numbers was too small (difference of exponents >=
3001 | FLT_MANT_DIG+2) we return the other (and now we don't have to
'
3002 | check for finiteness or zero).
3005 PICLEA SYM (_fpCCR),a0
3007 #ifndef __mcoldfire__
3008 moveml sp@+,d2-d7 | restore data registers
3011 | XXX if frame pointer is ever removed, stack pointer must
3014 unlk a6 | and return
3019 PICLEA SYM (_fpCCR),a0
3021 #ifndef __mcoldfire__
3022 moveml sp@+,d2-d7 | restore data registers
3025 | XXX if frame pointer is ever removed, stack pointer must
3028 unlk a6 | and return
3031 | If the numbers are denormalized remember to put exponent equal to 1.
3034 movel d5,d6 | d5 contains 0x01000000
3041 notl d4 | make d4 into a mask for the fraction
3042 | (this was not executed after the jump)
3045 | The rest is mainly code for the different results which can be
3046 | returned (checking always for +/-INFINITY and NaN).
3049 | Return b (if a is zero).
3051 cmpl IMM (0x80000000),d0 | Check if b is -0
3054 andl IMM (0x80000000),d7 | Use the sign of a
3058 | Return a (if b is zero).
3062 | We have to check for NaN and +/-infty.
3064 andl IMM (0x80000000),d7 | put sign in d7
3065 bclr IMM (31),d0 | clear sign
3066 cmpl IMM (INFINITY),d0 | check for infty or NaN
3068 movel d0,d0 | check for zero (we do this because we don't
'
3069 bne Laddsf$ret | want to return -0 by mistake
3070 bclr IMM (31),d7 | if zero be sure to clear sign
3071 bra Laddsf$ret | if everything OK just return
3073 | The value to be returned is either +/-infty or NaN
3074 andl IMM (0x007fffff),d0 | check for NaN
3075 bne Lf$inop | if mantissa not zero is NaN
3079 | Normal exit (a and b nonzero, result is not NaN nor +/-infty).
3080 | We have to clear the exception flags (just the exception type).
3081 PICLEA SYM (_fpCCR),a0
3083 orl d7,d0 | put sign bit
3084 #ifndef __mcoldfire__
3085 moveml sp@+,d2-d7 | restore data registers
3088 | XXX if frame pointer is ever removed, stack pointer must
3091 unlk a6 | and return
3095 | Return a denormalized number (for addition we don't signal underflow
) '
3096 lsrl IMM (1),d0 | remember to shift right back once
3097 bra Laddsf$ret | and return
3099 | Note: when adding two floats of the same sign if either one is
3100 | NaN we return NaN without regard to whether the other is finite or
3101 | not. When subtracting them (i.e., when adding two numbers of
3102 | opposite signs) things are more complicated: if both are INFINITY
3103 | we return NaN, if only one is INFINITY and the other is NaN we return
3104 | NaN, but if it is finite we return INFINITY with the corresponding sign.
3108 | This could be faster but it is not worth the effort, since it is not
3109 | executed very often. We sacrifice speed for clarity here.
3110 movel a6@(8),d0 | get the numbers back (remember that we
3111 movel a6@(12),d1 | did some processing already)
3112 movel IMM (INFINITY),d4 | useful constant (INFINITY)
3113 movel d0,d2 | save sign bits
3115 bclr IMM (31),d0 | clear sign bits
3117 | We know that one of them is either NaN of +/-INFINITY
3118 | Check for NaN (if either one is NaN return NaN)
3119 cmpl d4,d0 | check first a (d0)
3121 cmpl d4,d1 | check now b (d1)
3123 | Now comes the check for +/-INFINITY. We know that both are (maybe not
3124 | finite) numbers, but we have to check if both are infinite whether we
3125 | are adding or subtracting them.
3126 eorl d3,d2 | to check sign bits
3129 andl IMM (0x80000000),d7 | get (common) sign bit
3132 | We know one (or both) are infinite, so we test for equality between the
3133 | two numbers (if they are equal they have to be infinite both, so we
3135 cmpl d1,d0 | are both infinite?
3136 beq Lf$inop | if so return NaN
3139 andl IMM (0x80000000),d7 | get a's sign bit
'
3140 cmpl d4,d0 | test now for infinity
3141 beq Lf$infty | if a is INFINITY return with this sign
3142 bchg IMM (31),d7 | else we know b is INFINITY and has
3143 bra Lf$infty | the opposite sign
3145 |=============================================================================
3147 |=============================================================================
3149 | float __mulsf3(float, float);
3152 #ifndef __mcoldfire__
3159 movel a6@(8),d0 | get a into d0
3160 movel a6@(12),d1 | and b into d1
3161 movel d0,d7 | d7 will hold the sign of the product
3163 andl IMM (0x80000000),d7
3164 movel IMM (INFINITY),d6 | useful constant (+INFINITY)
3165 movel d6,d5 | another (mask for fraction)
3167 movel IMM (0x00800000),d4 | this is to put hidden bit back
3168 bclr IMM (31),d0 | get rid of a's sign bit
'
3170 beq Lmulsf$a$0 | branch if a is zero
3171 bclr IMM (31),d1 | get rid of b's sign bit
'
3173 beq Lmulsf$b$0 | branch if b is zero
3174 cmpl d6,d0 | is a big?
3175 bhi Lmulsf$inop | if a is NaN return NaN
3176 beq Lmulsf$inf | if a is INFINITY we have to check b
3177 cmpl d6,d1 | now compare b with INFINITY
3178 bhi Lmulsf$inop | is b NaN?
3179 beq Lmulsf$overflow | is b INFINITY?
3180 | Here we have both numbers finite and nonzero (and with no sign bit).
3181 | Now we get the exponents into d2 and d3.
3182 andl d6,d2 | and isolate exponent in d2
3183 beq Lmulsf$a$den | if exponent is zero we have a denormalized
3184 andl d5,d0 | and isolate fraction
3185 orl d4,d0 | and put hidden bit back
3186 swap d2 | I like exponents in the first byte
3187 #ifndef __mcoldfire__
3198 #ifndef __mcoldfire__
3204 #ifndef __mcoldfire__
3205 addw d3,d2 | add exponents
3206 subw IMM (F_BIAS+1),d2 | and subtract bias (plus one)
3208 addl d3,d2 | add exponents
3209 subl IMM (F_BIAS+1),d2 | and subtract bias (plus one)
3212 | We are now ready to do the multiplication. The situation is as follows:
3213 | both a and b have bit FLT_MANT_DIG-1 set (even if they were
3214 | denormalized to start with!), which means that in the product
3215 | bit 2*(FLT_MANT_DIG-1) (that is, bit 2*FLT_MANT_DIG-2-32 of the
3216 | high long) is set.
3218 | To do the multiplication let us move the number a little bit around ...
3219 movel d1,d6 | second operand in d6
3220 movel d0,d5 | first operand in d4-d5
3222 movel d4,d1 | the sums will go in d0-d1
3225 | now bit FLT_MANT_DIG-1 becomes bit 31:
3226 lsll IMM (31-FLT_MANT_DIG+1),d6
3228 | Start the loop (we loop #FLT_MANT_DIG times):
3229 moveq IMM (FLT_MANT_DIG-1),d3
3230 1: addl d1,d1 | shift sum
3232 lsll IMM (1),d6 | get bit bn
3233 bcc 2f | if not set skip sum
3237 #ifndef __mcoldfire__
3238 dbf d3,1b | loop back
3244 | Now we have the product in d0-d1, with bit (FLT_MANT_DIG - 1) + FLT_MANT_DIG
3245 | (mod 32) of d0 set. The first thing to do now is to normalize it so bit
3246 | FLT_MANT_DIG is set (to do the rounding).
3247 #ifndef __mcoldfire__
3251 andw IMM (0x03ff),d3
3252 andw IMM (0xfd00),d1
3261 andl IMM (0xfffffd00),d1
3266 #ifndef __mcoldfire__
3272 moveq IMM (MULTIPLY),d5
3274 btst IMM (FLT_MANT_DIG+1),d0
3276 #ifndef __mcoldfire__
3291 moveq IMM (MULTIPLY),d5
3295 moveq IMM (MULTIPLY),d5
3299 moveq IMM (MULTIPLY),d5
3300 | If either is NaN return NaN; else both are (maybe infinite) numbers, so
3301 | return INFINITY with the correct sign (which is in d7).
3302 cmpl d6,d1 | is b NaN?
3303 bhi Lf$inop | if so return NaN
3304 bra Lf$overflow | else return +/-INFINITY
3306 | If either number is zero return zero, unless the other is +/-INFINITY,
3307 | or NaN, in which case we return NaN.
3309 | Here d1 (==b) is zero.
3310 movel a6@(8),d1 | get a again to check for non-finiteness
3313 movel a6@(12),d1 | get b again to check for non-finiteness
3314 1: bclr IMM (31),d1 | clear sign bit
3315 cmpl IMM (INFINITY),d1 | and check for a large exponent
3316 bge Lf$inop | if b is +/-INFINITY or NaN return NaN
3317 movel d7,d0 | else return signed zero
3318 PICLEA SYM (_fpCCR),a0 |
3320 #ifndef __mcoldfire__
3324 | XXX if frame pointer is ever removed, stack pointer must
3330 | If a number is denormalized we put an exponent of 1 but do not put the
3331 | hidden bit back into the fraction; instead we shift left until bit 23
3332 | (the hidden bit) is set, adjusting the exponent accordingly. We do this
3333 | to ensure that the product of the fractions is close to 1.
3337 1: addl d0,d0 | shift a left (until bit 23 is set)
3338 #ifndef __mcoldfire__
3339 subw IMM (1),d2 | and adjust exponent
3341 subql IMM (1),d2 | and adjust exponent
3343 btst IMM (FLT_MANT_DIG-1),d0
3345 bra 1b | else loop back
3350 1: addl d1,d1 | shift b left until bit 23 is set
3351 #ifndef __mcoldfire__
3352 subw IMM (1),d3 | and adjust exponent
3354 subql IMM (1),d3 | and adjust exponent
3356 btst IMM (FLT_MANT_DIG-1),d1
3358 bra 1b | else loop back
3360 |=============================================================================
3362 |=============================================================================
3364 | float __divsf3(float, float);
3367 #ifndef __mcoldfire__
3374 movel a6@(8),d0 | get a into d0
3375 movel a6@(12),d1 | and b into d1
3376 movel d0,d7 | d7 will hold the sign of the result
3378 andl IMM (0x80000000),d7 |
3379 movel IMM (INFINITY),d6 | useful constant (+INFINITY)
3380 movel d6,d5 | another (mask for fraction)
3382 movel IMM (0x00800000),d4 | this is to put hidden bit back
3383 bclr IMM (31),d0 | get rid of a's sign bit
'
3385 beq Ldivsf$a$0 | branch if a is zero
3386 bclr IMM (31),d1 | get rid of b's sign bit
'
3388 beq Ldivsf$b$0 | branch if b is zero
3389 cmpl d6,d0 | is a big?
3390 bhi Ldivsf$inop | if a is NaN return NaN
3391 beq Ldivsf$inf | if a is INFINITY we have to check b
3392 cmpl d6,d1 | now compare b with INFINITY
3393 bhi Ldivsf$inop | if b is NaN return NaN
3394 beq Ldivsf$underflow
3395 | Here we have both numbers finite and nonzero (and with no sign bit).
3396 | Now we get the exponents into d2 and d3 and normalize the numbers to
3397 | ensure that the ratio of the fractions is close to 1. We do this by
3398 | making sure that bit #FLT_MANT_DIG-1 (hidden bit) is set.
3399 andl d6,d2 | and isolate exponent in d2
3400 beq Ldivsf$a$den | if exponent is zero we have a denormalized
3401 andl d5,d0 | and isolate fraction
3402 orl d4,d0 | and put hidden bit back
3403 swap d2 | I like exponents in the first byte
3404 #ifndef __mcoldfire__
3415 #ifndef __mcoldfire__
3421 #ifndef __mcoldfire__
3422 subw d3,d2 | subtract exponents
3423 addw IMM (F_BIAS),d2 | and add bias
3425 subl d3,d2 | subtract exponents
3426 addl IMM (F_BIAS),d2 | and add bias
3429 | We are now ready to do the division. We have prepared things in such a way
3430 | that the ratio of the fractions will be less than 2 but greater than 1/2.
3431 | At this point the registers in use are:
3432 | d0 holds a (first operand, bit FLT_MANT_DIG=0, bit FLT_MANT_DIG-1=1)
3433 | d1 holds b (second operand, bit FLT_MANT_DIG=1)
3434 | d2 holds the difference of the exponents, corrected by the bias
3435 | d7 holds the sign of the ratio
3436 | d4, d5, d6 hold some constants
3437 movel d7,a0 | d6-d7 will hold the ratio of the fractions
3441 moveq IMM (FLT_MANT_DIG+1),d3
3442 1: cmpl d0,d1 | is a < b?
3444 bset d3,d6 | set a bit in d6
3445 subl d1,d0 | if a >= b a <-- a-b
3446 beq 3f | if a is zero, exit
3447 2: addl d0,d0 | multiply a by 2
3448 #ifndef __mcoldfire__
3455 | Now we keep going to set the sticky bit ...
3456 moveq IMM (FLT_MANT_DIG),d3
3460 #ifndef __mcoldfire__
3469 #ifndef __mcoldfire__
3470 subw IMM (FLT_MANT_DIG),d3
3473 subl IMM (FLT_MANT_DIG),d3
3478 movel d6,d0 | put the ratio in d0-d1
3479 movel a0,d7 | get sign back
3481 | Because of the normalization we did before we are guaranteed that
3482 | d0 is smaller than 2^26 but larger than 2^24. Thus bit 26 is not set,
3483 | bit 25 could be set, and if it is not set then bit 24 is necessarily set.
3484 btst IMM (FLT_MANT_DIG+1),d0
3485 beq 1f | if it is not set, then bit 24 is set
3487 #ifndef __mcoldfire__
3493 | Now round, check for over- and underflow, and exit.
3494 moveq IMM (DIVIDE),d5
3498 moveq IMM (DIVIDE),d5
3502 moveq IMM (DIVIDE),d5
3506 moveq IMM (DIVIDE),d5
3510 moveq IMM (DIVIDE),d5
3511 | If a is zero check to see whether b is zero also. In that case return
3512 | NaN; then check if b is NaN, and return NaN also in that case. Else
3513 | return a properly signed zero.
3514 andl IMM (0x7fffffff),d1 | clear sign bit and test b
3515 beq Lf$inop | if b is also zero return NaN
3516 cmpl IMM (INFINITY),d1 | check for NaN
3518 movel d7,d0 | else return signed zero
3519 PICLEA SYM (_fpCCR),a0 |
3521 #ifndef __mcoldfire__
3525 | XXX if frame pointer is ever removed, stack pointer must
3532 moveq IMM (DIVIDE),d5
3533 | If we got here a is not zero. Check if a is NaN; in that case return NaN,
3534 | else return +/-INFINITY. Remember that a is in d0 with the sign bit
3536 cmpl IMM (INFINITY),d0 | compare d0 with INFINITY
3537 bhi Lf$inop | if larger it is NaN
3538 bra Lf$div$0 | else signal DIVIDE_BY_ZERO
3541 moveq IMM (DIVIDE),d5
3542 | If a is INFINITY we have to check b
3543 cmpl IMM (INFINITY),d1 | compare b with INFINITY
3544 bge Lf$inop | if b is NaN or INFINITY return NaN
3545 bra Lf$overflow | else return overflow
3547 | If a number is denormalized we put an exponent of 1 but do not put the
3548 | bit back into the fraction.
3552 1: addl d0,d0 | shift a left until bit FLT_MANT_DIG-1 is set
3553 #ifndef __mcoldfire__
3554 subw IMM (1),d2 | and adjust exponent
3556 subl IMM (1),d2 | and adjust exponent
3558 btst IMM (FLT_MANT_DIG-1),d0
3565 1: addl d1,d1 | shift b left until bit FLT_MANT_DIG is set
3566 #ifndef __mcoldfire__
3567 subw IMM (1),d3 | and adjust exponent
3569 subl IMM (1),d3 | and adjust exponent
3571 btst IMM (FLT_MANT_DIG-1),d1
3576 | This is a common exit point for __mulsf3 and __divsf3.
3578 | First check for underlow in the exponent:
3579 #ifndef __mcoldfire__
3580 cmpw IMM (-FLT_MANT_DIG-1),d2
3582 cmpl IMM (-FLT_MANT_DIG-1),d2
3585 | It could happen that the exponent is less than 1, in which case the
3586 | number is denormalized. In this case we shift right and adjust the
3587 | exponent until it becomes 1 or the fraction is zero (in the latter case
3588 | we signal underflow and return zero).
3589 movel IMM (0),d6 | d6 is used temporarily
3590 #ifndef __mcoldfire__
3591 cmpw IMM (1),d2 | if the exponent is less than 1 we
3593 cmpl IMM (1),d2 | if the exponent is less than 1 we
3595 bge 2f | have to shift right (denormalize)
3597 #ifndef __mcoldfire__
3598 addw IMM (1),d2 | adjust the exponent
3599 lsrl IMM (1),d0 | shift right once
3601 roxrl IMM (1),d6 | d6 collect bits we would lose otherwise
3602 cmpw IMM (1),d2 | is the exponent 1 already?
3604 addql IMM (1),d2 | adjust the exponent
3614 cmpl IMM (1),d2 | is the exponent 1 already?
3616 beq 2f | if not loop back
3618 bra Lf$underflow | safety check, shouldn't execute
'
3619 2: orl d6,d1 | this is a trick so we don't lose
'
3620 | the extra bits which were flushed right
3621 | Now call the rounding routine (which takes care of denormalized numbers):
3622 lea pc@(Lround$0),a0 | to return from rounding routine
3623 PICLEA SYM (_fpCCR),a1 | check the rounding mode
3624 #ifdef __mcoldfire__
3627 movew a1@(6),d6 | rounding mode in d6
3628 beq Lround$to$nearest
3629 #ifndef __mcoldfire__
3630 cmpw IMM (ROUND_TO_PLUS),d6
3632 cmpl IMM (ROUND_TO_PLUS),d6
3638 | Here we have a correctly rounded result (either normalized or denormalized).
3640 | Here we should have either a normalized number or a denormalized one, and
3641 | the exponent is necessarily larger or equal to 1 (so we don't have to
'
3642 | check again for underflow!). We have to check for overflow or for a
3643 | denormalized number (which also signals underflow).
3644 | Check for overflow (i.e., exponent >= 255).
3645 #ifndef __mcoldfire__
3646 cmpw IMM (0x00ff),d2
3648 cmpl IMM (0x00ff),d2
3651 | Now check for a denormalized number (exponent==0).
3655 | Put back the exponents and sign and return.
3656 #ifndef __mcoldfire__
3657 lslw IMM (7),d2 | exponent back to fourth byte
3659 lsll IMM (7),d2 | exponent back to fourth byte
3661 bclr IMM (FLT_MANT_DIG-1),d0
3662 swap d0 | and put back exponent
3663 #ifndef __mcoldfire__
3669 orl d7,d0 | and sign also
3671 PICLEA SYM (_fpCCR),a0
3673 #ifndef __mcoldfire__
3677 | XXX if frame pointer is ever removed, stack pointer must
3683 |=============================================================================
3685 |=============================================================================
3687 | This is trivial and could be shorter if we didn't bother checking for NaN
'
3690 | float __negsf2(float);
3693 #ifndef __mcoldfire__
3700 moveq IMM (NEGATE),d5
3701 movel a6@(8),d0 | get number to negate in d0
3702 bchg IMM (31),d0 | negate
3703 movel d0,d1 | make a positive copy
3705 tstl d1 | check for zero
3706 beq 2f | if zero (either sign) return +zero
3707 cmpl IMM (INFINITY),d1 | compare to +INFINITY
3709 bhi Lf$inop | if larger (fraction not zero) is NaN
3710 movel d0,d7 | else get sign and return INFINITY
3711 andl IMM (0x80000000),d7
3713 1: PICLEA SYM (_fpCCR),a0
3715 #ifndef __mcoldfire__
3719 | XXX if frame pointer is ever removed, stack pointer must
3727 |=============================================================================
3729 |=============================================================================
3735 | int __cmpsf2_internal(float, float, int);
3736 SYM (__cmpsf2_internal):
3737 #ifndef __mcoldfire__
3739 moveml d2-d7,sp@- | save registers
3744 moveq IMM (COMPARE),d5
3745 movel a6@(8),d0 | get first operand
3746 movel a6@(12),d1 | get second operand
3747 | Check if either is NaN, and in that case return garbage and signal
3748 | INVALID_OPERATION. Check also if either is zero, and clear the signs
3751 andl IMM (0x7fffffff),d0
3753 cmpl IMM (0x7f800000),d0
3757 andl IMM (0x7fffffff),d1
3759 cmpl IMM (0x7f800000),d1
3765 | If the signs are not equal check if a >= 0
3767 bpl Lcmpsf$a$gt$b | if (a >= 0 && b < 0) => a > b
3768 bmi Lcmpsf$b$gt$a | if (a < 0 && b >= 0) => a < b
3770 | If the signs are equal check for < 0
3773 | If both are negative exchange them
3774 #ifndef __mcoldfire__
3782 | Now that they are positive we just compare them as longs (does this also
3783 | work for denormalized numbers?).
3785 bhi Lcmpsf$b$gt$a | |b| > |a|
3786 bne Lcmpsf$a$gt$b | |b| < |a|
3787 | If we got here a == b.
3788 movel IMM (EQUAL),d0
3789 #ifndef __mcoldfire__
3790 moveml sp@+,d2-d7 | put back the registers
3797 movel IMM (GREATER),d0
3798 #ifndef __mcoldfire__
3799 moveml sp@+,d2-d7 | put back the registers
3802 | XXX if frame pointer is ever removed, stack pointer must
3809 #ifndef __mcoldfire__
3810 moveml sp@+,d2-d7 | put back the registers
3813 | XXX if frame pointer is ever removed, stack pointer must
3828 moveq IMM (INEXACT_RESULT+INVALID_OPERATION),d7
3829 moveq IMM (SINGLE_FLOAT),d6
3830 PICJUMP $_exception_handler
3832 | int __cmpsf2(float, float);
3839 PICCALL SYM (__cmpsf2_internal)
3843 |=============================================================================
3845 |=============================================================================
3847 | The rounding routines expect the number to be normalized in registers
3848 | d0-d1, with the exponent in register d2. They assume that the
3849 | exponent is larger or equal to 1. They return a properly normalized number
3850 | if possible, and a denormalized number otherwise. The exponent is returned
3854 | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
3855 | Here we assume that the exponent is not too small (this should be checked
3856 | before entering the rounding routine), but the number could be denormalized.
3858 | Check for denormalized numbers:
3859 1: btst IMM (FLT_MANT_DIG),d0
3860 bne 2f | if set the number is normalized
3861 | Normalize shifting left until bit #FLT_MANT_DIG is set or the exponent
3862 | is one (remember that a denormalized number corresponds to an
3863 | exponent of -F_BIAS+1).
3864 #ifndef __mcoldfire__
3865 cmpw IMM (1),d2 | remember that the exponent is at least one
3867 cmpl IMM (1),d2 | remember that the exponent is at least one
3869 beq 2f | an exponent of one means denormalized
3870 addl d1,d1 | else shift and adjust the exponent
3872 #ifndef __mcoldfire__
3879 | Now round: we do it as follows: after the shifting we can write the
3880 | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
3881 | If delta < 1, do nothing. If delta > 1, add 1 to f.
3882 | If delta == 1, we make sure the rounded number will be even (odd?)
3884 btst IMM (0),d0 | is delta < 1?
3885 beq 2f | if so, do not do anything
3886 tstl d1 | is delta == 1?
3887 bne 1f | if so round to even
3889 andl IMM (2),d1 | bit 1 is the last significant bit
3892 1: movel IMM (1),d1 | else add 1
3894 | Shift right once (because we used bit #FLT_MANT_DIG!).
3896 | Now check again bit #FLT_MANT_DIG (rounding could have produced a
3897 | 'fraction overflow
' ...).
3898 btst IMM (FLT_MANT_DIG),d0
3901 #ifndef __mcoldfire__
3907 | If bit #FLT_MANT_DIG-1 is clear we have a denormalized number, so we
3908 | have to put the exponent to zero and return a denormalized number.
3909 btst IMM (FLT_MANT_DIG-1),d0
3919 #endif /* L_float */
3921 | gcc expects the routines __eqdf2, __nedf2, __gtdf2, __gedf2,
3922 | __ledf2, __ltdf2 to all return the same value as a direct call to
3923 | __cmpdf2 would. In this implementation, each of these routines
3924 | simply calls __cmpdf2. It would be more efficient to give the
3925 | __cmpdf2 routine several names, but separating them out will make it
3926 | easier to write efficient versions of these routines someday.
3927 | If the operands recompare unordered unordered __gtdf2 and __gedf2 return -1.
3928 | The other routines return 1.
3933 .globl SYM (__eqdf2)
3941 PICCALL SYM (__cmpdf2_internal)
3944 #endif /* L_eqdf2 */
3949 .globl SYM (__nedf2)
3957 PICCALL SYM (__cmpdf2_internal)
3960 #endif /* L_nedf2 */
3965 .globl SYM (__gtdf2)
3973 PICCALL SYM (__cmpdf2_internal)
3976 #endif /* L_gtdf2 */
3981 .globl SYM (__gedf2)
3989 PICCALL SYM (__cmpdf2_internal)
3992 #endif /* L_gedf2 */
3997 .globl SYM (__ltdf2)
4005 PICCALL SYM (__cmpdf2_internal)
4008 #endif /* L_ltdf2 */
4013 .globl SYM (__ledf2)
4021 PICCALL SYM (__cmpdf2_internal)
4024 #endif /* L_ledf2 */
4026 | The comments above about __eqdf2, et. al., also apply to __eqsf2,
4027 | et. al., except that the latter call __cmpsf2 rather than __cmpdf2.
4032 .globl SYM (__eqsf2)
4038 PICCALL SYM (__cmpsf2_internal)
4041 #endif /* L_eqsf2 */
4046 .globl SYM (__nesf2)
4052 PICCALL SYM (__cmpsf2_internal)
4055 #endif /* L_nesf2 */
4060 .globl SYM (__gtsf2)
4066 PICCALL SYM (__cmpsf2_internal)
4069 #endif /* L_gtsf2 */
4074 .globl SYM (__gesf2)
4080 PICCALL SYM (__cmpsf2_internal)
4083 #endif /* L_gesf2 */
4088 .globl SYM (__ltsf2)
4094 PICCALL SYM (__cmpsf2_internal)
4097 #endif /* L_ltsf2 */
4102 .globl SYM (__lesf2)
4108 PICCALL SYM (__cmpsf2_internal)
4111 #endif /* L_lesf2 */
4113 #if defined (__ELF__) && defined (__linux__)
4114 /* Make stack non-executable for ELF linux targets. */
4115 .section .note.GNU-stack,"",@progbits