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