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