| blob | 22f2772520cf099b0a57c21bccd81aef612fad5c |
1 /* libgcc routines for 68000 w/o floating-point hardware.
2 Copyright (C) 1994-2024 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
9 later version.
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
29 */
31 /* These are predefined by new versions of GNU cpp. */
33 #ifndef __USER_LABEL_PREFIX__
34 #define __USER_LABEL_PREFIX__ _
35 #endif
37 #ifndef __REGISTER_PREFIX__
38 #define __REGISTER_PREFIX__
39 #endif
41 #ifndef __IMMEDIATE_PREFIX__
42 #define __IMMEDIATE_PREFIX__ #
43 #endif
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. */
56 #ifdef __ELF__
57 #define FUNC(x) .type SYM(x),function
58 #else
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
62 the need arises. */
63 #define FUNC(x) .proc
64 #endif
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)
74 #define d0 REG (d0)
75 #define d1 REG (d1)
76 #define d2 REG (d2)
77 #define d3 REG (d3)
78 #define d4 REG (d4)
79 #define d5 REG (d5)
80 #define d6 REG (d6)
81 #define d7 REG (d7)
82 #define a0 REG (a0)
83 #define a1 REG (a1)
84 #define a2 REG (a2)
85 #define a3 REG (a3)
86 #define a4 REG (a4)
87 #define a5 REG (a5)
88 #define a6 REG (a6)
89 #define fp REG (fp)
90 #define sp REG (sp)
91 #define pc REG (pc)
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).
99 */
100 #ifndef __PIC__
102 /* Non PIC (absolute/relocatable) versions */
104 .macro PICCALL addr
105 jbsr \addr
106 .endm
108 .macro PICJUMP addr
109 jmp \addr
110 .endm
112 .macro PICLEA sym, reg
113 lea \sym, \reg
114 .endm
116 .macro PICPEA sym, areg
117 pea \sym
118 .endm
120 #else /* __PIC__ */
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
133 .endm
135 .macro PICPEA sym, areg
136 movel a5@(_current_shared_library_a5_offset_), \areg
137 movel \sym@GOT(\areg), sp@-
138 .endm
140 .macro PICCALL addr
141 PICLEA \addr,a0
142 jsr a0@
143 .endm
145 .macro PICJUMP addr
146 PICLEA \addr,a0
147 jmp a0@
148 .endm
150 # else /* !__ID_SHARED_LIBRARY__ */
152 /* Versions for -msep-data */
154 .macro PICLEA sym, reg
155 movel \sym@GOT(a5), \reg
156 .endm
158 .macro PICPEA sym, areg
159 movel \sym@GOT(a5), sp@-
160 .endm
162 .macro PICCALL addr
163 #if defined (__mcoldfire__) && !defined (__mcfisab__) && !defined (__mcfisac__)
164 lea \addr-.-8,a0
165 jsr pc@(a0)
166 #else
167 jbsr \addr
168 #endif
169 .endm
171 .macro PICJUMP addr
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__)
176 lea \addr-.-8,a0
177 jmp pc@(a0)
178 #else
179 bra \addr
180 #endif
181 .endm
183 # endif
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
193 .endm
195 .macro PICPEA sym, areg
196 movel #_GLOBAL_OFFSET_TABLE_@GOTPC, \areg
197 lea (-6, pc, \areg), \areg
198 movel \sym@GOT(\areg), sp@-
199 .endm
201 .macro PICCALL addr
202 #if defined (__mcoldfire__) && !defined (__mcfisab__) && !defined (__mcfisac__)
203 lea \addr-.-8,a0
204 jsr pc@(a0)
205 #else
206 jbsr \addr
207 #endif
208 .endm
210 .macro PICJUMP addr
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__)
215 lea \addr-.-8,a0
216 jmp pc@(a0)
217 #else
218 bra \addr
219 #endif
220 .endm
222 # endif
223 #endif /* __PIC__ */
226 #ifdef L_floatex
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
239 | different cases:
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 |=============================================================================
249 | exceptions
250 |=============================================================================
252 | This is the floating point condition code register (_fpCCR):
253 |
254 | struct {
255 | short _exception_bits;
256 | short _trap_enable_bits;
257 | short _sticky_bits;
258 | short _rounding_mode;
259 | short _format;
260 | short _last_operation;
261 | union {
262 | float sf;
263 | double df;
264 | } _operand1;
265 | union {
266 | float sf;
267 | double df;
268 | } _operand2;
269 | } _fpCCR;
271 .data
272 .even
274 .globl SYM (_fpCCR)
276 SYM (_fpCCR):
277 __exception_bits:
278 .word 0
279 __trap_enable_bits:
280 .word 0
281 __sticky_bits:
282 .word 0
283 __rounding_mode:
284 .word ROUND_TO_NEAREST
285 __format:
286 .word NIL
287 __last_operation:
288 .word NOOP
289 __operand1:
290 .long 0
291 .long 0
292 __operand2:
293 .long 0
294 .long 0
296 | Offsets:
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
308 UNDERFLOW = 0x0002
309 OVERFLOW = 0x0004
310 DIVIDE_BY_ZERO = 0x0008
311 INVALID_OPERATION = 0x0010
313 | The allowed rounding modes are:
314 UNKNOWN = -1
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:
321 NIL = 0
322 SINGLE_FLOAT = 1
323 DOUBLE_FLOAT = 2
324 LONG_FLOAT = 3
326 | The allowed values for the last operation are:
327 NOOP = 0
328 ADD = 1
329 MULTIPLY = 2
330 DIVIDE = 3
331 NEGATE = 4
332 COMPARE = 5
333 EXTENDSFDF = 6
334 TRUNCDFSF = 7
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)
347 .text
348 .even
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)
355 #else
356 clr.w a0@(STICK)
357 #endif
358 rts
360 |=============================================================================
361 | $_exception_handler
362 |=============================================================================
364 .globl $_exception_handler
366 .text
367 .even
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).
382 FPTRAP = 15
384 $_exception_handler:
385 PICLEA SYM (_fpCCR),a0
386 movew d7,a0@(EBITS) | set __exception_bits
387 #ifndef __mcoldfire__
388 orw d7,a0@(STICK) | and __sticky_bits
389 #else
390 movew a0@(STICK),d4
391 orl d7,d4
392 movew d4,a0@(STICK)
393 #endif
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
400 #else
401 cmpl IMM (SINGLE_FLOAT),d6
402 #endif
403 beq 1f
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)
408 bra 2f
409 1: movel a6@(8),a0@(OPER1)
410 movel a6@(12),a0@(OPER2)
411 2:
412 | And check whether the exception is trap-enabled:
413 #ifndef __mcoldfire__
414 andw a0@(TRAPE),d7 | is exception trap-enabled?
415 #else
416 clrl d6
417 movew a0@(TRAPE),d6
418 andl d6,d7
419 #endif
420 beq 1f | no, exit
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
425 #else
426 1: moveml sp@,d2-d7
427 | XXX if frame pointer is ever removed, stack pointer must
428 | be adjusted here.
429 #endif
430 unlk a6 | and return
431 rts
432 #endif /* L_floatex */
434 #ifdef L_mulsi3
435 .text
436 FUNC(__mulsi3)
437 .globl SYM (__mulsi3)
438 .globl SYM (__mulsi3_internal)
439 .hidden SYM (__mulsi3_internal)
440 SYM (__mulsi3):
441 SYM (__mulsi3_internal):
442 movew sp@(4), d0 /* x0 -> d0 */
443 muluw sp@(10), d0 /* x0*y1 */
444 movew sp@(6), d1 /* x1 -> d1 */
445 muluw sp@(8), d1 /* x1*y0 */
446 #ifndef __mcoldfire__
447 addw d1, d0
448 #else
449 addl d1, d0
450 #endif
451 swap d0
452 clrw d0
453 movew sp@(6), d1 /* x1 -> d1 */
454 muluw sp@(10), d1 /* x1*y1 */
455 addl d1, d0
457 rts
458 #endif /* L_mulsi3 */
460 #ifdef L_udivsi3
461 .text
462 FUNC(__udivsi3)
463 .globl SYM (__udivsi3)
464 .globl SYM (__udivsi3_internal)
465 .hidden SYM (__udivsi3_internal)
466 SYM (__udivsi3):
467 SYM (__udivsi3_internal):
468 #ifndef __mcoldfire__
469 movel d2, sp@-
470 movel sp@(12), d1 /* d1 = divisor */
471 movel sp@(8), d0 /* d0 = dividend */
473 cmpl IMM (0x10000), d1 /* divisor >= 2 ^ 16 ? */
474 jcc L3 /* then try next algorithm */
475 movel d0, d2
476 clrw d2
477 swap d2
478 divu d1, d2 /* high quotient in lower word */
479 movew d2, d0 /* save high quotient */
480 swap d0
481 movew sp@(10), d2 /* get low dividend + high rest */
482 divu d1, d2 /* low quotient */
483 movew d2, d0
484 jra L6
486 L3: movel d1, d2 /* use d2 as divisor backup */
487 L4: lsrl IMM (1), d1 /* shift divisor */
488 lsrl IMM (1), d0 /* shift dividend */
489 cmpl IMM (0x10000), d1 /* still divisor >= 2 ^ 16 ? */
490 jcc L4
491 divu d1, d0 /* now we have 16-bit divisor */
492 andl IMM (0xffff), d0 /* mask out divisor, ignore remainder */
494 /* Multiply the 16-bit tentative quotient with the 32-bit divisor. Because of
495 the operand ranges, this might give a 33-bit product. If this product is
496 greater than the dividend, the tentative quotient was too large. */
497 movel d2, d1
498 mulu d0, d1 /* low part, 32 bits */
499 swap d2
500 mulu d0, d2 /* high part, at most 17 bits */
501 swap d2 /* align high part with low part */
502 tstw d2 /* high part 17 bits? */
503 jne L5 /* if 17 bits, quotient was too large */
504 addl d2, d1 /* add parts */
505 jcs L5 /* if sum is 33 bits, quotient was too large */
506 cmpl sp@(8), d1 /* compare the sum with the dividend */
507 jls L6 /* if sum > dividend, quotient was too large */
508 L5: subql IMM (1), d0 /* adjust quotient */
510 L6: movel sp@+, d2
511 rts
513 #else /* __mcoldfire__ */
515 /* ColdFire implementation of non-restoring division algorithm from
516 Hennessy & Patterson, Appendix A. */
517 link a6,IMM (-12)
518 moveml d2-d4,sp@
519 movel a6@(8),d0
520 movel a6@(12),d1
521 clrl d2 | clear p
522 moveq IMM (31),d4
523 L1: addl d0,d0 | shift reg pair (p,a) one bit left
524 addxl d2,d2
525 movl d2,d3 | subtract b from p, store in tmp.
526 subl d1,d3
527 jcs L2 | if no carry,
528 bset IMM (0),d0 | set the low order bit of a to 1,
529 movl d3,d2 | and store tmp in p.
530 L2: subql IMM (1),d4
531 jcc L1
532 moveml sp@,d2-d4 | restore data registers
533 unlk a6 | and return
534 rts
535 #endif /* __mcoldfire__ */
537 #endif /* L_udivsi3 */
539 #ifdef L_divsi3
540 .text
541 FUNC(__divsi3)
542 .globl SYM (__divsi3)
543 .globl SYM (__divsi3_internal)
544 .hidden SYM (__divsi3_internal)
545 SYM (__divsi3):
546 SYM (__divsi3_internal):
547 movel d2, sp@-
549 moveq IMM (1), d2 /* sign of result stored in d2 (=1 or =-1) */
550 movel sp@(12), d1 /* d1 = divisor */
551 jpl L1
552 negl d1
553 #ifndef __mcoldfire__
554 negb d2 /* change sign because divisor <0 */
555 #else
556 negl d2 /* change sign because divisor <0 */
557 #endif
558 L1: movel sp@(8), d0 /* d0 = dividend */
559 jpl L2
560 negl d0
561 #ifndef __mcoldfire__
562 negb d2
563 #else
564 negl d2
565 #endif
567 L2: movel d1, sp@-
568 movel d0, sp@-
569 PICCALL SYM (__udivsi3_internal) /* divide abs(dividend) by abs(divisor) */
570 addql IMM (8), sp
572 tstb d2
573 jpl L3
574 negl d0
576 L3: movel sp@+, d2
577 rts
578 #endif /* L_divsi3 */
580 #ifdef L_umodsi3
581 .text
582 FUNC(__umodsi3)
583 .globl SYM (__umodsi3)
584 SYM (__umodsi3):
585 movel sp@(8), d1 /* d1 = divisor */
586 movel sp@(4), d0 /* d0 = dividend */
587 movel d1, sp@-
588 movel d0, sp@-
589 PICCALL SYM (__udivsi3_internal)
590 addql IMM (8), sp
591 movel sp@(8), d1 /* d1 = divisor */
592 #ifndef __mcoldfire__
593 movel d1, sp@-
594 movel d0, sp@-
595 PICCALL SYM (__mulsi3_internal) /* d0 = (a/b)*b */
596 addql IMM (8), sp
597 #else
598 mulsl d1,d0
599 #endif
600 movel sp@(4), d1 /* d1 = dividend */
601 subl d0, d1 /* d1 = a - (a/b)*b */
602 movel d1, d0
603 rts
604 #endif /* L_umodsi3 */
606 #ifdef L_modsi3
607 .text
608 FUNC(__modsi3)
609 .globl SYM (__modsi3)
610 SYM (__modsi3):
611 movel sp@(8), d1 /* d1 = divisor */
612 movel sp@(4), d0 /* d0 = dividend */
613 movel d1, sp@-
614 movel d0, sp@-
615 PICCALL SYM (__divsi3_internal)
616 addql IMM (8), sp
617 movel sp@(8), d1 /* d1 = divisor */
618 #ifndef __mcoldfire__
619 movel d1, sp@-
620 movel d0, sp@-
621 PICCALL SYM (__mulsi3_internal) /* d0 = (a/b)*b */
622 addql IMM (8), sp
623 #else
624 mulsl d1,d0
625 #endif
626 movel sp@(4), d1 /* d1 = dividend */
627 subl d0, d1 /* d1 = a - (a/b)*b */
628 movel d1, d0
629 rts
630 #endif /* L_modsi3 */
633 #ifdef L_double
635 .globl SYM (_fpCCR)
636 .globl $_exception_handler
638 QUIET_NaN = 0xffffffff
640 D_MAX_EXP = 0x07ff
641 D_BIAS = 1022
642 DBL_MAX_EXP = D_MAX_EXP - D_BIAS
643 DBL_MIN_EXP = 1 - D_BIAS
644 DBL_MANT_DIG = 53
646 INEXACT_RESULT = 0x0001
647 UNDERFLOW = 0x0002
648 OVERFLOW = 0x0004
649 DIVIDE_BY_ZERO = 0x0008
650 INVALID_OPERATION = 0x0010
652 DOUBLE_FLOAT = 2
654 NOOP = 0
655 ADD = 1
656 MULTIPLY = 2
657 DIVIDE = 3
658 NEGATE = 4
659 COMPARE = 5
660 EXTENDSFDF = 6
661 TRUNCDFSF = 7
663 UNKNOWN = -1
664 ROUND_TO_NEAREST = 0 | round result to nearest representable value
665 ROUND_TO_ZERO = 1 | round result towards zero
666 ROUND_TO_PLUS = 2 | round result towards plus infinity
667 ROUND_TO_MINUS = 3 | round result towards minus infinity
669 | Entry points:
671 .globl SYM (__adddf3)
672 .globl SYM (__subdf3)
673 .globl SYM (__muldf3)
674 .globl SYM (__divdf3)
675 .globl SYM (__negdf2)
676 .globl SYM (__cmpdf2)
677 .globl SYM (__cmpdf2_internal)
678 .hidden SYM (__cmpdf2_internal)
680 .text
681 .even
683 | These are common routines to return and signal exceptions.
685 Ld$den:
686 | Return and signal a denormalized number
687 orl d7,d0
688 movew IMM (INEXACT_RESULT+UNDERFLOW),d7
689 moveq IMM (DOUBLE_FLOAT),d6
690 PICJUMP $_exception_handler
692 Ld$infty:
693 Ld$overflow:
694 | Return a properly signed INFINITY and set the exception flags
695 movel IMM (0x7ff00000),d0
696 movel IMM (0),d1
697 orl d7,d0
698 movew IMM (INEXACT_RESULT+OVERFLOW),d7
699 moveq IMM (DOUBLE_FLOAT),d6
700 PICJUMP $_exception_handler
702 Ld$underflow:
703 | Return 0 and set the exception flags
704 movel IMM (0),d0
705 movel d0,d1
706 movew IMM (INEXACT_RESULT+UNDERFLOW),d7
707 moveq IMM (DOUBLE_FLOAT),d6
708 PICJUMP $_exception_handler
710 Ld$inop:
711 | Return a quiet NaN and set the exception flags
712 movel IMM (QUIET_NaN),d0
713 movel d0,d1
714 movew IMM (INEXACT_RESULT+INVALID_OPERATION),d7
715 moveq IMM (DOUBLE_FLOAT),d6
716 PICJUMP $_exception_handler
718 Ld$div$0:
719 | Return a properly signed INFINITY and set the exception flags
720 movel IMM (0x7ff00000),d0
721 movel IMM (0),d1
722 orl d7,d0
723 movew IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
724 moveq IMM (DOUBLE_FLOAT),d6
725 PICJUMP $_exception_handler
727 |=============================================================================
728 |=============================================================================
729 | double precision routines
730 |=============================================================================
731 |=============================================================================
733 | A double precision floating point number (double) has the format:
734 |
735 | struct _double {
736 | unsigned int sign : 1; /* sign bit */
737 | unsigned int exponent : 11; /* exponent, shifted by 126 */
738 | unsigned int fraction : 52; /* fraction */
739 | } double;
740 |
741 | Thus sizeof(double) = 8 (64 bits).
742 |
743 | All the routines are callable from C programs, and return the result
744 | in the register pair d0-d1. They also preserve all registers except
745 | d0-d1 and a0-a1.
747 |=============================================================================
748 | __subdf3
749 |=============================================================================
751 | double __subdf3(double, double);
752 FUNC(__subdf3)
753 SYM (__subdf3):
754 bchg IMM (31),sp@(12) | change sign of second operand
755 | and fall through, so we always add
756 |=============================================================================
757 | __adddf3
758 |=============================================================================
760 | double __adddf3(double, double);
761 FUNC(__adddf3)
762 SYM (__adddf3):
763 #ifndef __mcoldfire__
764 link a6,IMM (0) | everything will be done in registers
765 moveml d2-d7,sp@- | save all data registers and a2 (but d0-d1)
766 #else
767 link a6,IMM (-24)
768 moveml d2-d7,sp@
769 #endif
770 movel a6@(8),d0 | get first operand
771 movel a6@(12),d1 |
772 movel a6@(16),d2 | get second operand
773 movel a6@(20),d3 |
775 movel d0,d7 | get d0's sign bit in d7 '
776 addl d1,d1 | check and clear sign bit of a, and gain one
777 addxl d0,d0 | bit of extra precision
778 beq Ladddf$b | if zero return second operand
780 movel d2,d6 | save sign in d6
781 addl d3,d3 | get rid of sign bit and gain one bit of
782 addxl d2,d2 | extra precision
783 beq Ladddf$a | if zero return first operand
785 andl IMM (0x80000000),d7 | isolate a's sign bit '
786 swap d6 | and also b's sign bit '
787 #ifndef __mcoldfire__
788 andw IMM (0x8000),d6 |
789 orw d6,d7 | and combine them into d7, so that a's sign '
790 | bit is in the high word and b's is in the '
791 | low word, so d6 is free to be used
792 #else
793 andl IMM (0x8000),d6
794 orl d6,d7
795 #endif
796 movel d7,a0 | now save d7 into a0, so d7 is free to
797 | be used also
799 | Get the exponents and check for denormalized and/or infinity.
801 movel IMM (0x001fffff),d6 | mask for the fraction
802 movel IMM (0x00200000),d7 | mask to put hidden bit back
804 movel d0,d4 |
805 andl d6,d0 | get fraction in d0
806 notl d6 | make d6 into mask for the exponent
807 andl d6,d4 | get exponent in d4
808 beq Ladddf$a$den | branch if a is denormalized
809 cmpl d6,d4 | check for INFINITY or NaN
810 beq Ladddf$nf |
811 orl d7,d0 | and put hidden bit back
812 Ladddf$1:
813 swap d4 | shift right exponent so that it starts
814 #ifndef __mcoldfire__
815 lsrw IMM (5),d4 | in bit 0 and not bit 20
816 #else
817 lsrl IMM (5),d4 | in bit 0 and not bit 20
818 #endif
819 | Now we have a's exponent in d4 and fraction in d0-d1 '
820 movel d2,d5 | save b to get exponent
821 andl d6,d5 | get exponent in d5
822 beq Ladddf$b$den | branch if b is denormalized
823 cmpl d6,d5 | check for INFINITY or NaN
824 beq Ladddf$nf
825 notl d6 | make d6 into mask for the fraction again
826 andl d6,d2 | and get fraction in d2
827 orl d7,d2 | and put hidden bit back
828 Ladddf$2:
829 swap d5 | shift right exponent so that it starts
830 #ifndef __mcoldfire__
831 lsrw IMM (5),d5 | in bit 0 and not bit 20
832 #else
833 lsrl IMM (5),d5 | in bit 0 and not bit 20
834 #endif
836 | Now we have b's exponent in d5 and fraction in d2-d3. '
838 | The situation now is as follows: the signs are combined in a0, the
839 | numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
840 | and d5 (b). To do the rounding correctly we need to keep all the
841 | bits until the end, so we need to use d0-d1-d2-d3 for the first number
842 | and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
843 | exponents in a2-a3.
845 #ifndef __mcoldfire__
846 moveml a2-a3,sp@- | save the address registers
847 #else
848 movel a2,sp@-
849 movel a3,sp@-
850 movel a4,sp@-
851 #endif
853 movel d4,a2 | save the exponents
854 movel d5,a3 |
856 movel IMM (0),d7 | and move the numbers around
857 movel d7,d6 |
858 movel d3,d5 |
859 movel d2,d4 |
860 movel d7,d3 |
861 movel d7,d2 |
863 | Here we shift the numbers until the exponents are the same, and put
864 | the largest exponent in a2.
865 #ifndef __mcoldfire__
866 exg d4,a2 | get exponents back
867 exg d5,a3 |
868 cmpw d4,d5 | compare the exponents
869 #else
870 movel d4,a4 | get exponents back
871 movel a2,d4
872 movel a4,a2
873 movel d5,a4
874 movel a3,d5
875 movel a4,a3
876 cmpl d4,d5 | compare the exponents
877 #endif
878 beq Ladddf$3 | if equal don't shift '
879 bhi 9f | branch if second exponent is higher
881 | Here we have a's exponent larger than b's, so we have to shift b. We do
882 | this by using as counter d2:
883 1: movew d4,d2 | move largest exponent to d2
884 #ifndef __mcoldfire__
885 subw d5,d2 | and subtract second exponent
886 exg d4,a2 | get back the longs we saved
887 exg d5,a3 |
888 #else
889 subl d5,d2 | and subtract second exponent
890 movel d4,a4 | get back the longs we saved
891 movel a2,d4
892 movel a4,a2
893 movel d5,a4
894 movel a3,d5
895 movel a4,a3
896 #endif
897 | if difference is too large we don't shift (actually, we can just exit) '
898 #ifndef __mcoldfire__
899 cmpw IMM (DBL_MANT_DIG+2),d2
900 #else
901 cmpl IMM (DBL_MANT_DIG+2),d2
902 #endif
903 bge Ladddf$b$small
904 #ifndef __mcoldfire__
905 cmpw IMM (32),d2 | if difference >= 32, shift by longs
906 #else
907 cmpl IMM (32),d2 | if difference >= 32, shift by longs
908 #endif
909 bge 5f
910 2:
911 #ifndef __mcoldfire__
912 cmpw IMM (16),d2 | if difference >= 16, shift by words
913 #else
914 cmpl IMM (16),d2 | if difference >= 16, shift by words
915 #endif
916 bge 6f
917 bra 3f | enter dbra loop
919 4:
920 #ifndef __mcoldfire__
921 lsrl IMM (1),d4
922 roxrl IMM (1),d5
923 roxrl IMM (1),d6
924 roxrl IMM (1),d7
925 #else
926 lsrl IMM (1),d7
927 btst IMM (0),d6
928 beq 10f
929 bset IMM (31),d7
930 10: lsrl IMM (1),d6
931 btst IMM (0),d5
932 beq 11f
933 bset IMM (31),d6
934 11: lsrl IMM (1),d5
935 btst IMM (0),d4
936 beq 12f
937 bset IMM (31),d5
938 12: lsrl IMM (1),d4
939 #endif
940 3:
941 #ifndef __mcoldfire__
942 dbra d2,4b
943 #else
944 subql IMM (1),d2
945 bpl 4b
946 #endif
947 movel IMM (0),d2
948 movel d2,d3
949 bra Ladddf$4
950 5:
951 movel d6,d7
952 movel d5,d6
953 movel d4,d5
954 movel IMM (0),d4
955 #ifndef __mcoldfire__
956 subw IMM (32),d2
957 #else
958 subl IMM (32),d2
959 #endif
960 bra 2b
961 6:
962 movew d6,d7
963 swap d7
964 movew d5,d6
965 swap d6
966 movew d4,d5
967 swap d5
968 movew IMM (0),d4
969 swap d4
970 #ifndef __mcoldfire__
971 subw IMM (16),d2
972 #else
973 subl IMM (16),d2
974 #endif
975 bra 3b
977 9:
978 #ifndef __mcoldfire__
979 exg d4,d5
980 movew d4,d6
981 subw d5,d6 | keep d5 (largest exponent) in d4
982 exg d4,a2
983 exg d5,a3
984 #else
985 movel d5,d6
986 movel d4,d5
987 movel d6,d4
988 subl d5,d6
989 movel d4,a4
990 movel a2,d4
991 movel a4,a2
992 movel d5,a4
993 movel a3,d5
994 movel a4,a3
995 #endif
996 | if difference is too large we don't shift (actually, we can just exit) '
997 #ifndef __mcoldfire__
998 cmpw IMM (DBL_MANT_DIG+2),d6
999 #else
1000 cmpl IMM (DBL_MANT_DIG+2),d6
1001 #endif
1002 bge Ladddf$a$small
1003 #ifndef __mcoldfire__
1004 cmpw IMM (32),d6 | if difference >= 32, shift by longs
1005 #else
1006 cmpl IMM (32),d6 | if difference >= 32, shift by longs
1007 #endif
1008 bge 5f
1009 2:
1010 #ifndef __mcoldfire__
1011 cmpw IMM (16),d6 | if difference >= 16, shift by words
1012 #else
1013 cmpl IMM (16),d6 | if difference >= 16, shift by words
1014 #endif
1015 bge 6f
1016 bra 3f | enter dbra loop
1018 4:
1019 #ifndef __mcoldfire__
1020 lsrl IMM (1),d0
1021 roxrl IMM (1),d1
1022 roxrl IMM (1),d2
1023 roxrl IMM (1),d3
1024 #else
1025 lsrl IMM (1),d3
1026 btst IMM (0),d2
1027 beq 10f
1028 bset IMM (31),d3
1029 10: lsrl IMM (1),d2
1030 btst IMM (0),d1
1031 beq 11f
1032 bset IMM (31),d2
1033 11: lsrl IMM (1),d1
1034 btst IMM (0),d0
1035 beq 12f
1036 bset IMM (31),d1
1037 12: lsrl IMM (1),d0
1038 #endif
1039 3:
1040 #ifndef __mcoldfire__
1041 dbra d6,4b
1042 #else
1043 subql IMM (1),d6
1044 bpl 4b
1045 #endif
1046 movel IMM (0),d7
1047 movel d7,d6
1048 bra Ladddf$4
1049 5:
1050 movel d2,d3
1051 movel d1,d2
1052 movel d0,d1
1053 movel IMM (0),d0
1054 #ifndef __mcoldfire__
1055 subw IMM (32),d6
1056 #else
1057 subl IMM (32),d6
1058 #endif
1059 bra 2b
1060 6:
1061 movew d2,d3
1062 swap d3
1063 movew d1,d2
1064 swap d2
1065 movew d0,d1
1066 swap d1
1067 movew IMM (0),d0
1068 swap d0
1069 #ifndef __mcoldfire__
1070 subw IMM (16),d6
1071 #else
1072 subl IMM (16),d6
1073 #endif
1074 bra 3b
1075 Ladddf$3:
1076 #ifndef __mcoldfire__
1077 exg d4,a2
1078 exg d5,a3
1079 #else
1080 movel d4,a4
1081 movel a2,d4
1082 movel a4,a2
1083 movel d5,a4
1084 movel a3,d5
1085 movel a4,a3
1086 #endif
1087 Ladddf$4:
1088 | Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
1089 | the signs in a4.
1091 | Here we have to decide whether to add or subtract the numbers:
1092 #ifndef __mcoldfire__
1093 exg d7,a0 | get the signs
1094 exg d6,a3 | a3 is free to be used
1095 #else
1096 movel d7,a4
1097 movel a0,d7
1098 movel a4,a0
1099 movel d6,a4
1100 movel a3,d6
1101 movel a4,a3
1102 #endif
1103 movel d7,d6 |
1104 movew IMM (0),d7 | get a's sign in d7 '
1105 swap d6 |
1106 movew IMM (0),d6 | and b's sign in d6 '
1107 eorl d7,d6 | compare the signs
1108 bmi Lsubdf$0 | if the signs are different we have
1109 | to subtract
1110 #ifndef __mcoldfire__
1111 exg d7,a0 | else we add the numbers
1112 exg d6,a3 |
1113 #else
1114 movel d7,a4
1115 movel a0,d7
1116 movel a4,a0
1117 movel d6,a4
1118 movel a3,d6
1119 movel a4,a3
1120 #endif
1121 addl d7,d3 |
1122 addxl d6,d2 |
1123 addxl d5,d1 |
1124 addxl d4,d0 |
1126 movel a2,d4 | return exponent to d4
1127 movel a0,d7 |
1128 andl IMM (0x80000000),d7 | d7 now has the sign
1130 #ifndef __mcoldfire__
1131 moveml sp@+,a2-a3
1132 #else
1133 movel sp@+,a4
1134 movel sp@+,a3
1135 movel sp@+,a2
1136 #endif
1138 | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
1139 | the case of denormalized numbers in the rounding routine itself).
1140 | As in the addition (not in the subtraction!) we could have set
1141 | one more bit we check this:
1142 btst IMM (DBL_MANT_DIG+1),d0
1143 beq 1f
1144 #ifndef __mcoldfire__
1145 lsrl IMM (1),d0
1146 roxrl IMM (1),d1
1147 roxrl IMM (1),d2
1148 roxrl IMM (1),d3
1149 addw IMM (1),d4
1150 #else
1151 lsrl IMM (1),d3
1152 btst IMM (0),d2
1153 beq 10f
1154 bset IMM (31),d3
1155 10: lsrl IMM (1),d2
1156 btst IMM (0),d1
1157 beq 11f
1158 bset IMM (31),d2
1159 11: lsrl IMM (1),d1
1160 btst IMM (0),d0
1161 beq 12f
1162 bset IMM (31),d1
1163 12: lsrl IMM (1),d0
1164 addl IMM (1),d4
1165 #endif
1166 1:
1167 lea pc@(Ladddf$5),a0 | to return from rounding routine
1168 PICLEA SYM (_fpCCR),a1 | check the rounding mode
1169 #ifdef __mcoldfire__
1170 clrl d6
1171 #endif
1172 movew a1@(6),d6 | rounding mode in d6
1173 beq Lround$to$nearest
1174 #ifndef __mcoldfire__
1175 cmpw IMM (ROUND_TO_PLUS),d6
1176 #else
1177 cmpl IMM (ROUND_TO_PLUS),d6
1178 #endif
1179 bhi Lround$to$minus
1180 blt Lround$to$zero
1181 bra Lround$to$plus
1182 Ladddf$5:
1183 | Put back the exponent and check for overflow
1184 #ifndef __mcoldfire__
1185 cmpw IMM (0x7ff),d4 | is the exponent big?
1186 #else
1187 cmpl IMM (0x7ff),d4 | is the exponent big?
1188 #endif
1189 bge 1f
1190 bclr IMM (DBL_MANT_DIG-1),d0
1191 #ifndef __mcoldfire__
1192 lslw IMM (4),d4 | put exponent back into position
1193 #else
1194 lsll IMM (4),d4 | put exponent back into position
1195 #endif
1196 swap d0 |
1197 #ifndef __mcoldfire__
1198 orw d4,d0 |
1199 #else
1200 orl d4,d0 |
1201 #endif
1202 swap d0 |
1203 bra Ladddf$ret
1204 1:
1205 moveq IMM (ADD),d5
1206 bra Ld$overflow
1208 Lsubdf$0:
1209 | Here we do the subtraction.
1210 #ifndef __mcoldfire__
1211 exg d7,a0 | put sign back in a0
1212 exg d6,a3 |
1213 #else
1214 movel d7,a4
1215 movel a0,d7
1216 movel a4,a0
1217 movel d6,a4
1218 movel a3,d6
1219 movel a4,a3
1220 #endif
1221 subl d7,d3 |
1222 subxl d6,d2 |
1223 subxl d5,d1 |
1224 subxl d4,d0 |
1225 beq Ladddf$ret$1 | if zero just exit
1226 bpl 1f | if positive skip the following
1227 movel a0,d7 |
1228 bchg IMM (31),d7 | change sign bit in d7
1229 movel d7,a0 |
1230 negl d3 |
1231 negxl d2 |
1232 negxl d1 | and negate result
1233 negxl d0 |
1234 1:
1235 movel a2,d4 | return exponent to d4
1236 movel a0,d7
1237 andl IMM (0x80000000),d7 | isolate sign bit
1238 #ifndef __mcoldfire__
1239 moveml sp@+,a2-a3 |
1240 #else
1241 movel sp@+,a4
1242 movel sp@+,a3
1243 movel sp@+,a2
1244 #endif
1246 | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
1247 | the case of denormalized numbers in the rounding routine itself).
1248 | As in the addition (not in the subtraction!) we could have set
1249 | one more bit we check this:
1250 btst IMM (DBL_MANT_DIG+1),d0
1251 beq 1f
1252 #ifndef __mcoldfire__
1253 lsrl IMM (1),d0
1254 roxrl IMM (1),d1
1255 roxrl IMM (1),d2
1256 roxrl IMM (1),d3
1257 addw IMM (1),d4
1258 #else
1259 lsrl IMM (1),d3
1260 btst IMM (0),d2
1261 beq 10f
1262 bset IMM (31),d3
1263 10: lsrl IMM (1),d2
1264 btst IMM (0),d1
1265 beq 11f
1266 bset IMM (31),d2
1267 11: lsrl IMM (1),d1
1268 btst IMM (0),d0
1269 beq 12f
1270 bset IMM (31),d1
1271 12: lsrl IMM (1),d0
1272 addl IMM (1),d4
1273 #endif
1274 1:
1275 lea pc@(Lsubdf$1),a0 | to return from rounding routine
1276 PICLEA SYM (_fpCCR),a1 | check the rounding mode
1277 #ifdef __mcoldfire__
1278 clrl d6
1279 #endif
1280 movew a1@(6),d6 | rounding mode in d6
1281 beq Lround$to$nearest
1282 #ifndef __mcoldfire__
1283 cmpw IMM (ROUND_TO_PLUS),d6
1284 #else
1285 cmpl IMM (ROUND_TO_PLUS),d6
1286 #endif
1287 bhi Lround$to$minus
1288 blt Lround$to$zero
1289 bra Lround$to$plus
1290 Lsubdf$1:
1291 | Put back the exponent and sign (we don't have overflow). '
1292 bclr IMM (DBL_MANT_DIG-1),d0
1293 #ifndef __mcoldfire__
1294 lslw IMM (4),d4 | put exponent back into position
1295 #else
1296 lsll IMM (4),d4 | put exponent back into position
1297 #endif
1298 swap d0 |
1299 #ifndef __mcoldfire__
1300 orw d4,d0 |
1301 #else
1302 orl d4,d0 |
1303 #endif
1304 swap d0 |
1305 bra Ladddf$ret
1307 | If one of the numbers was too small (difference of exponents >=
1308 | DBL_MANT_DIG+1) we return the other (and now we don't have to '
1309 | check for finiteness or zero).
1310 Ladddf$a$small:
1311 #ifndef __mcoldfire__
1312 moveml sp@+,a2-a3
1313 #else
1314 movel sp@+,a4
1315 movel sp@+,a3
1316 movel sp@+,a2
1317 #endif
1318 movel a6@(16),d0
1319 movel a6@(20),d1
1320 PICLEA SYM (_fpCCR),a0
1321 movew IMM (0),a0@
1322 #ifndef __mcoldfire__
1323 moveml sp@+,d2-d7 | restore data registers
1324 #else
1325 moveml sp@,d2-d7
1326 | XXX if frame pointer is ever removed, stack pointer must
1327 | be adjusted here.
1328 #endif
1329 unlk a6 | and return
1330 rts
1332 Ladddf$b$small:
1333 #ifndef __mcoldfire__
1334 moveml sp@+,a2-a3
1335 #else
1336 movel sp@+,a4
1337 movel sp@+,a3
1338 movel sp@+,a2
1339 #endif
1340 movel a6@(8),d0
1341 movel a6@(12),d1
1342 PICLEA SYM (_fpCCR),a0
1343 movew IMM (0),a0@
1344 #ifndef __mcoldfire__
1345 moveml sp@+,d2-d7 | restore data registers
1346 #else
1347 moveml sp@,d2-d7
1348 | XXX if frame pointer is ever removed, stack pointer must
1349 | be adjusted here.
1350 #endif
1351 unlk a6 | and return
1352 rts
1354 Ladddf$a$den:
1355 movel d7,d4 | d7 contains 0x00200000
1356 bra Ladddf$1
1358 Ladddf$b$den:
1359 movel d7,d5 | d7 contains 0x00200000
1360 notl d6
1361 bra Ladddf$2
1363 Ladddf$b:
1364 | Return b (if a is zero)
1365 movel d2,d0
1366 movel d3,d1
1367 bne 1f | Check if b is -0
1368 cmpl IMM (0x80000000),d0
1369 bne 1f
1370 andl IMM (0x80000000),d7 | Use the sign of a
1371 clrl d0
1372 bra Ladddf$ret
1373 Ladddf$a:
1374 movel a6@(8),d0
1375 movel a6@(12),d1
1376 1:
1377 moveq IMM (ADD),d5
1378 | Check for NaN and +/-INFINITY.
1379 movel d0,d7 |
1380 andl IMM (0x80000000),d7 |
1381 bclr IMM (31),d0 |
1382 cmpl IMM (0x7ff00000),d0 |
1383 bge 2f |
1384 movel d0,d0 | check for zero, since we don't '
1385 bne Ladddf$ret | want to return -0 by mistake
1386 movel d1,d1 |
1387 bne Ladddf$ret |
1388 bclr IMM (31),d7 |
1389 bra Ladddf$ret |
1390 2:
1391 andl IMM (0x000fffff),d0 | check for NaN (nonzero fraction)
1392 orl d1,d0 |
1393 bne Ld$inop |
1394 bra Ld$infty |
1396 Ladddf$ret$1:
1397 #ifndef __mcoldfire__
1398 moveml sp@+,a2-a3 | restore regs and exit
1399 #else
1400 movel sp@+,a4
1401 movel sp@+,a3
1402 movel sp@+,a2
1403 #endif
1405 Ladddf$ret:
1406 | Normal exit.
1407 PICLEA SYM (_fpCCR),a0
1408 movew IMM (0),a0@
1409 orl d7,d0 | put sign bit back
1410 #ifndef __mcoldfire__
1411 moveml sp@+,d2-d7
1412 #else
1413 moveml sp@,d2-d7
1414 | XXX if frame pointer is ever removed, stack pointer must
1415 | be adjusted here.
1416 #endif
1417 unlk a6
1418 rts
1420 Ladddf$ret$den:
1421 | Return a denormalized number.
1422 #ifndef __mcoldfire__
1423 lsrl IMM (1),d0 | shift right once more
1424 roxrl IMM (1),d1 |
1425 #else
1426 lsrl IMM (1),d1
1427 btst IMM (0),d0
1428 beq 10f
1429 bset IMM (31),d1
1430 10: lsrl IMM (1),d0
1431 #endif
1432 bra Ladddf$ret
1434 Ladddf$nf:
1435 moveq IMM (ADD),d5
1436 | This could be faster but it is not worth the effort, since it is not
1437 | executed very often. We sacrifice speed for clarity here.
1438 movel a6@(8),d0 | get the numbers back (remember that we
1439 movel a6@(12),d1 | did some processing already)
1440 movel a6@(16),d2 |
1441 movel a6@(20),d3 |
1442 movel IMM (0x7ff00000),d4 | useful constant (INFINITY)
1443 movel d0,d7 | save sign bits
1444 movel d2,d6 |
1445 bclr IMM (31),d0 | clear sign bits
1446 bclr IMM (31),d2 |
1447 | We know that one of them is either NaN of +/-INFINITY
1448 | Check for NaN (if either one is NaN return NaN)
1449 cmpl d4,d0 | check first a (d0)
1450 bhi Ld$inop | if d0 > 0x7ff00000 or equal and
1451 bne 2f
1452 tstl d1 | d1 > 0, a is NaN
1453 bne Ld$inop |
1454 2: cmpl d4,d2 | check now b (d1)
1455 bhi Ld$inop |
1456 bne 3f
1457 tstl d3 |
1458 bne Ld$inop |
1459 3:
1460 | Now comes the check for +/-INFINITY. We know that both are (maybe not
1461 | finite) numbers, but we have to check if both are infinite whether we
1462 | are adding or subtracting them.
1463 eorl d7,d6 | to check sign bits
1464 bmi 1f
1465 andl IMM (0x80000000),d7 | get (common) sign bit
1466 bra Ld$infty
1467 1:
1468 | We know one (or both) are infinite, so we test for equality between the
1469 | two numbers (if they are equal they have to be infinite both, so we
1470 | return NaN).
1471 cmpl d2,d0 | are both infinite?
1472 bne 1f | if d0 <> d2 they are not equal
1473 cmpl d3,d1 | if d0 == d2 test d3 and d1
1474 beq Ld$inop | if equal return NaN
1475 1:
1476 andl IMM (0x80000000),d7 | get a's sign bit '
1477 cmpl d4,d0 | test now for infinity
1478 beq Ld$infty | if a is INFINITY return with this sign
1479 bchg IMM (31),d7 | else we know b is INFINITY and has
1480 bra Ld$infty | the opposite sign
1482 |=============================================================================
1483 | __muldf3
1484 |=============================================================================
1486 | double __muldf3(double, double);
1487 FUNC(__muldf3)
1488 SYM (__muldf3):
1489 #ifndef __mcoldfire__
1490 link a6,IMM (0)
1491 moveml d2-d7,sp@-
1492 #else
1493 link a6,IMM (-24)
1494 moveml d2-d7,sp@
1495 #endif
1496 movel a6@(8),d0 | get a into d0-d1
1497 movel a6@(12),d1 |
1498 movel a6@(16),d2 | and b into d2-d3
1499 movel a6@(20),d3 |
1500 movel d0,d7 | d7 will hold the sign of the product
1501 eorl d2,d7 |
1502 andl IMM (0x80000000),d7 |
1503 movel d7,a0 | save sign bit into a0
1504 movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
1505 movel d7,d6 | another (mask for fraction)
1506 notl d6 |
1507 bclr IMM (31),d0 | get rid of a's sign bit '
1508 movel d0,d4 |
1509 orl d1,d4 |
1510 beq Lmuldf$a$0 | branch if a is zero
1511 movel d0,d4 |
1512 bclr IMM (31),d2 | get rid of b's sign bit '
1513 movel d2,d5 |
1514 orl d3,d5 |
1515 beq Lmuldf$b$0 | branch if b is zero
1516 movel d2,d5 |
1517 cmpl d7,d0 | is a big?
1518 bhi Lmuldf$inop | if a is NaN return NaN
1519 beq Lmuldf$a$nf | we still have to check d1 and b ...
1520 cmpl d7,d2 | now compare b with INFINITY
1521 bhi Lmuldf$inop | is b NaN?
1522 beq Lmuldf$b$nf | we still have to check d3 ...
1523 | Here we have both numbers finite and nonzero (and with no sign bit).
1524 | Now we get the exponents into d4 and d5.
1525 andl d7,d4 | isolate exponent in d4
1526 beq Lmuldf$a$den | if exponent zero, have denormalized
1527 andl d6,d0 | isolate fraction
1528 orl IMM (0x00100000),d0 | and put hidden bit back
1529 swap d4 | I like exponents in the first byte
1530 #ifndef __mcoldfire__
1531 lsrw IMM (4),d4 |
1532 #else
1533 lsrl IMM (4),d4 |
1534 #endif
1535 Lmuldf$1:
1536 andl d7,d5 |
1537 beq Lmuldf$b$den |
1538 andl d6,d2 |
1539 orl IMM (0x00100000),d2 | and put hidden bit back
1540 swap d5 |
1541 #ifndef __mcoldfire__
1542 lsrw IMM (4),d5 |
1543 #else
1544 lsrl IMM (4),d5 |
1545 #endif
1546 Lmuldf$2: |
1547 #ifndef __mcoldfire__
1548 addw d5,d4 | add exponents
1549 subw IMM (D_BIAS+1),d4 | and subtract bias (plus one)
1550 #else
1551 addl d5,d4 | add exponents
1552 subl IMM (D_BIAS+1),d4 | and subtract bias (plus one)
1553 #endif
1555 | We are now ready to do the multiplication. The situation is as follows:
1556 | both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were
1557 | denormalized to start with!), which means that in the product bit 104
1558 | (which will correspond to bit 8 of the fourth long) is set.
1560 | Here we have to do the product.
1561 | To do it we have to juggle the registers back and forth, as there are not
1562 | enough to keep everything in them. So we use the address registers to keep
1563 | some intermediate data.
1565 #ifndef __mcoldfire__
1566 moveml a2-a3,sp@- | save a2 and a3 for temporary use
1567 #else
1568 movel a2,sp@-
1569 movel a3,sp@-
1570 movel a4,sp@-
1571 #endif
1572 movel IMM (0),a2 | a2 is a null register
1573 movel d4,a3 | and a3 will preserve the exponent
1575 | First, shift d2-d3 so bit 20 becomes bit 31:
1576 #ifndef __mcoldfire__
1577 rorl IMM (5),d2 | rotate d2 5 places right
1578 swap d2 | and swap it
1579 rorl IMM (5),d3 | do the same thing with d3
1580 swap d3 |
1581 movew d3,d6 | get the rightmost 11 bits of d3
1582 andw IMM (0x07ff),d6 |
1583 orw d6,d2 | and put them into d2
1584 andw IMM (0xf800),d3 | clear those bits in d3
1585 #else
1586 moveq IMM (11),d7 | left shift d2 11 bits
1587 lsll d7,d2
1588 movel d3,d6 | get a copy of d3
1589 lsll d7,d3 | left shift d3 11 bits
1590 andl IMM (0xffe00000),d6 | get the top 11 bits of d3
1591 moveq IMM (21),d7 | right shift them 21 bits
1592 lsrl d7,d6
1593 orl d6,d2 | stick them at the end of d2
1594 #endif
1596 movel d2,d6 | move b into d6-d7
1597 movel d3,d7 | move a into d4-d5
1598 movel d0,d4 | and clear d0-d1-d2-d3 (to put result)
1599 movel d1,d5 |
1600 movel IMM (0),d3 |
1601 movel d3,d2 |
1602 movel d3,d1 |
1603 movel d3,d0 |
1605 | We use a1 as counter:
1606 movel IMM (DBL_MANT_DIG-1),a1
1607 #ifndef __mcoldfire__
1608 exg d7,a1
1609 #else
1610 movel d7,a4
1611 movel a1,d7
1612 movel a4,a1
1613 #endif
1615 1:
1616 #ifndef __mcoldfire__
1617 exg d7,a1 | put counter back in a1
1618 #else
1619 movel d7,a4
1620 movel a1,d7
1621 movel a4,a1
1622 #endif
1623 addl d3,d3 | shift sum once left
1624 addxl d2,d2 |
1625 addxl d1,d1 |
1626 addxl d0,d0 |
1627 addl d7,d7 |
1628 addxl d6,d6 |
1629 bcc 2f | if bit clear skip the following
1630 #ifndef __mcoldfire__
1631 exg d7,a2 |
1632 #else
1633 movel d7,a4
1634 movel a2,d7
1635 movel a4,a2
1636 #endif
1637 addl d5,d3 | else add a to the sum
1638 addxl d4,d2 |
1639 addxl d7,d1 |
1640 addxl d7,d0 |
1641 #ifndef __mcoldfire__
1642 exg d7,a2 |
1643 #else
1644 movel d7,a4
1645 movel a2,d7
1646 movel a4,a2
1647 #endif
1648 2:
1649 #ifndef __mcoldfire__
1650 exg d7,a1 | put counter in d7
1651 dbf d7,1b | decrement and branch
1652 #else
1653 movel d7,a4
1654 movel a1,d7
1655 movel a4,a1
1656 subql IMM (1),d7
1657 bpl 1b
1658 #endif
1660 movel a3,d4 | restore exponent
1661 #ifndef __mcoldfire__
1662 moveml sp@+,a2-a3
1663 #else
1664 movel sp@+,a4
1665 movel sp@+,a3
1666 movel sp@+,a2
1667 #endif
1669 | Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The
1670 | first thing to do now is to normalize it so bit 8 becomes bit
1671 | DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
1672 swap d0
1673 swap d1
1674 movew d1,d0
1675 swap d2
1676 movew d2,d1
1677 swap d3
1678 movew d3,d2
1679 movew IMM (0),d3
1680 #ifndef __mcoldfire__
1681 lsrl IMM (1),d0
1682 roxrl IMM (1),d1
1683 roxrl IMM (1),d2
1684 roxrl IMM (1),d3
1685 lsrl IMM (1),d0
1686 roxrl IMM (1),d1
1687 roxrl IMM (1),d2
1688 roxrl IMM (1),d3
1689 lsrl IMM (1),d0
1690 roxrl IMM (1),d1
1691 roxrl IMM (1),d2
1692 roxrl IMM (1),d3
1693 #else
1694 moveq IMM (29),d6
1695 lsrl IMM (3),d3
1696 movel d2,d7
1697 lsll d6,d7
1698 orl d7,d3
1699 lsrl IMM (3),d2
1700 movel d1,d7
1701 lsll d6,d7
1702 orl d7,d2
1703 lsrl IMM (3),d1
1704 movel d0,d7
1705 lsll d6,d7
1706 orl d7,d1
1707 lsrl IMM (3),d0
1708 #endif
1710 | Now round, check for over- and underflow, and exit.
1711 movel a0,d7 | get sign bit back into d7
1712 moveq IMM (MULTIPLY),d5
1714 btst IMM (DBL_MANT_DIG+1-32),d0
1715 beq Lround$exit
1716 #ifndef __mcoldfire__
1717 lsrl IMM (1),d0
1718 roxrl IMM (1),d1
1719 addw IMM (1),d4
1720 #else
1721 lsrl IMM (1),d1
1722 btst IMM (0),d0
1723 beq 10f
1724 bset IMM (31),d1
1725 10: lsrl IMM (1),d0
1726 addl IMM (1),d4
1727 #endif
1728 bra Lround$exit
1730 Lmuldf$inop:
1731 moveq IMM (MULTIPLY),d5
1732 bra Ld$inop
1734 Lmuldf$b$nf:
1735 moveq IMM (MULTIPLY),d5
1736 movel a0,d7 | get sign bit back into d7
1737 tstl d3 | we know d2 == 0x7ff00000, so check d3
1738 bne Ld$inop | if d3 <> 0 b is NaN
1739 bra Ld$overflow | else we have overflow (since a is finite)
1741 Lmuldf$a$nf:
1742 moveq IMM (MULTIPLY),d5
1743 movel a0,d7 | get sign bit back into d7
1744 tstl d1 | we know d0 == 0x7ff00000, so check d1
1745 bne Ld$inop | if d1 <> 0 a is NaN
1746 bra Ld$overflow | else signal overflow
1748 | If either number is zero return zero, unless the other is +/-INFINITY or
1749 | NaN, in which case we return NaN.
1750 Lmuldf$b$0:
1751 moveq IMM (MULTIPLY),d5
1752 #ifndef __mcoldfire__
1753 exg d2,d0 | put b (==0) into d0-d1
1754 exg d3,d1 | and a (with sign bit cleared) into d2-d3
1755 movel a0,d0 | set result sign
1756 #else
1757 movel d0,d2 | put a into d2-d3
1758 movel d1,d3
1759 movel a0,d0 | put result zero into d0-d1
1760 movq IMM(0),d1
1761 #endif
1762 bra 1f
1763 Lmuldf$a$0:
1764 movel a0,d0 | set result sign
1765 movel a6@(16),d2 | put b into d2-d3 again
1766 movel a6@(20),d3 |
1767 bclr IMM (31),d2 | clear sign bit
1768 1: cmpl IMM (0x7ff00000),d2 | check for non-finiteness
1769 bge Ld$inop | in case NaN or +/-INFINITY return NaN
1770 PICLEA SYM (_fpCCR),a0
1771 movew IMM (0),a0@
1772 #ifndef __mcoldfire__
1773 moveml sp@+,d2-d7
1774 #else
1775 moveml sp@,d2-d7
1776 | XXX if frame pointer is ever removed, stack pointer must
1777 | be adjusted here.
1778 #endif
1779 unlk a6
1780 rts
1782 | If a number is denormalized we put an exponent of 1 but do not put the
1783 | hidden bit back into the fraction; instead we shift left until bit 21
1784 | (the hidden bit) is set, adjusting the exponent accordingly. We do this
1785 | to ensure that the product of the fractions is close to 1.
1786 Lmuldf$a$den:
1787 movel IMM (1),d4
1788 andl d6,d0
1789 1: addl d1,d1 | shift a left until bit 20 is set
1790 addxl d0,d0 |
1791 #ifndef __mcoldfire__
1792 subw IMM (1),d4 | and adjust exponent
1793 #else
1794 subl IMM (1),d4 | and adjust exponent
1795 #endif
1796 btst IMM (20),d0 |
1797 bne Lmuldf$1 |
1798 bra 1b
1800 Lmuldf$b$den:
1801 movel IMM (1),d5
1802 andl d6,d2
1803 1: addl d3,d3 | shift b left until bit 20 is set
1804 addxl d2,d2 |
1805 #ifndef __mcoldfire__
1806 subw IMM (1),d5 | and adjust exponent
1807 #else
1808 subql IMM (1),d5 | and adjust exponent
1809 #endif
1810 btst IMM (20),d2 |
1811 bne Lmuldf$2 |
1812 bra 1b
1815 |=============================================================================
1816 | __divdf3
1817 |=============================================================================
1819 | double __divdf3(double, double);
1820 FUNC(__divdf3)
1821 SYM (__divdf3):
1822 #ifndef __mcoldfire__
1823 link a6,IMM (0)
1824 moveml d2-d7,sp@-
1825 #else
1826 link a6,IMM (-24)
1827 moveml d2-d7,sp@
1828 #endif
1829 movel a6@(8),d0 | get a into d0-d1
1830 movel a6@(12),d1 |
1831 movel a6@(16),d2 | and b into d2-d3
1832 movel a6@(20),d3 |
1833 movel d0,d7 | d7 will hold the sign of the result
1834 eorl d2,d7 |
1835 andl IMM (0x80000000),d7
1836 movel d7,a0 | save sign into a0
1837 movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
1838 movel d7,d6 | another (mask for fraction)
1839 notl d6 |
1840 bclr IMM (31),d0 | get rid of a's sign bit '
1841 movel d0,d4 |
1842 orl d1,d4 |
1843 beq Ldivdf$a$0 | branch if a is zero
1844 movel d0,d4 |
1845 bclr IMM (31),d2 | get rid of b's sign bit '
1846 movel d2,d5 |
1847 orl d3,d5 |
1848 beq Ldivdf$b$0 | branch if b is zero
1849 movel d2,d5
1850 cmpl d7,d0 | is a big?
1851 bhi Ldivdf$inop | if a is NaN return NaN
1852 beq Ldivdf$a$nf | if d0 == 0x7ff00000 we check d1
1853 cmpl d7,d2 | now compare b with INFINITY
1854 bhi Ldivdf$inop | if b is NaN return NaN
1855 beq Ldivdf$b$nf | if d2 == 0x7ff00000 we check d3
1856 | Here we have both numbers finite and nonzero (and with no sign bit).
1857 | Now we get the exponents into d4 and d5 and normalize the numbers to
1858 | ensure that the ratio of the fractions is around 1. We do this by
1859 | making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
1860 | set, even if they were denormalized to start with.
1861 | Thus, the result will satisfy: 2 > result > 1/2.
1862 andl d7,d4 | and isolate exponent in d4
1863 beq Ldivdf$a$den | if exponent is zero we have a denormalized
1864 andl d6,d0 | and isolate fraction
1865 orl IMM (0x00100000),d0 | and put hidden bit back
1866 swap d4 | I like exponents in the first byte
1867 #ifndef __mcoldfire__
1868 lsrw IMM (4),d4 |
1869 #else
1870 lsrl IMM (4),d4 |
1871 #endif
1872 Ldivdf$1: |
1873 andl d7,d5 |
1874 beq Ldivdf$b$den |
1875 andl d6,d2 |
1876 orl IMM (0x00100000),d2
1877 swap d5 |
1878 #ifndef __mcoldfire__
1879 lsrw IMM (4),d5 |
1880 #else
1881 lsrl IMM (4),d5 |
1882 #endif
1883 Ldivdf$2: |
1884 #ifndef __mcoldfire__
1885 subw d5,d4 | subtract exponents
1886 addw IMM (D_BIAS),d4 | and add bias
1887 #else
1888 subl d5,d4 | subtract exponents
1889 addl IMM (D_BIAS),d4 | and add bias
1890 #endif
1892 | We are now ready to do the division. We have prepared things in such a way
1893 | that the ratio of the fractions will be less than 2 but greater than 1/2.
1894 | At this point the registers in use are:
1895 | d0-d1 hold a (first operand, bit DBL_MANT_DIG-32=0, bit
1896 | DBL_MANT_DIG-1-32=1)
1897 | d2-d3 hold b (second operand, bit DBL_MANT_DIG-32=1)
1898 | d4 holds the difference of the exponents, corrected by the bias
1899 | a0 holds the sign of the ratio
1901 | To do the rounding correctly we need to keep information about the
1902 | nonsignificant bits. One way to do this would be to do the division
1903 | using four registers; another is to use two registers (as originally
1904 | I did), but use a sticky bit to preserve information about the
1905 | fractional part. Note that we can keep that info in a1, which is not
1906 | used.
1907 movel IMM (0),d6 | d6-d7 will hold the result
1908 movel d6,d7 |
1909 movel IMM (0),a1 | and a1 will hold the sticky bit
1911 movel IMM (DBL_MANT_DIG-32+1),d5
1913 1: cmpl d0,d2 | is a < b?
1914 bhi 3f | if b > a skip the following
1915 beq 4f | if d0==d2 check d1 and d3
1916 2: subl d3,d1 |
1917 subxl d2,d0 | a <-- a - b
1918 bset d5,d6 | set the corresponding bit in d6
1919 3: addl d1,d1 | shift a by 1
1920 addxl d0,d0 |
1921 #ifndef __mcoldfire__
1922 dbra d5,1b | and branch back
1923 #else
1924 subql IMM (1), d5
1925 bpl 1b
1926 #endif
1927 bra 5f
1928 4: cmpl d1,d3 | here d0==d2, so check d1 and d3
1929 bhi 3b | if d1 > d2 skip the subtraction
1930 bra 2b | else go do it
1931 5:
1932 | Here we have to start setting the bits in the second long.
1933 movel IMM (31),d5 | again d5 is counter
1935 1: cmpl d0,d2 | is a < b?
1936 bhi 3f | if b > a skip the following
1937 beq 4f | if d0==d2 check d1 and d3
1938 2: subl d3,d1 |
1939 subxl d2,d0 | a <-- a - b
1940 bset d5,d7 | set the corresponding bit in d7
1941 3: addl d1,d1 | shift a by 1
1942 addxl d0,d0 |
1943 #ifndef __mcoldfire__
1944 dbra d5,1b | and branch back
1945 #else
1946 subql IMM (1), d5
1947 bpl 1b
1948 #endif
1949 bra 5f
1950 4: cmpl d1,d3 | here d0==d2, so check d1 and d3
1951 bhi 3b | if d1 > d2 skip the subtraction
1952 bra 2b | else go do it
1953 5:
1954 | Now go ahead checking until we hit a one, which we store in d2.
1955 movel IMM (DBL_MANT_DIG),d5
1956 1: cmpl d2,d0 | is a < b?
1957 bhi 4f | if b < a, exit
1958 beq 3f | if d0==d2 check d1 and d3
1959 2: addl d1,d1 | shift a by 1
1960 addxl d0,d0 |
1961 #ifndef __mcoldfire__
1962 dbra d5,1b | and branch back
1963 #else
1964 subql IMM (1), d5
1965 bpl 1b
1966 #endif
1967 movel IMM (0),d2 | here no sticky bit was found
1968 movel d2,d3
1969 bra 5f
1970 3: cmpl d1,d3 | here d0==d2, so check d1 and d3
1971 bhi 2b | if d1 > d2 go back
1972 4:
1973 | Here put the sticky bit in d2-d3 (in the position which actually corresponds
1974 | to it; if you don't do this the algorithm loses in some cases). '
1975 movel IMM (0),d2
1976 movel d2,d3
1977 #ifndef __mcoldfire__
1978 subw IMM (DBL_MANT_DIG),d5
1979 addw IMM (63),d5
1980 cmpw IMM (31),d5
1981 #else
1982 subl IMM (DBL_MANT_DIG),d5
1983 addl IMM (63),d5
1984 cmpl IMM (31),d5
1985 #endif
1986 bhi 2f
1987 1: bset d5,d3
1988 bra 5f
1989 #ifndef __mcoldfire__
1990 subw IMM (32),d5
1991 #else
1992 subl IMM (32),d5
1993 #endif
1994 2: bset d5,d2
1995 5:
1996 | Finally we are finished! Move the longs in the address registers to
1997 | their final destination:
1998 movel d6,d0
1999 movel d7,d1
2000 movel IMM (0),d3
2002 | Here we have finished the division, with the result in d0-d1-d2-d3, with
2003 | 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
2004 | If it is not, then definitely bit 21 is set. Normalize so bit 22 is
2005 | not set:
2006 btst IMM (DBL_MANT_DIG-32+1),d0
2007 beq 1f
2008 #ifndef __mcoldfire__
2009 lsrl IMM (1),d0
2010 roxrl IMM (1),d1
2011 roxrl IMM (1),d2
2012 roxrl IMM (1),d3
2013 addw IMM (1),d4
2014 #else
2015 lsrl IMM (1),d3
2016 btst IMM (0),d2
2017 beq 10f
2018 bset IMM (31),d3
2019 10: lsrl IMM (1),d2
2020 btst IMM (0),d1
2021 beq 11f
2022 bset IMM (31),d2
2023 11: lsrl IMM (1),d1
2024 btst IMM (0),d0
2025 beq 12f
2026 bset IMM (31),d1
2027 12: lsrl IMM (1),d0
2028 addl IMM (1),d4
2029 #endif
2030 1:
2031 | Now round, check for over- and underflow, and exit.
2032 movel a0,d7 | restore sign bit to d7
2033 moveq IMM (DIVIDE),d5
2034 bra Lround$exit
2036 Ldivdf$inop:
2037 moveq IMM (DIVIDE),d5
2038 bra Ld$inop
2040 Ldivdf$a$0:
2041 | If a is zero check to see whether b is zero also. In that case return
2042 | NaN; then check if b is NaN, and return NaN also in that case. Else
2043 | return a properly signed zero.
2044 moveq IMM (DIVIDE),d5
2045 bclr IMM (31),d2 |
2046 movel d2,d4 |
2047 orl d3,d4 |
2048 beq Ld$inop | if b is also zero return NaN
2049 cmpl IMM (0x7ff00000),d2 | check for NaN
2050 bhi Ld$inop |
2051 blt 1f |
2052 tstl d3 |
2053 bne Ld$inop |
2054 1: movel a0,d0 | else return signed zero
2055 moveq IMM(0),d1 |
2056 PICLEA SYM (_fpCCR),a0 | clear exception flags
2057 movew IMM (0),a0@ |
2058 #ifndef __mcoldfire__
2059 moveml sp@+,d2-d7 |
2060 #else
2061 moveml sp@,d2-d7 |
2062 | XXX if frame pointer is ever removed, stack pointer must
2063 | be adjusted here.
2064 #endif
2065 unlk a6 |
2066 rts |
2068 Ldivdf$b$0:
2069 moveq IMM (DIVIDE),d5
2070 | If we got here a is not zero. Check if a is NaN; in that case return NaN,
2071 | else return +/-INFINITY. Remember that a is in d0 with the sign bit
2072 | cleared already.
2073 movel a0,d7 | put a's sign bit back in d7 '
2074 cmpl IMM (0x7ff00000),d0 | compare d0 with INFINITY
2075 bhi Ld$inop | if larger it is NaN
2076 tstl d1 |
2077 bne Ld$inop |
2078 bra Ld$div$0 | else signal DIVIDE_BY_ZERO
2080 Ldivdf$b$nf:
2081 moveq IMM (DIVIDE),d5
2082 | If d2 == 0x7ff00000 we have to check d3.
2083 tstl d3 |
2084 bne Ld$inop | if d3 <> 0, b is NaN
2085 bra Ld$underflow | else b is +/-INFINITY, so signal underflow
2087 Ldivdf$a$nf:
2088 moveq IMM (DIVIDE),d5
2089 | If d0 == 0x7ff00000 we have to check d1.
2090 tstl d1 |
2091 bne Ld$inop | if d1 <> 0, a is NaN
2092 | If a is INFINITY we have to check b
2093 cmpl d7,d2 | compare b with INFINITY
2094 bge Ld$inop | if b is NaN or INFINITY return NaN
2095 movl a0,d7 | restore sign bit to d7
2096 bra Ld$overflow | else return overflow
2098 | If a number is denormalized we put an exponent of 1 but do not put the
2099 | bit back into the fraction.
2100 Ldivdf$a$den:
2101 movel IMM (1),d4
2102 andl d6,d0
2103 1: addl d1,d1 | shift a left until bit 20 is set
2104 addxl d0,d0
2105 #ifndef __mcoldfire__
2106 subw IMM (1),d4 | and adjust exponent
2107 #else
2108 subl IMM (1),d4 | and adjust exponent
2109 #endif
2110 btst IMM (DBL_MANT_DIG-32-1),d0
2111 bne Ldivdf$1
2112 bra 1b
2114 Ldivdf$b$den:
2115 movel IMM (1),d5
2116 andl d6,d2
2117 1: addl d3,d3 | shift b left until bit 20 is set
2118 addxl d2,d2
2119 #ifndef __mcoldfire__
2120 subw IMM (1),d5 | and adjust exponent
2121 #else
2122 subql IMM (1),d5 | and adjust exponent
2123 #endif
2124 btst IMM (DBL_MANT_DIG-32-1),d2
2125 bne Ldivdf$2
2126 bra 1b
2128 Lround$exit:
2129 | This is a common exit point for __muldf3 and __divdf3. When they enter
2130 | this point the sign of the result is in d7, the result in d0-d1, normalized
2131 | so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.
2133 | First check for underlow in the exponent:
2134 #ifndef __mcoldfire__
2135 cmpw IMM (-DBL_MANT_DIG-1),d4
2136 #else
2137 cmpl IMM (-DBL_MANT_DIG-1),d4
2138 #endif
2139 blt Ld$underflow
2140 | It could happen that the exponent is less than 1, in which case the
2141 | number is denormalized. In this case we shift right and adjust the
2142 | exponent until it becomes 1 or the fraction is zero (in the latter case
2143 | we signal underflow and return zero).
2144 movel d7,a0 |
2145 movel IMM (0),d6 | use d6-d7 to collect bits flushed right
2146 movel d6,d7 | use d6-d7 to collect bits flushed right
2147 #ifndef __mcoldfire__
2148 cmpw IMM (1),d4 | if the exponent is less than 1 we
2149 #else
2150 cmpl IMM (1),d4 | if the exponent is less than 1 we
2151 #endif
2152 bge 2f | have to shift right (denormalize)
2153 1:
2154 #ifndef __mcoldfire__
2155 addw IMM (1),d4 | adjust the exponent
2156 lsrl IMM (1),d0 | shift right once
2157 roxrl IMM (1),d1 |
2158 roxrl IMM (1),d2 |
2159 roxrl IMM (1),d3 |
2160 roxrl IMM (1),d6 |
2161 roxrl IMM (1),d7 |
2162 cmpw IMM (1),d4 | is the exponent 1 already?
2163 #else
2164 addl IMM (1),d4 | adjust the exponent
2165 lsrl IMM (1),d7
2166 btst IMM (0),d6
2167 beq 13f
2168 bset IMM (31),d7
2169 13: lsrl IMM (1),d6
2170 btst IMM (0),d3
2171 beq 14f
2172 bset IMM (31),d6
2173 14: lsrl IMM (1),d3
2174 btst IMM (0),d2
2175 beq 10f
2176 bset IMM (31),d3
2177 10: lsrl IMM (1),d2
2178 btst IMM (0),d1
2179 beq 11f
2180 bset IMM (31),d2
2181 11: lsrl IMM (1),d1
2182 btst IMM (0),d0
2183 beq 12f
2184 bset IMM (31),d1
2185 12: lsrl IMM (1),d0
2186 cmpl IMM (1),d4 | is the exponent 1 already?
2187 #endif
2188 beq 2f | if not loop back
2189 bra 1b |
2190 bra Ld$underflow | safety check, shouldn't execute '
2191 2: orl d6,d2 | this is a trick so we don't lose '
2192 orl d7,d3 | the bits which were flushed right
2193 movel a0,d7 | get back sign bit into d7
2194 | Now call the rounding routine (which takes care of denormalized numbers):
2195 lea pc@(Lround$0),a0 | to return from rounding routine
2196 PICLEA SYM (_fpCCR),a1 | check the rounding mode
2197 #ifdef __mcoldfire__
2198 clrl d6
2199 #endif
2200 movew a1@(6),d6 | rounding mode in d6
2201 beq Lround$to$nearest
2202 #ifndef __mcoldfire__
2203 cmpw IMM (ROUND_TO_PLUS),d6
2204 #else
2205 cmpl IMM (ROUND_TO_PLUS),d6
2206 #endif
2207 bhi Lround$to$minus
2208 blt Lround$to$zero
2209 bra Lround$to$plus
2210 Lround$0:
2211 | Here we have a correctly rounded result (either normalized or denormalized).
2213 | Here we should have either a normalized number or a denormalized one, and
2214 | the exponent is necessarily larger or equal to 1 (so we don't have to '
2215 | check again for underflow!). We have to check for overflow or for a
2216 | denormalized number (which also signals underflow).
2217 | Check for overflow (i.e., exponent >= 0x7ff).
2218 #ifndef __mcoldfire__
2219 cmpw IMM (0x07ff),d4
2220 #else
2221 cmpl IMM (0x07ff),d4
2222 #endif
2223 bge Ld$overflow
2224 | Now check for a denormalized number (exponent==0):
2225 movew d4,d4
2226 beq Ld$den
2227 1:
2228 | Put back the exponents and sign and return.
2229 #ifndef __mcoldfire__
2230 lslw IMM (4),d4 | exponent back to fourth byte
2231 #else
2232 lsll IMM (4),d4 | exponent back to fourth byte
2233 #endif
2234 bclr IMM (DBL_MANT_DIG-32-1),d0
2235 swap d0 | and put back exponent
2236 #ifndef __mcoldfire__
2237 orw d4,d0 |
2238 #else
2239 orl d4,d0 |
2240 #endif
2241 swap d0 |
2242 orl d7,d0 | and sign also
2244 PICLEA SYM (_fpCCR),a0
2245 movew IMM (0),a0@
2246 #ifndef __mcoldfire__
2247 moveml sp@+,d2-d7
2248 #else
2249 moveml sp@,d2-d7
2250 | XXX if frame pointer is ever removed, stack pointer must
2251 | be adjusted here.
2252 #endif
2253 unlk a6
2254 rts
2256 |=============================================================================
2257 | __negdf2
2258 |=============================================================================
2260 | double __negdf2(double, double);
2261 FUNC(__negdf2)
2262 SYM (__negdf2):
2263 #ifndef __mcoldfire__
2264 link a6,IMM (0)
2265 moveml d2-d7,sp@-
2266 #else
2267 link a6,IMM (-24)
2268 moveml d2-d7,sp@
2269 #endif
2270 moveq IMM (NEGATE),d5
2271 movel a6@(8),d0 | get number to negate in d0-d1
2272 movel a6@(12),d1 |
2273 bchg IMM (31),d0 | negate
2274 movel d0,d2 | make a positive copy (for the tests)
2275 bclr IMM (31),d2 |
2276 movel d2,d4 | check for zero
2277 orl d1,d4 |
2278 beq 2f | if zero (either sign) return +zero
2279 cmpl IMM (0x7ff00000),d2 | compare to +INFINITY
2280 blt 1f | if finite, return
2281 bhi Ld$inop | if larger (fraction not zero) is NaN
2282 tstl d1 | if d2 == 0x7ff00000 check d1
2283 bne Ld$inop |
2284 movel d0,d7 | else get sign and return INFINITY
2285 andl IMM (0x80000000),d7
2286 bra Ld$infty
2287 1: PICLEA SYM (_fpCCR),a0
2288 movew IMM (0),a0@
2289 #ifndef __mcoldfire__
2290 moveml sp@+,d2-d7
2291 #else
2292 moveml sp@,d2-d7
2293 | XXX if frame pointer is ever removed, stack pointer must
2294 | be adjusted here.
2295 #endif
2296 unlk a6
2297 rts
2298 2: bclr IMM (31),d0
2299 bra 1b
2301 |=============================================================================
2302 | __cmpdf2
2303 |=============================================================================
2305 GREATER = 1
2306 LESS = -1
2307 EQUAL = 0
2309 | int __cmpdf2_internal(double, double, int);
2310 SYM (__cmpdf2_internal):
2311 #ifndef __mcoldfire__
2312 link a6,IMM (0)
2313 moveml d2-d7,sp@- | save registers
2314 #else
2315 link a6,IMM (-24)
2316 moveml d2-d7,sp@
2317 #endif
2318 moveq IMM (COMPARE),d5
2319 movel a6@(8),d0 | get first operand
2320 movel a6@(12),d1 |
2321 movel a6@(16),d2 | get second operand
2322 movel a6@(20),d3 |
2323 | First check if a and/or b are (+/-) zero and in that case clear
2324 | the sign bit.
2325 movel d0,d6 | copy signs into d6 (a) and d7(b)
2326 bclr IMM (31),d0 | and clear signs in d0 and d2
2327 movel d2,d7 |
2328 bclr IMM (31),d2 |
2329 cmpl IMM (0x7ff00000),d0 | check for a == NaN
2330 bhi Lcmpd$inop | if d0 > 0x7ff00000, a is NaN
2331 beq Lcmpdf$a$nf | if equal can be INFINITY, so check d1
2332 movel d0,d4 | copy into d4 to test for zero
2333 orl d1,d4 |
2334 beq Lcmpdf$a$0 |
2335 Lcmpdf$0:
2336 cmpl IMM (0x7ff00000),d2 | check for b == NaN
2337 bhi Lcmpd$inop | if d2 > 0x7ff00000, b is NaN
2338 beq Lcmpdf$b$nf | if equal can be INFINITY, so check d3
2339 movel d2,d4 |
2340 orl d3,d4 |
2341 beq Lcmpdf$b$0 |
2342 Lcmpdf$1:
2343 | Check the signs
2344 eorl d6,d7
2345 bpl 1f
2346 | If the signs are not equal check if a >= 0
2347 tstl d6
2348 bpl Lcmpdf$a$gt$b | if (a >= 0 && b < 0) => a > b
2349 bmi Lcmpdf$b$gt$a | if (a < 0 && b >= 0) => a < b
2350 1:
2351 | If the signs are equal check for < 0
2352 tstl d6
2353 bpl 1f
2354 | If both are negative exchange them
2355 #ifndef __mcoldfire__
2356 exg d0,d2
2357 exg d1,d3
2358 #else
2359 movel d0,d7
2360 movel d2,d0
2361 movel d7,d2
2362 movel d1,d7
2363 movel d3,d1
2364 movel d7,d3
2365 #endif
2366 1:
2367 | Now that they are positive we just compare them as longs (does this also
2368 | work for denormalized numbers?).
2369 cmpl d0,d2
2370 bhi Lcmpdf$b$gt$a | |b| > |a|
2371 bne Lcmpdf$a$gt$b | |b| < |a|
2372 | If we got here d0 == d2, so we compare d1 and d3.
2373 cmpl d1,d3
2374 bhi Lcmpdf$b$gt$a | |b| > |a|
2375 bne Lcmpdf$a$gt$b | |b| < |a|
2376 | If we got here a == b.
2377 movel IMM (EQUAL),d0
2378 #ifndef __mcoldfire__
2379 moveml sp@+,d2-d7 | put back the registers
2380 #else
2381 moveml sp@,d2-d7
2382 | XXX if frame pointer is ever removed, stack pointer must
2383 | be adjusted here.
2384 #endif
2385 unlk a6
2386 rts
2387 Lcmpdf$a$gt$b:
2388 movel IMM (GREATER),d0
2389 #ifndef __mcoldfire__
2390 moveml sp@+,d2-d7 | put back the registers
2391 #else
2392 moveml sp@,d2-d7
2393 | XXX if frame pointer is ever removed, stack pointer must
2394 | be adjusted here.
2395 #endif
2396 unlk a6
2397 rts
2398 Lcmpdf$b$gt$a:
2399 movel IMM (LESS),d0
2400 #ifndef __mcoldfire__
2401 moveml sp@+,d2-d7 | put back the registers
2402 #else
2403 moveml sp@,d2-d7
2404 | XXX if frame pointer is ever removed, stack pointer must
2405 | be adjusted here.
2406 #endif
2407 unlk a6
2408 rts
2410 Lcmpdf$a$0:
2411 bclr IMM (31),d6
2412 bra Lcmpdf$0
2413 Lcmpdf$b$0:
2414 bclr IMM (31),d7
2415 bra Lcmpdf$1
2417 Lcmpdf$a$nf:
2418 tstl d1
2419 bne Ld$inop
2420 bra Lcmpdf$0
2422 Lcmpdf$b$nf:
2423 tstl d3
2424 bne Ld$inop
2425 bra Lcmpdf$1
2427 Lcmpd$inop:
2428 movl a6@(24),d0
2429 moveq IMM (INEXACT_RESULT+INVALID_OPERATION),d7
2430 moveq IMM (DOUBLE_FLOAT),d6
2431 PICJUMP $_exception_handler
2433 | int __cmpdf2(double, double);
2434 FUNC(__cmpdf2)
2435 SYM (__cmpdf2):
2436 link a6,IMM (0)
2437 pea 1
2438 movl a6@(20),sp@-
2439 movl a6@(16),sp@-
2440 movl a6@(12),sp@-
2441 movl a6@(8),sp@-
2442 PICCALL SYM (__cmpdf2_internal)
2443 unlk a6
2444 rts
2446 |=============================================================================
2447 | rounding routines
2448 |=============================================================================
2450 | The rounding routines expect the number to be normalized in registers
2451 | d0-d1-d2-d3, with the exponent in register d4. They assume that the
2452 | exponent is larger or equal to 1. They return a properly normalized number
2453 | if possible, and a denormalized number otherwise. The exponent is returned
2454 | in d4.
2456 Lround$to$nearest:
2457 | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
2458 | Here we assume that the exponent is not too small (this should be checked
2459 | before entering the rounding routine), but the number could be denormalized.
2461 | Check for denormalized numbers:
2462 1: btst IMM (DBL_MANT_DIG-32),d0
2463 bne 2f | if set the number is normalized
2464 | Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent
2465 | is one (remember that a denormalized number corresponds to an
2466 | exponent of -D_BIAS+1).
2467 #ifndef __mcoldfire__
2468 cmpw IMM (1),d4 | remember that the exponent is at least one
2469 #else
2470 cmpl IMM (1),d4 | remember that the exponent is at least one
2471 #endif
2472 beq 2f | an exponent of one means denormalized
2473 addl d3,d3 | else shift and adjust the exponent
2474 addxl d2,d2 |
2475 addxl d1,d1 |
2476 addxl d0,d0 |
2477 #ifndef __mcoldfire__
2478 dbra d4,1b |
2479 #else
2480 subql IMM (1), d4
2481 bpl 1b
2482 #endif
2483 2:
2484 | Now round: we do it as follows: after the shifting we can write the
2485 | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
2486 | If delta < 1, do nothing. If delta > 1, add 1 to f.
2487 | If delta == 1, we make sure the rounded number will be even (odd?)
2488 | (after shifting).
2489 btst IMM (0),d1 | is delta < 1?
2490 beq 2f | if so, do not do anything
2491 orl d2,d3 | is delta == 1?
2492 bne 1f | if so round to even
2493 movel d1,d3 |
2494 andl IMM (2),d3 | bit 1 is the last significant bit
2495 movel IMM (0),d2 |
2496 addl d3,d1 |
2497 addxl d2,d0 |
2498 bra 2f |
2499 1: movel IMM (1),d3 | else add 1
2500 movel IMM (0),d2 |
2501 addl d3,d1 |
2502 addxl d2,d0
2503 | Shift right once (because we used bit #DBL_MANT_DIG-32!).
2504 2:
2505 #ifndef __mcoldfire__
2506 lsrl IMM (1),d0
2507 roxrl IMM (1),d1
2508 #else
2509 lsrl IMM (1),d1
2510 btst IMM (0),d0
2511 beq 10f
2512 bset IMM (31),d1
2513 10: lsrl IMM (1),d0
2514 #endif
2516 | Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
2517 | 'fraction overflow' ...).
2518 btst IMM (DBL_MANT_DIG-32),d0
2519 beq 1f
2520 #ifndef __mcoldfire__
2521 lsrl IMM (1),d0
2522 roxrl IMM (1),d1
2523 addw IMM (1),d4
2524 #else
2525 lsrl IMM (1),d1
2526 btst IMM (0),d0
2527 beq 10f
2528 bset IMM (31),d1
2529 10: lsrl IMM (1),d0
2530 addl IMM (1),d4
2531 #endif
2532 1:
2533 | If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we
2534 | have to put the exponent to zero and return a denormalized number.
2535 btst IMM (DBL_MANT_DIG-32-1),d0
2536 beq 1f
2537 jmp a0@
2538 1: movel IMM (0),d4
2539 jmp a0@
2541 Lround$to$zero:
2542 Lround$to$plus:
2543 Lround$to$minus:
2544 jmp a0@
2545 #endif /* L_double */
2547 #ifdef L_float
2549 .globl SYM (_fpCCR)
2550 .globl $_exception_handler
2552 QUIET_NaN = 0xffffffff
2553 SIGNL_NaN = 0x7f800001
2554 INFINITY = 0x7f800000
2556 F_MAX_EXP = 0xff
2557 F_BIAS = 126
2558 FLT_MAX_EXP = F_MAX_EXP - F_BIAS
2559 FLT_MIN_EXP = 1 - F_BIAS
2560 FLT_MANT_DIG = 24
2562 INEXACT_RESULT = 0x0001
2563 UNDERFLOW = 0x0002
2564 OVERFLOW = 0x0004
2565 DIVIDE_BY_ZERO = 0x0008
2566 INVALID_OPERATION = 0x0010
2568 SINGLE_FLOAT = 1
2570 NOOP = 0
2571 ADD = 1
2572 MULTIPLY = 2
2573 DIVIDE = 3
2574 NEGATE = 4
2575 COMPARE = 5
2576 EXTENDSFDF = 6
2577 TRUNCDFSF = 7
2579 UNKNOWN = -1
2580 ROUND_TO_NEAREST = 0 | round result to nearest representable value
2581 ROUND_TO_ZERO = 1 | round result towards zero
2582 ROUND_TO_PLUS = 2 | round result towards plus infinity
2583 ROUND_TO_MINUS = 3 | round result towards minus infinity
2585 | Entry points:
2587 .globl SYM (__addsf3)
2588 .globl SYM (__subsf3)
2589 .globl SYM (__mulsf3)
2590 .globl SYM (__divsf3)
2591 .globl SYM (__negsf2)
2592 .globl SYM (__cmpsf2)
2593 .globl SYM (__cmpsf2_internal)
2594 .hidden SYM (__cmpsf2_internal)
2596 | These are common routines to return and signal exceptions.
2598 .text
2599 .even
2601 Lf$den:
2602 | Return and signal a denormalized number
2603 orl d7,d0
2604 moveq IMM (INEXACT_RESULT+UNDERFLOW),d7
2605 moveq IMM (SINGLE_FLOAT),d6
2606 PICJUMP $_exception_handler
2608 Lf$infty:
2609 Lf$overflow:
2610 | Return a properly signed INFINITY and set the exception flags
2611 movel IMM (INFINITY),d0
2612 orl d7,d0
2613 moveq IMM (INEXACT_RESULT+OVERFLOW),d7
2614 moveq IMM (SINGLE_FLOAT),d6
2615 PICJUMP $_exception_handler
2617 Lf$underflow:
2618 | Return 0 and set the exception flags
2619 moveq IMM (0),d0
2620 moveq IMM (INEXACT_RESULT+UNDERFLOW),d7
2621 moveq IMM (SINGLE_FLOAT),d6
2622 PICJUMP $_exception_handler
2624 Lf$inop:
2625 | Return a quiet NaN and set the exception flags
2626 movel IMM (QUIET_NaN),d0
2627 moveq IMM (INEXACT_RESULT+INVALID_OPERATION),d7
2628 moveq IMM (SINGLE_FLOAT),d6
2629 PICJUMP $_exception_handler
2631 Lf$div$0:
2632 | Return a properly signed INFINITY and set the exception flags
2633 movel IMM (INFINITY),d0
2634 orl d7,d0
2635 moveq IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
2636 moveq IMM (SINGLE_FLOAT),d6
2637 PICJUMP $_exception_handler
2639 |=============================================================================
2640 |=============================================================================
2641 | single precision routines
2642 |=============================================================================
2643 |=============================================================================
2645 | A single precision floating point number (float) has the format:
2646 |
2647 | struct _float {
2648 | unsigned int sign : 1; /* sign bit */
2649 | unsigned int exponent : 8; /* exponent, shifted by 126 */
2650 | unsigned int fraction : 23; /* fraction */
2651 | } float;
2652 |
2653 | Thus sizeof(float) = 4 (32 bits).
2654 |
2655 | All the routines are callable from C programs, and return the result
2656 | in the single register d0. They also preserve all registers except
2657 | d0-d1 and a0-a1.
2659 |=============================================================================
2660 | __subsf3
2661 |=============================================================================
2663 | float __subsf3(float, float);
2664 FUNC(__subsf3)
2665 SYM (__subsf3):
2666 bchg IMM (31),sp@(8) | change sign of second operand
2667 | and fall through
2668 |=============================================================================
2669 | __addsf3
2670 |=============================================================================
2672 | float __addsf3(float, float);
2673 FUNC(__addsf3)
2674 SYM (__addsf3):
2675 #ifndef __mcoldfire__
2676 link a6,IMM (0) | everything will be done in registers
2677 moveml d2-d7,sp@- | save all data registers but d0-d1
2678 #else
2679 link a6,IMM (-24)
2680 moveml d2-d7,sp@
2681 #endif
2682 movel a6@(8),d0 | get first operand
2683 movel a6@(12),d1 | get second operand
2684 movel d0,a0 | get d0's sign bit '
2685 addl d0,d0 | check and clear sign bit of a
2686 beq Laddsf$b | if zero return second operand
2687 movel d1,a1 | save b's sign bit '
2688 addl d1,d1 | get rid of sign bit
2689 beq Laddsf$a | if zero return first operand
2691 | Get the exponents and check for denormalized and/or infinity.
2693 movel IMM (0x00ffffff),d4 | mask to get fraction
2694 movel IMM (0x01000000),d5 | mask to put hidden bit back
2696 movel d0,d6 | save a to get exponent
2697 andl d4,d0 | get fraction in d0
2698 notl d4 | make d4 into a mask for the exponent
2699 andl d4,d6 | get exponent in d6
2700 beq Laddsf$a$den | branch if a is denormalized
2701 cmpl d4,d6 | check for INFINITY or NaN
2702 beq Laddsf$nf
2703 swap d6 | put exponent into first word
2704 orl d5,d0 | and put hidden bit back
2705 Laddsf$1:
2706 | Now we have a's exponent in d6 (second byte) and the mantissa in d0. '
2707 movel d1,d7 | get exponent in d7
2708 andl d4,d7 |
2709 beq Laddsf$b$den | branch if b is denormalized
2710 cmpl d4,d7 | check for INFINITY or NaN
2711 beq Laddsf$nf
2712 swap d7 | put exponent into first word
2713 notl d4 | make d4 into a mask for the fraction
2714 andl d4,d1 | get fraction in d1
2715 orl d5,d1 | and put hidden bit back
2716 Laddsf$2:
2717 | Now we have b's exponent in d7 (second byte) and the mantissa in d1. '
2719 | Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we
2720 | shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
2721 | bit).
2723 movel d1,d2 | move b to d2, since we want to use
2724 | two registers to do the sum
2725 movel IMM (0),d1 | and clear the new ones
2726 movel d1,d3 |
2728 | Here we shift the numbers in registers d0 and d1 so the exponents are the
2729 | same, and put the largest exponent in d6. Note that we are using two
2730 | registers for each number (see the discussion by D. Knuth in "Seminumerical
2731 | Algorithms").
2732 #ifndef __mcoldfire__
2733 cmpw d6,d7 | compare exponents
2734 #else
2735 cmpl d6,d7 | compare exponents
2736 #endif
2737 beq Laddsf$3 | if equal don't shift '
2738 bhi 5f | branch if second exponent largest
2739 1:
2740 subl d6,d7 | keep the largest exponent
2741 negl d7
2742 #ifndef __mcoldfire__
2743 lsrw IMM (8),d7 | put difference in lower byte
2744 #else
2745 lsrl IMM (8),d7 | put difference in lower byte
2746 #endif
2747 | if difference is too large we don't shift (actually, we can just exit) '
2748 #ifndef __mcoldfire__
2749 cmpw IMM (FLT_MANT_DIG+2),d7
2750 #else
2751 cmpl IMM (FLT_MANT_DIG+2),d7
2752 #endif
2753 bge Laddsf$b$small
2754 #ifndef __mcoldfire__
2755 cmpw IMM (16),d7 | if difference >= 16 swap
2756 #else
2757 cmpl IMM (16),d7 | if difference >= 16 swap
2758 #endif
2759 bge 4f
2760 2:
2761 #ifndef __mcoldfire__
2762 subw IMM (1),d7
2763 #else
2764 subql IMM (1), d7
2765 #endif
2766 3:
2767 #ifndef __mcoldfire__
2768 lsrl IMM (1),d2 | shift right second operand
2769 roxrl IMM (1),d3
2770 dbra d7,3b
2771 #else
2772 lsrl IMM (1),d3
2773 btst IMM (0),d2
2774 beq 10f
2775 bset IMM (31),d3
2776 10: lsrl IMM (1),d2
2777 subql IMM (1), d7
2778 bpl 3b
2779 #endif
2780 bra Laddsf$3
2781 4:
2782 movew d2,d3
2783 swap d3
2784 movew d3,d2
2785 swap d2
2786 #ifndef __mcoldfire__
2787 subw IMM (16),d7
2788 #else
2789 subl IMM (16),d7
2790 #endif
2791 bne 2b | if still more bits, go back to normal case
2792 bra Laddsf$3
2793 5:
2794 #ifndef __mcoldfire__
2795 exg d6,d7 | exchange the exponents
2796 #else
2797 eorl d6,d7
2798 eorl d7,d6
2799 eorl d6,d7
2800 #endif
2801 subl d6,d7 | keep the largest exponent
2802 negl d7 |
2803 #ifndef __mcoldfire__
2804 lsrw IMM (8),d7 | put difference in lower byte
2805 #else
2806 lsrl IMM (8),d7 | put difference in lower byte
2807 #endif
2808 | if difference is too large we don't shift (and exit!) '
2809 #ifndef __mcoldfire__
2810 cmpw IMM (FLT_MANT_DIG+2),d7
2811 #else
2812 cmpl IMM (FLT_MANT_DIG+2),d7
2813 #endif
2814 bge Laddsf$a$small
2815 #ifndef __mcoldfire__
2816 cmpw IMM (16),d7 | if difference >= 16 swap
2817 #else
2818 cmpl IMM (16),d7 | if difference >= 16 swap
2819 #endif
2820 bge 8f
2821 6:
2822 #ifndef __mcoldfire__
2823 subw IMM (1),d7
2824 #else
2825 subl IMM (1),d7
2826 #endif
2827 7:
2828 #ifndef __mcoldfire__
2829 lsrl IMM (1),d0 | shift right first operand
2830 roxrl IMM (1),d1
2831 dbra d7,7b
2832 #else
2833 lsrl IMM (1),d1
2834 btst IMM (0),d0
2835 beq 10f
2836 bset IMM (31),d1
2837 10: lsrl IMM (1),d0
2838 subql IMM (1),d7
2839 bpl 7b
2840 #endif
2841 bra Laddsf$3
2842 8:
2843 movew d0,d1
2844 swap d1
2845 movew d1,d0
2846 swap d0
2847 #ifndef __mcoldfire__
2848 subw IMM (16),d7
2849 #else
2850 subl IMM (16),d7
2851 #endif
2852 bne 6b | if still more bits, go back to normal case
2853 | otherwise we fall through
2855 | Now we have a in d0-d1, b in d2-d3, and the largest exponent in d6 (the
2856 | signs are stored in a0 and a1).
2858 Laddsf$3:
2859 | Here we have to decide whether to add or subtract the numbers
2860 #ifndef __mcoldfire__
2861 exg d6,a0 | get signs back
2862 exg d7,a1 | and save the exponents
2863 #else
2864 movel d6,d4
2865 movel a0,d6
2866 movel d4,a0
2867 movel d7,d4
2868 movel a1,d7
2869 movel d4,a1
2870 #endif
2871 eorl d6,d7 | combine sign bits
2872 bmi Lsubsf$0 | if negative a and b have opposite
2873 | sign so we actually subtract the
2874 | numbers
2876 | Here we have both positive or both negative
2877 #ifndef __mcoldfire__
2878 exg d6,a0 | now we have the exponent in d6
2879 #else
2880 movel d6,d4
2881 movel a0,d6
2882 movel d4,a0
2883 #endif
2884 movel a0,d7 | and sign in d7
2885 andl IMM (0x80000000),d7
2886 | Here we do the addition.
2887 addl d3,d1
2888 addxl d2,d0
2889 | Note: now we have d2, d3, d4 and d5 to play with!
2891 | Put the exponent, in the first byte, in d2, to use the "standard" rounding
2892 | routines:
2893 movel d6,d2
2894 #ifndef __mcoldfire__
2895 lsrw IMM (8),d2
2896 #else
2897 lsrl IMM (8),d2
2898 #endif
2900 | Before rounding normalize so bit #FLT_MANT_DIG is set (we will consider
2901 | the case of denormalized numbers in the rounding routine itself).
2902 | As in the addition (not in the subtraction!) we could have set
2903 | one more bit we check this:
2904 btst IMM (FLT_MANT_DIG+1),d0
2905 beq 1f
2906 #ifndef __mcoldfire__
2907 lsrl IMM (1),d0
2908 roxrl IMM (1),d1
2909 #else
2910 lsrl IMM (1),d1
2911 btst IMM (0),d0
2912 beq 10f
2913 bset IMM (31),d1
2914 10: lsrl IMM (1),d0
2915 #endif
2916 addl IMM (1),d2
2917 1:
2918 lea pc@(Laddsf$4),a0 | to return from rounding routine
2919 PICLEA SYM (_fpCCR),a1 | check the rounding mode
2920 #ifdef __mcoldfire__
2921 clrl d6
2922 #endif
2923 movew a1@(6),d6 | rounding mode in d6
2924 beq Lround$to$nearest
2925 #ifndef __mcoldfire__
2926 cmpw IMM (ROUND_TO_PLUS),d6
2927 #else
2928 cmpl IMM (ROUND_TO_PLUS),d6
2929 #endif
2930 bhi Lround$to$minus
2931 blt Lround$to$zero
2932 bra Lround$to$plus
2933 Laddsf$4:
2934 | Put back the exponent, but check for overflow.
2935 #ifndef __mcoldfire__
2936 cmpw IMM (0xff),d2
2937 #else
2938 cmpl IMM (0xff),d2
2939 #endif
2940 bge 1f
2941 bclr IMM (FLT_MANT_DIG-1),d0
2942 #ifndef __mcoldfire__
2943 lslw IMM (7),d2
2944 #else
2945 lsll IMM (7),d2
2946 #endif
2947 swap d2
2948 orl d2,d0
2949 bra Laddsf$ret
2950 1:
2951 moveq IMM (ADD),d5
2952 bra Lf$overflow
2954 Lsubsf$0:
2955 | We are here if a > 0 and b < 0 (sign bits cleared).
2956 | Here we do the subtraction.
2957 movel d6,d7 | put sign in d7
2958 andl IMM (0x80000000),d7
2960 subl d3,d1 | result in d0-d1
2961 subxl d2,d0 |
2962 beq Laddsf$ret | if zero just exit
2963 bpl 1f | if positive skip the following
2964 bchg IMM (31),d7 | change sign bit in d7
2965 negl d1
2966 negxl d0
2967 1:
2968 #ifndef __mcoldfire__
2969 exg d2,a0 | now we have the exponent in d2
2970 lsrw IMM (8),d2 | put it in the first byte
2971 #else
2972 movel d2,d4
2973 movel a0,d2
2974 movel d4,a0
2975 lsrl IMM (8),d2 | put it in the first byte
2976 #endif
2978 | Now d0-d1 is positive and the sign bit is in d7.
2980 | Note that we do not have to normalize, since in the subtraction bit
2981 | #FLT_MANT_DIG+1 is never set, and denormalized numbers are handled by
2982 | the rounding routines themselves.
2983 lea pc@(Lsubsf$1),a0 | to return from rounding routine
2984 PICLEA SYM (_fpCCR),a1 | check the rounding mode
2985 #ifdef __mcoldfire__
2986 clrl d6
2987 #endif
2988 movew a1@(6),d6 | rounding mode in d6
2989 beq Lround$to$nearest
2990 #ifndef __mcoldfire__
2991 cmpw IMM (ROUND_TO_PLUS),d6
2992 #else
2993 cmpl IMM (ROUND_TO_PLUS),d6
2994 #endif
2995 bhi Lround$to$minus
2996 blt Lround$to$zero
2997 bra Lround$to$plus
2998 Lsubsf$1:
2999 | Put back the exponent (we can't have overflow!). '
3000 bclr IMM (FLT_MANT_DIG-1),d0
3001 #ifndef __mcoldfire__
3002 lslw IMM (7),d2
3003 #else
3004 lsll IMM (7),d2
3005 #endif
3006 swap d2
3007 orl d2,d0
3008 bra Laddsf$ret
3010 | If one of the numbers was too small (difference of exponents >=
3011 | FLT_MANT_DIG+2) we return the other (and now we don't have to '
3012 | check for finiteness or zero).
3013 Laddsf$a$small:
3014 movel a6@(12),d0
3015 PICLEA SYM (_fpCCR),a0
3016 movew IMM (0),a0@
3017 #ifndef __mcoldfire__
3018 moveml sp@+,d2-d7 | restore data registers
3019 #else
3020 moveml sp@,d2-d7
3021 | XXX if frame pointer is ever removed, stack pointer must
3022 | be adjusted here.
3023 #endif
3024 unlk a6 | and return
3025 rts
3027 Laddsf$b$small:
3028 movel a6@(8),d0
3029 PICLEA SYM (_fpCCR),a0
3030 movew IMM (0),a0@
3031 #ifndef __mcoldfire__
3032 moveml sp@+,d2-d7 | restore data registers
3033 #else
3034 moveml sp@,d2-d7
3035 | XXX if frame pointer is ever removed, stack pointer must
3036 | be adjusted here.
3037 #endif
3038 unlk a6 | and return
3039 rts
3041 | If the numbers are denormalized remember to put exponent equal to 1.
3043 Laddsf$a$den:
3044 movel d5,d6 | d5 contains 0x01000000
3045 swap d6
3046 bra Laddsf$1
3048 Laddsf$b$den:
3049 movel d5,d7
3050 swap d7
3051 notl d4 | make d4 into a mask for the fraction
3052 | (this was not executed after the jump)
3053 bra Laddsf$2
3055 | The rest is mainly code for the different results which can be
3056 | returned (checking always for +/-INFINITY and NaN).
3058 Laddsf$b:
3059 | Return b (if a is zero).
3060 movel a6@(12),d0
3061 cmpl IMM (0x80000000),d0 | Check if b is -0
3062 bne 1f
3063 movel a0,d7
3064 andl IMM (0x80000000),d7 | Use the sign of a
3065 clrl d0
3066 bra Laddsf$ret
3067 Laddsf$a:
3068 | Return a (if b is zero).
3069 movel a6@(8),d0
3070 1:
3071 moveq IMM (ADD),d5
3072 | We have to check for NaN and +/-infty.
3073 movel d0,d7
3074 andl IMM (0x80000000),d7 | put sign in d7
3075 bclr IMM (31),d0 | clear sign
3076 cmpl IMM (INFINITY),d0 | check for infty or NaN
3077 bge 2f
3078 movel d0,d0 | check for zero (we do this because we don't '
3079 bne Laddsf$ret | want to return -0 by mistake
3080 bclr IMM (31),d7 | if zero be sure to clear sign
3081 bra Laddsf$ret | if everything OK just return
3082 2:
3083 | The value to be returned is either +/-infty or NaN
3084 andl IMM (0x007fffff),d0 | check for NaN
3085 bne Lf$inop | if mantissa not zero is NaN
3086 bra Lf$infty
3088 Laddsf$ret:
3089 | Normal exit (a and b nonzero, result is not NaN nor +/-infty).
3090 | We have to clear the exception flags (just the exception type).
3091 PICLEA SYM (_fpCCR),a0
3092 movew IMM (0),a0@
3093 orl d7,d0 | put sign bit
3094 #ifndef __mcoldfire__
3095 moveml sp@+,d2-d7 | restore data registers
3096 #else
3097 moveml sp@,d2-d7
3098 | XXX if frame pointer is ever removed, stack pointer must
3099 | be adjusted here.
3100 #endif
3101 unlk a6 | and return
3102 rts
3104 Laddsf$ret$den:
3105 | Return a denormalized number (for addition we don't signal underflow) '
3106 lsrl IMM (1),d0 | remember to shift right back once
3107 bra Laddsf$ret | and return
3109 | Note: when adding two floats of the same sign if either one is
3110 | NaN we return NaN without regard to whether the other is finite or
3111 | not. When subtracting them (i.e., when adding two numbers of
3112 | opposite signs) things are more complicated: if both are INFINITY
3113 | we return NaN, if only one is INFINITY and the other is NaN we return
3114 | NaN, but if it is finite we return INFINITY with the corresponding sign.
3116 Laddsf$nf:
3117 moveq IMM (ADD),d5
3118 | This could be faster but it is not worth the effort, since it is not
3119 | executed very often. We sacrifice speed for clarity here.
3120 movel a6@(8),d0 | get the numbers back (remember that we
3121 movel a6@(12),d1 | did some processing already)
3122 movel IMM (INFINITY),d4 | useful constant (INFINITY)
3123 movel d0,d2 | save sign bits
3124 movel d0,d7 | into d7 as well as we may need the sign
3125 | bit before jumping to LfSinfty
3126 movel d1,d3
3127 bclr IMM (31),d0 | clear sign bits
3128 bclr IMM (31),d1
3129 | We know that one of them is either NaN of +/-INFINITY
3130 | Check for NaN (if either one is NaN return NaN)
3131 cmpl d4,d0 | check first a (d0)
3132 bhi Lf$inop
3133 cmpl d4,d1 | check now b (d1)
3134 bhi Lf$inop
3135 | Now comes the check for +/-INFINITY. We know that both are (maybe not
3136 | finite) numbers, but we have to check if both are infinite whether we
3137 | are adding or subtracting them.
3138 eorl d3,d2 | to check sign bits
3139 bmi 1f
3140 andl IMM (0x80000000),d7 | get (common) sign bit
3141 bra Lf$infty
3142 1:
3143 | We know one (or both) are infinite, so we test for equality between the
3144 | two numbers (if they are equal they have to be infinite both, so we
3145 | return NaN).
3146 cmpl d1,d0 | are both infinite?
3147 beq Lf$inop | if so return NaN
3149 andl IMM (0x80000000),d7 | get a's sign bit '
3150 cmpl d4,d0 | test now for infinity
3151 beq Lf$infty | if a is INFINITY return with this sign
3152 bchg IMM (31),d7 | else we know b is INFINITY and has
3153 bra Lf$infty | the opposite sign
3155 |=============================================================================
3156 | __mulsf3
3157 |=============================================================================
3159 | float __mulsf3(float, float);
3160 FUNC(__mulsf3)
3161 SYM (__mulsf3):
3162 #ifndef __mcoldfire__
3163 link a6,IMM (0)
3164 moveml d2-d7,sp@-
3165 #else
3166 link a6,IMM (-24)
3167 moveml d2-d7,sp@
3168 #endif
3169 movel a6@(8),d0 | get a into d0
3170 movel a6@(12),d1 | and b into d1
3171 movel d0,d7 | d7 will hold the sign of the product
3172 eorl d1,d7 |
3173 andl IMM (0x80000000),d7
3174 movel IMM (INFINITY),d6 | useful constant (+INFINITY)
3175 movel d6,d5 | another (mask for fraction)
3176 notl d5 |
3177 movel IMM (0x00800000),d4 | this is to put hidden bit back
3178 bclr IMM (31),d0 | get rid of a's sign bit '
3179 movel d0,d2 |
3180 beq Lmulsf$a$0 | branch if a is zero
3181 bclr IMM (31),d1 | get rid of b's sign bit '
3182 movel d1,d3 |
3183 beq Lmulsf$b$0 | branch if b is zero
3184 cmpl d6,d0 | is a big?
3185 bhi Lmulsf$inop | if a is NaN return NaN
3186 beq Lmulsf$inf | if a is INFINITY we have to check b
3187 cmpl d6,d1 | now compare b with INFINITY
3188 bhi Lmulsf$inop | is b NaN?
3189 beq Lmulsf$overflow | is b INFINITY?
3190 | Here we have both numbers finite and nonzero (and with no sign bit).
3191 | Now we get the exponents into d2 and d3.
3192 andl d6,d2 | and isolate exponent in d2
3193 beq Lmulsf$a$den | if exponent is zero we have a denormalized
3194 andl d5,d0 | and isolate fraction
3195 orl d4,d0 | and put hidden bit back
3196 swap d2 | I like exponents in the first byte
3197 #ifndef __mcoldfire__
3198 lsrw IMM (7),d2 |
3199 #else
3200 lsrl IMM (7),d2 |
3201 #endif
3202 Lmulsf$1: | number
3203 andl d6,d3 |
3204 beq Lmulsf$b$den |
3205 andl d5,d1 |
3206 orl d4,d1 |
3207 swap d3 |
3208 #ifndef __mcoldfire__
3209 lsrw IMM (7),d3 |
3210 #else
3211 lsrl IMM (7),d3 |
3212 #endif
3213 Lmulsf$2: |
3214 #ifndef __mcoldfire__
3215 addw d3,d2 | add exponents
3216 subw IMM (F_BIAS+1),d2 | and subtract bias (plus one)
3217 #else
3218 addl d3,d2 | add exponents
3219 subl IMM (F_BIAS+1),d2 | and subtract bias (plus one)
3220 #endif
3222 | We are now ready to do the multiplication. The situation is as follows:
3223 | both a and b have bit FLT_MANT_DIG-1 set (even if they were
3224 | denormalized to start with!), which means that in the product
3225 | bit 2*(FLT_MANT_DIG-1) (that is, bit 2*FLT_MANT_DIG-2-32 of the
3226 | high long) is set.
3228 | To do the multiplication let us move the number a little bit around ...
3229 movel d1,d6 | second operand in d6
3230 movel d0,d5 | first operand in d4-d5
3231 movel IMM (0),d4
3232 movel d4,d1 | the sums will go in d0-d1
3233 movel d4,d0
3235 | now bit FLT_MANT_DIG-1 becomes bit 31:
3236 lsll IMM (31-FLT_MANT_DIG+1),d6
3238 | Start the loop (we loop #FLT_MANT_DIG times):
3239 moveq IMM (FLT_MANT_DIG-1),d3
3240 1: addl d1,d1 | shift sum
3241 addxl d0,d0
3242 lsll IMM (1),d6 | get bit bn
3243 bcc 2f | if not set skip sum
3244 addl d5,d1 | add a
3245 addxl d4,d0
3246 2:
3247 #ifndef __mcoldfire__
3248 dbf d3,1b | loop back
3249 #else
3250 subql IMM (1),d3
3251 bpl 1b
3252 #endif
3254 | Now we have the product in d0-d1, with bit (FLT_MANT_DIG - 1) + FLT_MANT_DIG
3255 | (mod 32) of d0 set. The first thing to do now is to normalize it so bit
3256 | FLT_MANT_DIG is set (to do the rounding).
3257 #ifndef __mcoldfire__
3258 rorl IMM (6),d1
3259 swap d1
3260 movew d1,d3
3261 andw IMM (0x03ff),d3
3262 andw IMM (0xfd00),d1
3263 #else
3264 movel d1,d3
3265 lsll IMM (8),d1
3266 addl d1,d1
3267 addl d1,d1
3268 moveq IMM (22),d5
3269 lsrl d5,d3
3270 orl d3,d1
3271 andl IMM (0xfffffd00),d1
3272 #endif
3273 lsll IMM (8),d0
3274 addl d0,d0
3275 addl d0,d0
3276 #ifndef __mcoldfire__
3277 orw d3,d0
3278 #else
3279 orl d3,d0
3280 #endif
3282 moveq IMM (MULTIPLY),d5
3284 btst IMM (FLT_MANT_DIG+1),d0
3285 beq Lround$exit
3286 #ifndef __mcoldfire__
3287 lsrl IMM (1),d0
3288 roxrl IMM (1),d1
3289 addw IMM (1),d2
3290 #else
3291 lsrl IMM (1),d1
3292 btst IMM (0),d0
3293 beq 10f
3294 bset IMM (31),d1
3295 10: lsrl IMM (1),d0
3296 addql IMM (1),d2
3297 #endif
3298 bra Lround$exit
3300 Lmulsf$inop:
3301 moveq IMM (MULTIPLY),d5
3302 bra Lf$inop
3304 Lmulsf$overflow:
3305 moveq IMM (MULTIPLY),d5
3306 bra Lf$overflow
3308 Lmulsf$inf:
3309 moveq IMM (MULTIPLY),d5
3310 | If either is NaN return NaN; else both are (maybe infinite) numbers, so
3311 | return INFINITY with the correct sign (which is in d7).
3312 cmpl d6,d1 | is b NaN?
3313 bhi Lf$inop | if so return NaN
3314 bra Lf$overflow | else return +/-INFINITY
3316 | If either number is zero return zero, unless the other is +/-INFINITY,
3317 | or NaN, in which case we return NaN.
3318 Lmulsf$b$0:
3319 | Here d1 (==b) is zero.
3320 movel a6@(8),d1 | get a again to check for non-finiteness
3321 bra 1f
3322 Lmulsf$a$0:
3323 movel a6@(12),d1 | get b again to check for non-finiteness
3324 1: bclr IMM (31),d1 | clear sign bit
3325 cmpl IMM (INFINITY),d1 | and check for a large exponent
3326 bge Lf$inop | if b is +/-INFINITY or NaN return NaN
3327 movel d7,d0 | else return signed zero
3328 PICLEA SYM (_fpCCR),a0 |
3329 movew IMM (0),a0@ |
3330 #ifndef __mcoldfire__
3331 moveml sp@+,d2-d7 |
3332 #else
3333 moveml sp@,d2-d7
3334 | XXX if frame pointer is ever removed, stack pointer must
3335 | be adjusted here.
3336 #endif
3337 unlk a6 |
3338 rts |
3340 | If a number is denormalized we put an exponent of 1 but do not put the
3341 | hidden bit back into the fraction; instead we shift left until bit 23
3342 | (the hidden bit) is set, adjusting the exponent accordingly. We do this
3343 | to ensure that the product of the fractions is close to 1.
3344 Lmulsf$a$den:
3345 movel IMM (1),d2
3346 andl d5,d0
3347 1: addl d0,d0 | shift a left (until bit 23 is set)
3348 #ifndef __mcoldfire__
3349 subw IMM (1),d2 | and adjust exponent
3350 #else
3351 subql IMM (1),d2 | and adjust exponent
3352 #endif
3353 btst IMM (FLT_MANT_DIG-1),d0
3354 bne Lmulsf$1 |
3355 bra 1b | else loop back
3357 Lmulsf$b$den:
3358 movel IMM (1),d3
3359 andl d5,d1
3360 1: addl d1,d1 | shift b left until bit 23 is set
3361 #ifndef __mcoldfire__
3362 subw IMM (1),d3 | and adjust exponent
3363 #else
3364 subql IMM (1),d3 | and adjust exponent
3365 #endif
3366 btst IMM (FLT_MANT_DIG-1),d1
3367 bne Lmulsf$2 |
3368 bra 1b | else loop back
3370 |=============================================================================
3371 | __divsf3
3372 |=============================================================================
3374 | float __divsf3(float, float);
3375 FUNC(__divsf3)
3376 SYM (__divsf3):
3377 #ifndef __mcoldfire__
3378 link a6,IMM (0)
3379 moveml d2-d7,sp@-
3380 #else
3381 link a6,IMM (-24)
3382 moveml d2-d7,sp@
3383 #endif
3384 movel a6@(8),d0 | get a into d0
3385 movel a6@(12),d1 | and b into d1
3386 movel d0,d7 | d7 will hold the sign of the result
3387 eorl d1,d7 |
3388 andl IMM (0x80000000),d7 |
3389 movel IMM (INFINITY),d6 | useful constant (+INFINITY)
3390 movel d6,d5 | another (mask for fraction)
3391 notl d5 |
3392 movel IMM (0x00800000),d4 | this is to put hidden bit back
3393 bclr IMM (31),d0 | get rid of a's sign bit '
3394 movel d0,d2 |
3395 beq Ldivsf$a$0 | branch if a is zero
3396 bclr IMM (31),d1 | get rid of b's sign bit '
3397 movel d1,d3 |
3398 beq Ldivsf$b$0 | branch if b is zero
3399 cmpl d6,d0 | is a big?
3400 bhi Ldivsf$inop | if a is NaN return NaN
3401 beq Ldivsf$inf | if a is INFINITY we have to check b
3402 cmpl d6,d1 | now compare b with INFINITY
3403 bhi Ldivsf$inop | if b is NaN return NaN
3404 beq Ldivsf$underflow
3405 | Here we have both numbers finite and nonzero (and with no sign bit).
3406 | Now we get the exponents into d2 and d3 and normalize the numbers to
3407 | ensure that the ratio of the fractions is close to 1. We do this by
3408 | making sure that bit #FLT_MANT_DIG-1 (hidden bit) is set.
3409 andl d6,d2 | and isolate exponent in d2
3410 beq Ldivsf$a$den | if exponent is zero we have a denormalized
3411 andl d5,d0 | and isolate fraction
3412 orl d4,d0 | and put hidden bit back
3413 swap d2 | I like exponents in the first byte
3414 #ifndef __mcoldfire__
3415 lsrw IMM (7),d2 |
3416 #else
3417 lsrl IMM (7),d2 |
3418 #endif
3419 Ldivsf$1: |
3420 andl d6,d3 |
3421 beq Ldivsf$b$den |
3422 andl d5,d1 |
3423 orl d4,d1 |
3424 swap d3 |
3425 #ifndef __mcoldfire__
3426 lsrw IMM (7),d3 |
3427 #else
3428 lsrl IMM (7),d3 |
3429 #endif
3430 Ldivsf$2: |
3431 #ifndef __mcoldfire__
3432 subw d3,d2 | subtract exponents
3433 addw IMM (F_BIAS),d2 | and add bias
3434 #else
3435 subl d3,d2 | subtract exponents
3436 addl IMM (F_BIAS),d2 | and add bias
3437 #endif
3439 | We are now ready to do the division. We have prepared things in such a way
3440 | that the ratio of the fractions will be less than 2 but greater than 1/2.
3441 | At this point the registers in use are:
3442 | d0 holds a (first operand, bit FLT_MANT_DIG=0, bit FLT_MANT_DIG-1=1)
3443 | d1 holds b (second operand, bit FLT_MANT_DIG=1)
3444 | d2 holds the difference of the exponents, corrected by the bias
3445 | d7 holds the sign of the ratio
3446 | d4, d5, d6 hold some constants
3447 movel d7,a0 | d6-d7 will hold the ratio of the fractions
3448 movel IMM (0),d6 |
3449 movel d6,d7
3451 moveq IMM (FLT_MANT_DIG+1),d3
3452 1: cmpl d0,d1 | is a < b?
3453 bhi 2f |
3454 bset d3,d6 | set a bit in d6
3455 subl d1,d0 | if a >= b a <-- a-b
3456 beq 3f | if a is zero, exit
3457 2: addl d0,d0 | multiply a by 2
3458 #ifndef __mcoldfire__
3459 dbra d3,1b
3460 #else
3461 subql IMM (1),d3
3462 bpl 1b
3463 #endif
3465 | Now we keep going to set the sticky bit ...
3466 moveq IMM (FLT_MANT_DIG),d3
3467 1: cmpl d0,d1
3468 ble 2f
3469 addl d0,d0
3470 #ifndef __mcoldfire__
3471 dbra d3,1b
3472 #else
3473 subql IMM(1),d3
3474 bpl 1b
3475 #endif
3476 movel IMM (0),d1
3477 bra 3f
3478 2: movel IMM (0),d1
3479 #ifndef __mcoldfire__
3480 subw IMM (FLT_MANT_DIG),d3
3481 addw IMM (31),d3
3482 #else
3483 subl IMM (FLT_MANT_DIG),d3
3484 addl IMM (31),d3
3485 #endif
3486 bset d3,d1
3487 3:
3488 movel d6,d0 | put the ratio in d0-d1
3489 movel a0,d7 | get sign back
3491 | Because of the normalization we did before we are guaranteed that
3492 | d0 is smaller than 2^26 but larger than 2^24. Thus bit 26 is not set,
3493 | bit 25 could be set, and if it is not set then bit 24 is necessarily set.
3494 btst IMM (FLT_MANT_DIG+1),d0
3495 beq 1f | if it is not set, then bit 24 is set
3496 lsrl IMM (1),d0 |
3497 #ifndef __mcoldfire__
3498 addw IMM (1),d2 |
3499 #else
3500 addl IMM (1),d2 |
3501 #endif
3502 1:
3503 | Now round, check for over- and underflow, and exit.
3504 moveq IMM (DIVIDE),d5
3505 bra Lround$exit
3507 Ldivsf$inop:
3508 moveq IMM (DIVIDE),d5
3509 bra Lf$inop
3511 Ldivsf$overflow:
3512 moveq IMM (DIVIDE),d5
3513 bra Lf$overflow
3515 Ldivsf$underflow:
3516 moveq IMM (DIVIDE),d5
3517 bra Lf$underflow
3519 Ldivsf$a$0:
3520 moveq IMM (DIVIDE),d5
3521 | If a is zero check to see whether b is zero also. In that case return
3522 | NaN; then check if b is NaN, and return NaN also in that case. Else
3523 | return a properly signed zero.
3524 andl IMM (0x7fffffff),d1 | clear sign bit and test b
3525 beq Lf$inop | if b is also zero return NaN
3526 cmpl IMM (INFINITY),d1 | check for NaN
3527 bhi Lf$inop |
3528 movel d7,d0 | else return signed zero
3529 PICLEA SYM (_fpCCR),a0 |
3530 movew IMM (0),a0@ |
3531 #ifndef __mcoldfire__
3532 moveml sp@+,d2-d7 |
3533 #else
3534 moveml sp@,d2-d7 |
3535 | XXX if frame pointer is ever removed, stack pointer must
3536 | be adjusted here.
3537 #endif
3538 unlk a6 |
3539 rts |
3541 Ldivsf$b$0:
3542 moveq IMM (DIVIDE),d5
3543 | If we got here a is not zero. Check if a is NaN; in that case return NaN,
3544 | else return +/-INFINITY. Remember that a is in d0 with the sign bit
3545 | cleared already.
3546 cmpl IMM (INFINITY),d0 | compare d0 with INFINITY
3547 bhi Lf$inop | if larger it is NaN
3548 bra Lf$div$0 | else signal DIVIDE_BY_ZERO
3550 Ldivsf$inf:
3551 moveq IMM (DIVIDE),d5
3552 | If a is INFINITY we have to check b
3553 cmpl IMM (INFINITY),d1 | compare b with INFINITY
3554 bge Lf$inop | if b is NaN or INFINITY return NaN
3555 bra Lf$overflow | else return overflow
3557 | If a number is denormalized we put an exponent of 1 but do not put the
3558 | bit back into the fraction.
3559 Ldivsf$a$den:
3560 movel IMM (1),d2
3561 andl d5,d0
3562 1: addl d0,d0 | shift a left until bit FLT_MANT_DIG-1 is set
3563 #ifndef __mcoldfire__
3564 subw IMM (1),d2 | and adjust exponent
3565 #else
3566 subl IMM (1),d2 | and adjust exponent
3567 #endif
3568 btst IMM (FLT_MANT_DIG-1),d0
3569 bne Ldivsf$1
3570 bra 1b
3572 Ldivsf$b$den:
3573 movel IMM (1),d3
3574 andl d5,d1
3575 1: addl d1,d1 | shift b left until bit FLT_MANT_DIG is set
3576 #ifndef __mcoldfire__
3577 subw IMM (1),d3 | and adjust exponent
3578 #else
3579 subl IMM (1),d3 | and adjust exponent
3580 #endif
3581 btst IMM (FLT_MANT_DIG-1),d1
3582 bne Ldivsf$2
3583 bra 1b
3585 Lround$exit:
3586 | This is a common exit point for __mulsf3 and __divsf3.
3588 | First check for underlow in the exponent:
3589 #ifndef __mcoldfire__
3590 cmpw IMM (-FLT_MANT_DIG-1),d2
3591 #else
3592 cmpl IMM (-FLT_MANT_DIG-1),d2
3593 #endif
3594 blt Lf$underflow
3595 | It could happen that the exponent is less than 1, in which case the
3596 | number is denormalized. In this case we shift right and adjust the
3597 | exponent until it becomes 1 or the fraction is zero (in the latter case
3598 | we signal underflow and return zero).
3599 movel IMM (0),d6 | d6 is used temporarily
3600 #ifndef __mcoldfire__
3601 cmpw IMM (1),d2 | if the exponent is less than 1 we
3602 #else
3603 cmpl IMM (1),d2 | if the exponent is less than 1 we
3604 #endif
3605 bge 2f | have to shift right (denormalize)
3606 1:
3607 #ifndef __mcoldfire__
3608 addw IMM (1),d2 | adjust the exponent
3609 lsrl IMM (1),d0 | shift right once
3610 roxrl IMM (1),d1 |
3611 roxrl IMM (1),d6 | d6 collect bits we would lose otherwise
3612 cmpw IMM (1),d2 | is the exponent 1 already?
3613 #else
3614 addql IMM (1),d2 | adjust the exponent
3615 lsrl IMM (1),d6
3616 btst IMM (0),d1
3617 beq 11f
3618 bset IMM (31),d6
3619 11: lsrl IMM (1),d1
3620 btst IMM (0),d0
3621 beq 10f
3622 bset IMM (31),d1
3623 10: lsrl IMM (1),d0
3624 cmpl IMM (1),d2 | is the exponent 1 already?
3625 #endif
3626 beq 2f | if not loop back
3627 bra 1b |
3628 bra Lf$underflow | safety check, shouldn't execute '
3629 2: orl d6,d1 | this is a trick so we don't lose '
3630 | the extra bits which were flushed right
3631 | Now call the rounding routine (which takes care of denormalized numbers):
3632 lea pc@(Lround$0),a0 | to return from rounding routine
3633 PICLEA SYM (_fpCCR),a1 | check the rounding mode
3634 #ifdef __mcoldfire__
3635 clrl d6
3636 #endif
3637 movew a1@(6),d6 | rounding mode in d6
3638 beq Lround$to$nearest
3639 #ifndef __mcoldfire__
3640 cmpw IMM (ROUND_TO_PLUS),d6
3641 #else
3642 cmpl IMM (ROUND_TO_PLUS),d6
3643 #endif
3644 bhi Lround$to$minus
3645 blt Lround$to$zero
3646 bra Lround$to$plus
3647 Lround$0:
3648 | Here we have a correctly rounded result (either normalized or denormalized).
3650 | Here we should have either a normalized number or a denormalized one, and
3651 | the exponent is necessarily larger or equal to 1 (so we don't have to '
3652 | check again for underflow!). We have to check for overflow or for a
3653 | denormalized number (which also signals underflow).
3654 | Check for overflow (i.e., exponent >= 255).
3655 #ifndef __mcoldfire__
3656 cmpw IMM (0x00ff),d2
3657 #else
3658 cmpl IMM (0x00ff),d2
3659 #endif
3660 bge Lf$overflow
3661 | Now check for a denormalized number (exponent==0).
3662 movew d2,d2
3663 beq Lf$den
3664 1:
3665 | Put back the exponents and sign and return.
3666 #ifndef __mcoldfire__
3667 lslw IMM (7),d2 | exponent back to fourth byte
3668 #else
3669 lsll IMM (7),d2 | exponent back to fourth byte
3670 #endif
3671 bclr IMM (FLT_MANT_DIG-1),d0
3672 swap d0 | and put back exponent
3673 #ifndef __mcoldfire__
3674 orw d2,d0 |
3675 #else
3676 orl d2,d0
3677 #endif
3678 swap d0 |
3679 orl d7,d0 | and sign also
3681 PICLEA SYM (_fpCCR),a0
3682 movew IMM (0),a0@
3683 #ifndef __mcoldfire__
3684 moveml sp@+,d2-d7
3685 #else
3686 moveml sp@,d2-d7
3687 | XXX if frame pointer is ever removed, stack pointer must
3688 | be adjusted here.
3689 #endif
3690 unlk a6
3691 rts
3693 |=============================================================================
3694 | __negsf2
3695 |=============================================================================
3697 | This is trivial and could be shorter if we didn't bother checking for NaN '
3698 | and +/-INFINITY.
3700 | float __negsf2(float);
3701 FUNC(__negsf2)
3702 SYM (__negsf2):
3703 #ifndef __mcoldfire__
3704 link a6,IMM (0)
3705 moveml d2-d7,sp@-
3706 #else
3707 link a6,IMM (-24)
3708 moveml d2-d7,sp@
3709 #endif
3710 moveq IMM (NEGATE),d5
3711 movel a6@(8),d0 | get number to negate in d0
3712 bchg IMM (31),d0 | negate
3713 movel d0,d1 | make a positive copy
3714 bclr IMM (31),d1 |
3715 tstl d1 | check for zero
3716 beq 2f | if zero (either sign) return +zero
3717 cmpl IMM (INFINITY),d1 | compare to +INFINITY
3718 blt 1f |
3719 bhi Lf$inop | if larger (fraction not zero) is NaN
3720 movel d0,d7 | else get sign and return INFINITY
3721 andl IMM (0x80000000),d7
3722 bra Lf$infty
3723 1: PICLEA SYM (_fpCCR),a0
3724 movew IMM (0),a0@
3725 #ifndef __mcoldfire__
3726 moveml sp@+,d2-d7
3727 #else
3728 moveml sp@,d2-d7
3729 | XXX if frame pointer is ever removed, stack pointer must
3730 | be adjusted here.
3731 #endif
3732 unlk a6
3733 rts
3734 2: bclr IMM (31),d0
3735 bra 1b
3737 |=============================================================================
3738 | __cmpsf2
3739 |=============================================================================
3741 GREATER = 1
3742 LESS = -1
3743 EQUAL = 0
3745 | int __cmpsf2_internal(float, float, int);
3746 SYM (__cmpsf2_internal):
3747 #ifndef __mcoldfire__
3748 link a6,IMM (0)
3749 moveml d2-d7,sp@- | save registers
3750 #else
3751 link a6,IMM (-24)
3752 moveml d2-d7,sp@
3753 #endif
3754 moveq IMM (COMPARE),d5
3755 movel a6@(8),d0 | get first operand
3756 movel a6@(12),d1 | get second operand
3757 | Check if either is NaN, and in that case return garbage and signal
3758 | INVALID_OPERATION. Check also if either is zero, and clear the signs
3759 | if necessary.
3760 movel d0,d6
3761 andl IMM (0x7fffffff),d0
3762 beq Lcmpsf$a$0
3763 cmpl IMM (0x7f800000),d0
3764 bhi Lcmpf$inop
3765 Lcmpsf$1:
3766 movel d1,d7
3767 andl IMM (0x7fffffff),d1
3768 beq Lcmpsf$b$0
3769 cmpl IMM (0x7f800000),d1
3770 bhi Lcmpf$inop
3771 Lcmpsf$2:
3772 | Check the signs
3773 eorl d6,d7
3774 bpl 1f
3775 | If the signs are not equal check if a >= 0
3776 tstl d6
3777 bpl Lcmpsf$a$gt$b | if (a >= 0 && b < 0) => a > b
3778 bmi Lcmpsf$b$gt$a | if (a < 0 && b >= 0) => a < b
3779 1:
3780 | If the signs are equal check for < 0
3781 tstl d6
3782 bpl 1f
3783 | If both are negative exchange them
3784 #ifndef __mcoldfire__
3785 exg d0,d1
3786 #else
3787 movel d0,d7
3788 movel d1,d0
3789 movel d7,d1
3790 #endif
3791 1:
3792 | Now that they are positive we just compare them as longs (does this also
3793 | work for denormalized numbers?).
3794 cmpl d0,d1
3795 bhi Lcmpsf$b$gt$a | |b| > |a|
3796 bne Lcmpsf$a$gt$b | |b| < |a|
3797 | If we got here a == b.
3798 movel IMM (EQUAL),d0
3799 #ifndef __mcoldfire__
3800 moveml sp@+,d2-d7 | put back the registers
3801 #else
3802 moveml sp@,d2-d7
3803 #endif
3804 unlk a6
3805 rts
3806 Lcmpsf$a$gt$b:
3807 movel IMM (GREATER),d0
3808 #ifndef __mcoldfire__
3809 moveml sp@+,d2-d7 | put back the registers
3810 #else
3811 moveml sp@,d2-d7
3812 | XXX if frame pointer is ever removed, stack pointer must
3813 | be adjusted here.
3814 #endif
3815 unlk a6
3816 rts
3817 Lcmpsf$b$gt$a:
3818 movel IMM (LESS),d0
3819 #ifndef __mcoldfire__
3820 moveml sp@+,d2-d7 | put back the registers
3821 #else
3822 moveml sp@,d2-d7
3823 | XXX if frame pointer is ever removed, stack pointer must
3824 | be adjusted here.
3825 #endif
3826 unlk a6
3827 rts
3829 Lcmpsf$a$0:
3830 bclr IMM (31),d6
3831 bra Lcmpsf$1
3832 Lcmpsf$b$0:
3833 bclr IMM (31),d7
3834 bra Lcmpsf$2
3836 Lcmpf$inop:
3837 movl a6@(16),d0
3838 moveq IMM (INEXACT_RESULT+INVALID_OPERATION),d7
3839 moveq IMM (SINGLE_FLOAT),d6
3840 PICJUMP $_exception_handler
3842 | int __cmpsf2(float, float);
3843 FUNC(__cmpsf2)
3844 SYM (__cmpsf2):
3845 link a6,IMM (0)
3846 pea 1
3847 movl a6@(12),sp@-
3848 movl a6@(8),sp@-
3849 PICCALL SYM (__cmpsf2_internal)
3850 unlk a6
3851 rts
3853 |=============================================================================
3854 | rounding routines
3855 |=============================================================================
3857 | The rounding routines expect the number to be normalized in registers
3858 | d0-d1, with the exponent in register d2. They assume that the
3859 | exponent is larger or equal to 1. They return a properly normalized number
3860 | if possible, and a denormalized number otherwise. The exponent is returned
3861 | in d2.
3863 Lround$to$nearest:
3864 | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
3865 | Here we assume that the exponent is not too small (this should be checked
3866 | before entering the rounding routine), but the number could be denormalized.
3868 | Check for denormalized numbers:
3869 1: btst IMM (FLT_MANT_DIG),d0
3870 bne 2f | if set the number is normalized
3871 | Normalize shifting left until bit #FLT_MANT_DIG is set or the exponent
3872 | is one (remember that a denormalized number corresponds to an
3873 | exponent of -F_BIAS+1).
3874 #ifndef __mcoldfire__
3875 cmpw IMM (1),d2 | remember that the exponent is at least one
3876 #else
3877 cmpl IMM (1),d2 | remember that the exponent is at least one
3878 #endif
3879 beq 2f | an exponent of one means denormalized
3880 addl d1,d1 | else shift and adjust the exponent
3881 addxl d0,d0 |
3882 #ifndef __mcoldfire__
3883 dbra d2,1b |
3884 #else
3885 subql IMM (1),d2
3886 bpl 1b
3887 #endif
3888 2:
3889 | Now round: we do it as follows: after the shifting we can write the
3890 | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
3891 | If delta < 1, do nothing. If delta > 1, add 1 to f.
3892 | If delta == 1, we make sure the rounded number will be even (odd?)
3893 | (after shifting).
3894 btst IMM (0),d0 | is delta < 1?
3895 beq 2f | if so, do not do anything
3896 tstl d1 | is delta == 1?
3897 bne 1f | if so round to even
3898 movel d0,d1 |
3899 andl IMM (2),d1 | bit 1 is the last significant bit
3900 addl d1,d0 |
3901 bra 2f |
3902 1: movel IMM (1),d1 | else add 1
3903 addl d1,d0 |
3904 | Shift right once (because we used bit #FLT_MANT_DIG!).
3905 2: lsrl IMM (1),d0
3906 | Now check again bit #FLT_MANT_DIG (rounding could have produced a
3907 | 'fraction overflow' ...).
3908 btst IMM (FLT_MANT_DIG),d0
3909 beq 1f
3910 lsrl IMM (1),d0
3911 #ifndef __mcoldfire__
3912 addw IMM (1),d2
3913 #else
3914 addql IMM (1),d2
3915 #endif
3916 1:
3917 | If bit #FLT_MANT_DIG-1 is clear we have a denormalized number, so we
3918 | have to put the exponent to zero and return a denormalized number.
3919 btst IMM (FLT_MANT_DIG-1),d0
3920 beq 1f
3921 jmp a0@
3922 1: movel IMM (0),d2
3923 jmp a0@
3925 Lround$to$zero:
3926 Lround$to$plus:
3927 Lround$to$minus:
3928 jmp a0@
3929 #endif /* L_float */
3931 | gcc expects the routines __eqdf2, __nedf2, __gtdf2, __gedf2,
3932 | __ledf2, __ltdf2 to all return the same value as a direct call to
3933 | __cmpdf2 would. In this implementation, each of these routines
3934 | simply calls __cmpdf2. It would be more efficient to give the
3935 | __cmpdf2 routine several names, but separating them out will make it
3936 | easier to write efficient versions of these routines someday.
3937 | If the operands recompare unordered unordered __gtdf2 and __gedf2 return -1.
3938 | The other routines return 1.
3940 #ifdef L_eqdf2
3941 .text
3942 FUNC(__eqdf2)
3943 .globl SYM (__eqdf2)
3944 SYM (__eqdf2):
3945 link a6,IMM (0)
3946 pea 1
3947 movl a6@(20),sp@-
3948 movl a6@(16),sp@-
3949 movl a6@(12),sp@-
3950 movl a6@(8),sp@-
3951 PICCALL SYM (__cmpdf2_internal)
3952 unlk a6
3953 rts
3954 #endif /* L_eqdf2 */
3956 #ifdef L_nedf2
3957 .text
3958 FUNC(__nedf2)
3959 .globl SYM (__nedf2)
3960 SYM (__nedf2):
3961 link a6,IMM (0)
3962 pea 1
3963 movl a6@(20),sp@-
3964 movl a6@(16),sp@-
3965 movl a6@(12),sp@-
3966 movl a6@(8),sp@-
3967 PICCALL SYM (__cmpdf2_internal)
3968 unlk a6
3969 rts
3970 #endif /* L_nedf2 */
3972 #ifdef L_gtdf2
3973 .text
3974 FUNC(__gtdf2)
3975 .globl SYM (__gtdf2)
3976 SYM (__gtdf2):
3977 link a6,IMM (0)
3978 pea -1
3979 movl a6@(20),sp@-
3980 movl a6@(16),sp@-
3981 movl a6@(12),sp@-
3982 movl a6@(8),sp@-
3983 PICCALL SYM (__cmpdf2_internal)
3984 unlk a6
3985 rts
3986 #endif /* L_gtdf2 */
3988 #ifdef L_gedf2
3989 .text
3990 FUNC(__gedf2)
3991 .globl SYM (__gedf2)
3992 SYM (__gedf2):
3993 link a6,IMM (0)
3994 pea -1
3995 movl a6@(20),sp@-
3996 movl a6@(16),sp@-
3997 movl a6@(12),sp@-
3998 movl a6@(8),sp@-
3999 PICCALL SYM (__cmpdf2_internal)
4000 unlk a6
4001 rts
4002 #endif /* L_gedf2 */
4004 #ifdef L_ltdf2
4005 .text
4006 FUNC(__ltdf2)
4007 .globl SYM (__ltdf2)
4008 SYM (__ltdf2):
4009 link a6,IMM (0)
4010 pea 1
4011 movl a6@(20),sp@-
4012 movl a6@(16),sp@-
4013 movl a6@(12),sp@-
4014 movl a6@(8),sp@-
4015 PICCALL SYM (__cmpdf2_internal)
4016 unlk a6
4017 rts
4018 #endif /* L_ltdf2 */
4020 #ifdef L_ledf2
4021 .text
4022 FUNC(__ledf2)
4023 .globl SYM (__ledf2)
4024 SYM (__ledf2):
4025 link a6,IMM (0)
4026 pea 1
4027 movl a6@(20),sp@-
4028 movl a6@(16),sp@-
4029 movl a6@(12),sp@-
4030 movl a6@(8),sp@-
4031 PICCALL SYM (__cmpdf2_internal)
4032 unlk a6
4033 rts
4034 #endif /* L_ledf2 */
4036 | The comments above about __eqdf2, et. al., also apply to __eqsf2,
4037 | et. al., except that the latter call __cmpsf2 rather than __cmpdf2.
4039 #ifdef L_eqsf2
4040 .text
4041 FUNC(__eqsf2)
4042 .globl SYM (__eqsf2)
4043 SYM (__eqsf2):
4044 link a6,IMM (0)
4045 pea 1
4046 movl a6@(12),sp@-
4047 movl a6@(8),sp@-
4048 PICCALL SYM (__cmpsf2_internal)
4049 unlk a6
4050 rts
4051 #endif /* L_eqsf2 */
4053 #ifdef L_nesf2
4054 .text
4055 FUNC(__nesf2)
4056 .globl SYM (__nesf2)
4057 SYM (__nesf2):
4058 link a6,IMM (0)
4059 pea 1
4060 movl a6@(12),sp@-
4061 movl a6@(8),sp@-
4062 PICCALL SYM (__cmpsf2_internal)
4063 unlk a6
4064 rts
4065 #endif /* L_nesf2 */
4067 #ifdef L_gtsf2
4068 .text
4069 FUNC(__gtsf2)
4070 .globl SYM (__gtsf2)
4071 SYM (__gtsf2):
4072 link a6,IMM (0)
4073 pea -1
4074 movl a6@(12),sp@-
4075 movl a6@(8),sp@-
4076 PICCALL SYM (__cmpsf2_internal)
4077 unlk a6
4078 rts
4079 #endif /* L_gtsf2 */
4081 #ifdef L_gesf2
4082 .text
4083 FUNC(__gesf2)
4084 .globl SYM (__gesf2)
4085 SYM (__gesf2):
4086 link a6,IMM (0)
4087 pea -1
4088 movl a6@(12),sp@-
4089 movl a6@(8),sp@-
4090 PICCALL SYM (__cmpsf2_internal)
4091 unlk a6
4092 rts
4093 #endif /* L_gesf2 */
4095 #ifdef L_ltsf2
4096 .text
4097 FUNC(__ltsf2)
4098 .globl SYM (__ltsf2)
4099 SYM (__ltsf2):
4100 link a6,IMM (0)
4101 pea 1
4102 movl a6@(12),sp@-
4103 movl a6@(8),sp@-
4104 PICCALL SYM (__cmpsf2_internal)
4105 unlk a6
4106 rts
4107 #endif /* L_ltsf2 */
4109 #ifdef L_lesf2
4110 .text
4111 FUNC(__lesf2)
4112 .globl SYM (__lesf2)
4113 SYM (__lesf2):
4114 link a6,IMM (0)
4115 pea 1
4116 movl a6@(12),sp@-
4117 movl a6@(8),sp@-
4118 PICCALL SYM (__cmpsf2_internal)
4119 unlk a6
4120 rts
4121 #endif /* L_lesf2 */
4123 #if defined (__ELF__) && defined (__linux__)
4124 /* Make stack non-executable for ELF linux targets. */
4125 .section .note.GNU-stack,"",@progbits
4126 #endif