![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to lb1sf68.asm CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [GCC 2.x] / gcc / config / m68k|
Return to lb1sf68.asm CVS log | Up to [GCC 2.x] / gcc / config / m68k |
1.1 root 1: /* libgcc1 routines for 68000 w/o floating-point hardware. */
2: /* Copyright (C) 1994 Free Software Foundation, Inc.
3:
4: This file is free software; you can redistribute it and/or modify it
5: under the terms of the GNU General Public License as published by the
6: Free Software Foundation; either version 2, or (at your option) any
7: later version.
8:
9: In addition to the permissions in the GNU General Public License, the
10: Free Software Foundation gives you unlimited permission to link the
11: compiled version of this file with other programs, and to distribute
12: those programs without any restriction coming from the use of this
13: file. (The General Public License restrictions do apply in other
14: respects; for example, they cover modification of the file, and
15: distribution when not linked into another program.)
16:
17: This file is distributed in the hope that it will be useful, but
18: WITHOUT ANY WARRANTY; without even the implied warranty of
19: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20: General Public License for more details.
21:
22: You should have received a copy of the GNU General Public License
23: along with this program; see the file COPYING. If not, write to
24: the Free Software Foundation, 59 Temple Place - Suite 330,
25: Boston, MA 02111-1307, USA. */
26:
27: /* As a special exception, if you link this library with files
28: compiled with GCC to produce an executable, this does not cause
29: the resulting executable to be covered by the GNU General Public License.
30: This exception does not however invalidate any other reasons why
31: the executable file might be covered by the GNU General Public License. */
32:
33: /* Use this one for any 680x0; assumes no floating point hardware.
34: The trailing " '" appearing on some lines is for ANSI preprocessors. Yuk.
35: Some of this code comes from MINIX, via the folks at ericsson.
36: D. V. Henkel-Wallace ([email protected]) Fete Bastille, 1992
37: */
38:
39: /* These are predefined by new versions of GNU cpp. */
40:
41: #ifndef __USER_LABEL_PREFIX__
42: #define __USER_LABEL_PREFIX__ _
43: #endif
44:
45: #ifndef __REGISTER_PREFIX__
46: #define __REGISTER_PREFIX__
47: #endif
48:
49: #ifndef __IMMEDIATE_PREFIX__
50: #define __IMMEDIATE_PREFIX__ #
51: #endif
52:
53: /* ANSI concatenation macros. */
54:
55: #define CONCAT1(a, b) CONCAT2(a, b)
56: #define CONCAT2(a, b) a ## b
57:
58: /* Use the right prefix for global labels. */
59:
60: #define SYM(x) CONCAT1 (__USER_LABEL_PREFIX__, x)
61:
62: /* Use the right prefix for registers. */
63:
64: #define REG(x) CONCAT1 (__REGISTER_PREFIX__, x)
65:
66: /* Use the right prefix for immediate values. */
67:
68: #define IMM(x) CONCAT1 (__IMMEDIATE_PREFIX__, x)
69:
70: #define d0 REG (d0)
71: #define d1 REG (d1)
72: #define d2 REG (d2)
73: #define d3 REG (d3)
74: #define d4 REG (d4)
75: #define d5 REG (d5)
76: #define d6 REG (d6)
77: #define d7 REG (d7)
78: #define a0 REG (a0)
79: #define a1 REG (a1)
80: #define a2 REG (a2)
81: #define a3 REG (a3)
82: #define a4 REG (a4)
83: #define a5 REG (a5)
84: #define a6 REG (a6)
85: #define fp REG (fp)
86: #define sp REG (sp)
87:
88: #ifdef L_floatex
89:
90: | This is an attempt at a decent floating point (single, double and
91: | extended double) code for the GNU C compiler. It should be easy to
92: | adapt to other compilers (but beware of the local labels!).
93:
94: | Starting date: 21 October, 1990
95:
96: | It is convenient to introduce the notation (s,e,f) for a floating point
97: | number, where s=sign, e=exponent, f=fraction. We will call a floating
98: | point number fpn to abbreviate, independently of the precision.
99: | Let MAX_EXP be in each case the maximum exponent (255 for floats, 1023
100: | for doubles and 16383 for long doubles). We then have the following
101: | different cases:
102: | 1. Normalized fpns have 0 < e < MAX_EXP. They correspond to
103: | (-1)^s x 1.f x 2^(e-bias-1).
104: | 2. Denormalized fpns have e=0. They correspond to numbers of the form
105: | (-1)^s x 0.f x 2^(-bias).
106: | 3. +/-INFINITY have e=MAX_EXP, f=0.
107: | 4. Quiet NaN (Not a Number) have all bits set.
108: | 5. Signaling NaN (Not a Number) have s=0, e=MAX_EXP, f=1.
109:
110: |=============================================================================
111: | exceptions
112: |=============================================================================
113:
114: | This is the floating point condition code register (_fpCCR):
115: |
116: | struct {
117: | short _exception_bits;
118: | short _trap_enable_bits;
119: | short _sticky_bits;
120: | short _rounding_mode;
121: | short _format;
122: | short _last_operation;
123: | union {
124: | float sf;
125: | double df;
126: | } _operand1;
127: | union {
128: | float sf;
129: | double df;
130: | } _operand2;
131: | } _fpCCR;
132:
133: .data
134: .even
135:
136: .globl SYM (_fpCCR)
137:
138: SYM (_fpCCR):
139: __exception_bits:
140: .word 0
141: __trap_enable_bits:
142: .word 0
143: __sticky_bits:
144: .word 0
145: __rounding_mode:
146: .word ROUND_TO_NEAREST
147: __format:
148: .word NIL
149: __last_operation:
150: .word NOOP
151: __operand1:
152: .long 0
153: .long 0
154: __operand2:
155: .long 0
156: .long 0
157:
158: | Offsets:
159: EBITS = __exception_bits - SYM (_fpCCR)
160: TRAPE = __trap_enable_bits - SYM (_fpCCR)
161: STICK = __sticky_bits - SYM (_fpCCR)
162: ROUND = __rounding_mode - SYM (_fpCCR)
163: FORMT = __format - SYM (_fpCCR)
164: LASTO = __last_operation - SYM (_fpCCR)
165: OPER1 = __operand1 - SYM (_fpCCR)
166: OPER2 = __operand2 - SYM (_fpCCR)
167:
168: | The following exception types are supported:
169: INEXACT_RESULT = 0x0001
170: UNDERFLOW = 0x0002
171: OVERFLOW = 0x0004
172: DIVIDE_BY_ZERO = 0x0008
173: INVALID_OPERATION = 0x0010
174:
175: | The allowed rounding modes are:
176: UNKNOWN = -1
177: ROUND_TO_NEAREST = 0 | round result to nearest representable value
178: ROUND_TO_ZERO = 1 | round result towards zero
179: ROUND_TO_PLUS = 2 | round result towards plus infinity
180: ROUND_TO_MINUS = 3 | round result towards minus infinity
181:
182: | The allowed values of format are:
183: NIL = 0
184: SINGLE_FLOAT = 1
185: DOUBLE_FLOAT = 2
186: LONG_FLOAT = 3
187:
188: | The allowed values for the last operation are:
189: NOOP = 0
190: ADD = 1
191: MULTIPLY = 2
192: DIVIDE = 3
193: NEGATE = 4
194: COMPARE = 5
195: EXTENDSFDF = 6
196: TRUNCDFSF = 7
197:
198: |=============================================================================
199: | __clear_sticky_bits
200: |=============================================================================
201:
202: | The sticky bits are normally not cleared (thus the name), whereas the
203: | exception type and exception value reflect the last computation.
204: | This routine is provided to clear them (you can also write to _fpCCR,
205: | since it is globally visible).
206:
207: .globl SYM (__clear_sticky_bit)
208:
209: .text
210: .even
211:
212: | void __clear_sticky_bits(void);
213: SYM (__clear_sticky_bit):
214: lea SYM (_fpCCR),a0
215: movew IMM (0),a0@(STICK)
216: rts
217:
218: |=============================================================================
219: | $_exception_handler
220: |=============================================================================
221:
222: .globl $_exception_handler
223:
224: .text
225: .even
226:
227: | This is the common exit point if an exception occurs.
228: | NOTE: it is NOT callable from C!
229: | It expects the exception type in d7, the format (SINGLE_FLOAT,
230: | DOUBLE_FLOAT or LONG_FLOAT) in d6, and the last operation code in d5.
231: | It sets the corresponding exception and sticky bits, and the format.
232: | Depending on the format if fills the corresponding slots for the
233: | operands which produced the exception (all this information is provided
234: | so if you write your own exception handlers you have enough information
235: | to deal with the problem).
236: | Then checks to see if the corresponding exception is trap-enabled,
237: | in which case it pushes the address of _fpCCR and traps through
238: | trap FPTRAP (15 for the moment).
239:
240: FPTRAP = 15
241:
242: $_exception_handler:
243: lea SYM (_fpCCR),a0
244: movew d7,a0@(EBITS) | set __exception_bits
245: orw d7,a0@(STICK) | and __sticky_bits
246: movew d6,a0@(FORMT) | and __format
247: movew d5,a0@(LASTO) | and __last_operation
248:
249: | Now put the operands in place:
250: cmpw IMM (SINGLE_FLOAT),d6
251: beq 1f
252: movel a6@(8),a0@(OPER1)
253: movel a6@(12),a0@(OPER1+4)
254: movel a6@(16),a0@(OPER2)
255: movel a6@(20),a0@(OPER2+4)
256: bra 2f
257: 1: movel a6@(8),a0@(OPER1)
258: movel a6@(12),a0@(OPER2)
259: 2:
260: | And check whether the exception is trap-enabled:
261: andw a0@(TRAPE),d7 | is exception trap-enabled?
262: beq 1f | no, exit
263: pea SYM (_fpCCR) | yes, push address of _fpCCR
264: trap IMM (FPTRAP) | and trap
265: 1: moveml sp@+,d2-d7 | restore data registers
266: unlk a6 | and return
267: rts
268: #endif /* L_floatex */
269:
270: #ifdef L_mulsi3
271: .text
272: .proc
273: .globl SYM (__mulsi3)
274: SYM (__mulsi3):
275: movew sp@(4), d0 /* x0 -> d0 */
276: muluw sp@(10), d0 /* x0*y1 */
277: movew sp@(6), d1 /* x1 -> d1 */
278: muluw sp@(8), d1 /* x1*y0 */
279: addw d1, d0
280: swap d0
281: clrw d0
282: movew sp@(6), d1 /* x1 -> d1 */
283: muluw sp@(10), d1 /* x1*y1 */
284: addl d1, d0
285:
286: rts
287: #endif /* L_mulsi3 */
288:
289: #ifdef L_udivsi3
290: .text
291: .proc
292: .globl SYM (__udivsi3)
293: SYM (__udivsi3):
294: movel d2, sp@-
295: movel sp@(12), d1 /* d1 = divisor */
296: movel sp@(8), d0 /* d0 = dividend */
297:
298: cmpl IMM (0x10000), d1 /* divisor >= 2 ^ 16 ? */
299: jcc L3 /* then try next algorithm */
300: movel d0, d2
301: clrw d2
302: swap d2
303: divu d1, d2 /* high quotient in lower word */
304: movew d2, d0 /* save high quotient */
305: swap d0
306: movew sp@(10), d2 /* get low dividend + high rest */
307: divu d1, d2 /* low quotient */
308: movew d2, d0
309: jra L6
310:
311: L3: movel d1, d2 /* use d2 as divisor backup */
312: L4: lsrl IMM (1), d1 /* shift divisor */
313: lsrl IMM (1), d0 /* shift dividend */
314: cmpl IMM (0x10000), d1 /* still divisor >= 2 ^ 16 ? */
315: jcc L4
316: divu d1, d0 /* now we have 16 bit divisor */
317: andl IMM (0xffff), d0 /* mask out divisor, ignore remainder */
318:
319: /* Multiply the 16 bit tentative quotient with the 32 bit divisor. Because of
320: the operand ranges, this might give a 33 bit product. If this product is
321: greater than the dividend, the tentative quotient was too large. */
322: movel d2, d1
323: mulu d0, d1 /* low part, 32 bits */
324: swap d2
325: mulu d0, d2 /* high part, at most 17 bits */
326: swap d2 /* align high part with low part */
327: btst IMM (0), d2 /* high part 17 bits? */
328: jne L5 /* if 17 bits, quotient was too large */
329: addl d2, d1 /* add parts */
330: jcs L5 /* if sum is 33 bits, quotient was too large */
331: cmpl sp@(8), d1 /* compare the sum with the dividend */
332: jls L6 /* if sum > dividend, quotient was too large */
333: L5: subql IMM (1), d0 /* adjust quotient */
334:
335: L6: movel sp@+, d2
336: rts
337: #endif /* L_udivsi3 */
338:
339: #ifdef L_divsi3
340: .text
341: .proc
342: .globl SYM (__divsi3)
343: SYM (__divsi3):
344: movel d2, sp@-
345:
346: moveb IMM (1), d2 /* sign of result stored in d2 (=1 or =-1) */
347: movel sp@(12), d1 /* d1 = divisor */
348: jpl L1
349: negl d1
350: negb d2 /* change sign because divisor <0 */
351: L1: movel sp@(8), d0 /* d0 = dividend */
352: jpl L2
353: negl d0
354: negb d2
355:
356: L2: movel d1, sp@-
357: movel d0, sp@-
358: jbsr SYM (__udivsi3) /* divide abs(dividend) by abs(divisor) */
359: addql IMM (8), sp
360:
361: tstb d2
362: jpl L3
363: negl d0
364:
365: L3: movel sp@+, d2
366: rts
367: #endif /* L_divsi3 */
368:
369: #ifdef L_umodsi3
370: .text
371: .proc
372: .globl SYM (__umodsi3)
373: SYM (__umodsi3):
374: movel sp@(8), d1 /* d1 = divisor */
375: movel sp@(4), d0 /* d0 = dividend */
376: movel d1, sp@-
377: movel d0, sp@-
378: jbsr SYM (__udivsi3)
379: addql IMM (8), sp
380: movel sp@(8), d1 /* d1 = divisor */
381: movel d1, sp@-
382: movel d0, sp@-
383: jbsr SYM (__mulsi3) /* d0 = (a/b)*b */
384: addql IMM (8), sp
385: movel sp@(4), d1 /* d1 = dividend */
386: subl d0, d1 /* d1 = a - (a/b)*b */
387: movel d1, d0
388: rts
389: #endif /* L_umodsi3 */
390:
391: #ifdef L_modsi3
392: .text
393: .proc
394: .globl SYM (__modsi3)
395: SYM (__modsi3):
396: movel sp@(8), d1 /* d1 = divisor */
397: movel sp@(4), d0 /* d0 = dividend */
398: movel d1, sp@-
399: movel d0, sp@-
400: jbsr SYM (__divsi3)
401: addql IMM (8), sp
402: movel sp@(8), d1 /* d1 = divisor */
403: movel d1, sp@-
404: movel d0, sp@-
405: jbsr SYM (__mulsi3) /* d0 = (a/b)*b */
406: addql IMM (8), sp
407: movel sp@(4), d1 /* d1 = dividend */
408: subl d0, d1 /* d1 = a - (a/b)*b */
409: movel d1, d0
410: rts
411: #endif /* L_modsi3 */
412:
413:
414: #ifdef L_double
415:
416: .globl SYM (_fpCCR)
417: .globl $_exception_handler
418:
419: QUIET_NaN = 0xffffffff
420:
421: D_MAX_EXP = 0x07ff
422: D_BIAS = 1022
423: DBL_MAX_EXP = D_MAX_EXP - D_BIAS
424: DBL_MIN_EXP = 1 - D_BIAS
425: DBL_MANT_DIG = 53
426:
427: INEXACT_RESULT = 0x0001
428: UNDERFLOW = 0x0002
429: OVERFLOW = 0x0004
430: DIVIDE_BY_ZERO = 0x0008
431: INVALID_OPERATION = 0x0010
432:
433: DOUBLE_FLOAT = 2
434:
435: NOOP = 0
436: ADD = 1
437: MULTIPLY = 2
438: DIVIDE = 3
439: NEGATE = 4
440: COMPARE = 5
441: EXTENDSFDF = 6
442: TRUNCDFSF = 7
443:
444: UNKNOWN = -1
445: ROUND_TO_NEAREST = 0 | round result to nearest representable value
446: ROUND_TO_ZERO = 1 | round result towards zero
447: ROUND_TO_PLUS = 2 | round result towards plus infinity
448: ROUND_TO_MINUS = 3 | round result towards minus infinity
449:
450: | Entry points:
451:
452: .globl SYM (__adddf3)
453: .globl SYM (__subdf3)
454: .globl SYM (__muldf3)
455: .globl SYM (__divdf3)
456: .globl SYM (__negdf2)
457: .globl SYM (__cmpdf2)
458:
459: .text
460: .even
461:
462: | These are common routines to return and signal exceptions.
463:
464: Ld$den:
465: | Return and signal a denormalized number
466: orl d7,d0
467: movew IMM (UNDERFLOW),d7
468: orw IMM (INEXACT_RESULT),d7
469: movew IMM (DOUBLE_FLOAT),d6
470: jmp $_exception_handler
471:
472: Ld$infty:
473: Ld$overflow:
474: | Return a properly signed INFINITY and set the exception flags
475: movel IMM (0x7ff00000),d0
476: movel IMM (0),d1
477: orl d7,d0
478: movew IMM (OVERFLOW),d7
479: orw IMM (INEXACT_RESULT),d7
480: movew IMM (DOUBLE_FLOAT),d6
481: jmp $_exception_handler
482:
483: Ld$underflow:
484: | Return 0 and set the exception flags
485: movel IMM (0),d0
486: movel d0,d1
487: movew IMM (UNDERFLOW),d7
488: orw IMM (INEXACT_RESULT),d7
489: movew IMM (DOUBLE_FLOAT),d6
490: jmp $_exception_handler
491:
492: Ld$inop:
493: | Return a quiet NaN and set the exception flags
494: movel IMM (QUIET_NaN),d0
495: movel d0,d1
496: movew IMM (INVALID_OPERATION),d7
497: orw IMM (INEXACT_RESULT),d7
498: movew IMM (DOUBLE_FLOAT),d6
499: jmp $_exception_handler
500:
501: Ld$div$0:
502: | Return a properly signed INFINITY and set the exception flags
503: movel IMM (0x7ff00000),d0
504: movel IMM (0),d1
505: orl d7,d0
506: movew IMM (DIVIDE_BY_ZERO),d7
507: orw IMM (INEXACT_RESULT),d7
508: movew IMM (DOUBLE_FLOAT),d6
509: jmp $_exception_handler
510:
511: |=============================================================================
512: |=============================================================================
513: | double precision routines
514: |=============================================================================
515: |=============================================================================
516:
517: | A double precision floating point number (double) has the format:
518: |
519: | struct _double {
520: | unsigned int sign : 1; /* sign bit */
521: | unsigned int exponent : 11; /* exponent, shifted by 126 */
522: | unsigned int fraction : 52; /* fraction */
523: | } double;
524: |
525: | Thus sizeof(double) = 8 (64 bits).
526: |
527: | All the routines are callable from C programs, and return the result
528: | in the register pair d0-d1. They also preserve all registers except
529: | d0-d1 and a0-a1.
530:
531: |=============================================================================
532: | __subdf3
533: |=============================================================================
534:
535: | double __subdf3(double, double);
536: SYM (__subdf3):
537: bchg IMM (31),sp@(12) | change sign of second operand
538: | and fall through, so we always add
539: |=============================================================================
540: | __adddf3
541: |=============================================================================
542:
543: | double __adddf3(double, double);
544: SYM (__adddf3):
545: link a6,IMM (0) | everything will be done in registers
546: moveml d2-d7,sp@- | save all data registers and a2 (but d0-d1)
547: movel a6@(8),d0 | get first operand
548: movel a6@(12),d1 |
549: movel a6@(16),d2 | get second operand
550: movel a6@(20),d3 |
551:
552: movel d0,d7 | get d0's sign bit in d7 '
553: addl d1,d1 | check and clear sign bit of a, and gain one
554: addxl d0,d0 | bit of extra precision
555: beq Ladddf$b | if zero return second operand
556:
557: movel d2,d6 | save sign in d6
558: addl d3,d3 | get rid of sign bit and gain one bit of
559: addxl d2,d2 | extra precision
560: beq Ladddf$a | if zero return first operand
561:
562: andl IMM (0x80000000),d7 | isolate a's sign bit '
563: swap d6 | and also b's sign bit '
564: andw IMM (0x8000),d6 |
565: orw d6,d7 | and combine them into d7, so that a's sign '
566: | bit is in the high word and b's is in the '
567: | low word, so d6 is free to be used
568: movel d7,a0 | now save d7 into a0, so d7 is free to
569: | be used also
570:
571: | Get the exponents and check for denormalized and/or infinity.
572:
573: movel IMM (0x001fffff),d6 | mask for the fraction
574: movel IMM (0x00200000),d7 | mask to put hidden bit back
575:
576: movel d0,d4 |
577: andl d6,d0 | get fraction in d0
578: notl d6 | make d6 into mask for the exponent
579: andl d6,d4 | get exponent in d4
580: beq Ladddf$a$den | branch if a is denormalized
581: cmpl d6,d4 | check for INFINITY or NaN
582: beq Ladddf$nf |
583: orl d7,d0 | and put hidden bit back
584: Ladddf$1:
585: swap d4 | shift right exponent so that it starts
586: lsrw IMM (5),d4 | in bit 0 and not bit 20
587: | Now we have a's exponent in d4 and fraction in d0-d1 '
588: movel d2,d5 | save b to get exponent
589: andl d6,d5 | get exponent in d5
590: beq Ladddf$b$den | branch if b is denormalized
591: cmpl d6,d5 | check for INFINITY or NaN
592: beq Ladddf$nf
593: notl d6 | make d6 into mask for the fraction again
594: andl d6,d2 | and get fraction in d2
595: orl d7,d2 | and put hidden bit back
596: Ladddf$2:
597: swap d5 | shift right exponent so that it starts
598: lsrw IMM (5),d5 | in bit 0 and not bit 20
599:
600: | Now we have b's exponent in d5 and fraction in d2-d3. '
601:
602: | The situation now is as follows: the signs are combined in a0, the
603: | numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
604: | and d5 (b). To do the rounding correctly we need to keep all the
605: | bits until the end, so we need to use d0-d1-d2-d3 for the first number
606: | and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
607: | exponents in a2-a3.
608:
609: moveml a2-a3,sp@- | save the address registers
610:
611: movel d4,a2 | save the exponents
612: movel d5,a3 |
613:
614: movel IMM (0),d7 | and move the numbers around
615: movel d7,d6 |
616: movel d3,d5 |
617: movel d2,d4 |
618: movel d7,d3 |
619: movel d7,d2 |
620:
621: | Here we shift the numbers until the exponents are the same, and put
622: | the largest exponent in a2.
623: exg d4,a2 | get exponents back
624: exg d5,a3 |
625: cmpw d4,d5 | compare the exponents
626: beq Ladddf$3 | if equal don't shift '
627: bhi 9f | branch if second exponent is higher
628:
629: | Here we have a's exponent larger than b's, so we have to shift b. We do
630: | this by using as counter d2:
631: 1: movew d4,d2 | move largest exponent to d2
632: subw d5,d2 | and subtract second exponent
633: exg d4,a2 | get back the longs we saved
634: exg d5,a3 |
635: | if difference is too large we don't shift (actually, we can just exit) '
636: cmpw IMM (DBL_MANT_DIG+2),d2
637: bge Ladddf$b$small
638: cmpw IMM (32),d2 | if difference >= 32, shift by longs
639: bge 5f
640: 2: cmpw IMM (16),d2 | if difference >= 16, shift by words
641: bge 6f
642: bra 3f | enter dbra loop
643:
644: 4: lsrl IMM (1),d4
645: roxrl IMM (1),d5
646: roxrl IMM (1),d6
647: roxrl IMM (1),d7
648: 3: dbra d2,4b
649: movel IMM (0),d2
650: movel d2,d3
651: bra Ladddf$4
652: 5:
653: movel d6,d7
654: movel d5,d6
655: movel d4,d5
656: movel IMM (0),d4
657: subw IMM (32),d2
658: bra 2b
659: 6:
660: movew d6,d7
661: swap d7
662: movew d5,d6
663: swap d6
664: movew d4,d5
665: swap d5
666: movew IMM (0),d4
667: swap d4
668: subw IMM (16),d2
669: bra 3b
670:
671: 9: exg d4,d5
672: movew d4,d6
673: subw d5,d6 | keep d5 (largest exponent) in d4
674: exg d4,a2
675: exg d5,a3
676: | if difference is too large we don't shift (actually, we can just exit) '
677: cmpw IMM (DBL_MANT_DIG+2),d6
678: bge Ladddf$a$small
679: cmpw IMM (32),d6 | if difference >= 32, shift by longs
680: bge 5f
681: 2: cmpw IMM (16),d6 | if difference >= 16, shift by words
682: bge 6f
683: bra 3f | enter dbra loop
684:
685: 4: lsrl IMM (1),d0
686: roxrl IMM (1),d1
687: roxrl IMM (1),d2
688: roxrl IMM (1),d3
689: 3: dbra d6,4b
690: movel IMM (0),d7
691: movel d7,d6
692: bra Ladddf$4
693: 5:
694: movel d2,d3
695: movel d1,d2
696: movel d0,d1
697: movel IMM (0),d0
698: subw IMM (32),d6
699: bra 2b
700: 6:
701: movew d2,d3
702: swap d3
703: movew d1,d2
704: swap d2
705: movew d0,d1
706: swap d1
707: movew IMM (0),d0
708: swap d0
709: subw IMM (16),d6
710: bra 3b
711: Ladddf$3:
712: exg d4,a2
713: exg d5,a3
714: Ladddf$4:
715: | Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
716: | the signs in a4.
717:
718: | Here we have to decide whether to add or subtract the numbers:
719: exg d7,a0 | get the signs
720: exg d6,a3 | a3 is free to be used
721: movel d7,d6 |
722: movew IMM (0),d7 | get a's sign in d7 '
723: swap d6 |
724: movew IMM (0),d6 | and b's sign in d6 '
725: eorl d7,d6 | compare the signs
726: bmi Lsubdf$0 | if the signs are different we have
727: | to subtract
728: exg d7,a0 | else we add the numbers
729: exg d6,a3 |
730: addl d7,d3 |
731: addxl d6,d2 |
732: addxl d5,d1 |
733: addxl d4,d0 |
734:
735: movel a2,d4 | return exponent to d4
736: movel a0,d7 |
737: andl IMM (0x80000000),d7 | d7 now has the sign
738:
739: moveml sp@+,a2-a3
740:
741: | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
742: | the case of denormalized numbers in the rounding routine itself).
743: | As in the addition (not in the subtraction!) we could have set
744: | one more bit we check this:
745: btst IMM (DBL_MANT_DIG+1),d0
746: beq 1f
747: lsrl IMM (1),d0
748: roxrl IMM (1),d1
749: roxrl IMM (1),d2
750: roxrl IMM (1),d3
751: addw IMM (1),d4
752: 1:
753: lea Ladddf$5,a0 | to return from rounding routine
754: lea SYM (_fpCCR),a1 | check the rounding mode
755: movew a1@(6),d6 | rounding mode in d6
756: beq Lround$to$nearest
757: cmpw IMM (ROUND_TO_PLUS),d6
758: bhi Lround$to$minus
759: blt Lround$to$zero
760: bra Lround$to$plus
761: Ladddf$5:
762: | Put back the exponent and check for overflow
763: cmpw IMM (0x7ff),d4 | is the exponent big?
764: bge 1f
765: bclr IMM (DBL_MANT_DIG-1),d0
766: lslw IMM (4),d4 | put exponent back into position
767: swap d0 |
768: orw d4,d0 |
769: swap d0 |
770: bra Ladddf$ret
771: 1:
772: movew IMM (ADD),d5
773: bra Ld$overflow
774:
775: Lsubdf$0:
776: | Here we do the subtraction.
777: exg d7,a0 | put sign back in a0
778: exg d6,a3 |
779: subl d7,d3 |
780: subxl d6,d2 |
781: subxl d5,d1 |
782: subxl d4,d0 |
783: beq Ladddf$ret$1 | if zero just exit
784: bpl 1f | if positive skip the following
785: exg d7,a0 |
786: bchg IMM (31),d7 | change sign bit in d7
787: exg d7,a0 |
788: negl d3 |
789: negxl d2 |
790: negxl d1 | and negate result
791: negxl d0 |
792: 1:
793: movel a2,d4 | return exponent to d4
794: movel a0,d7
795: andl IMM (0x80000000),d7 | isolate sign bit
796: moveml sp@+,a2-a3 |
797:
798: | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
799: | the case of denormalized numbers in the rounding routine itself).
800: | As in the addition (not in the subtraction!) we could have set
801: | one more bit we check this:
802: btst IMM (DBL_MANT_DIG+1),d0
803: beq 1f
804: lsrl IMM (1),d0
805: roxrl IMM (1),d1
806: roxrl IMM (1),d2
807: roxrl IMM (1),d3
808: addw IMM (1),d4
809: 1:
810: lea Lsubdf$1,a0 | to return from rounding routine
811: lea SYM (_fpCCR),a1 | check the rounding mode
812: movew a1@(6),d6 | rounding mode in d6
813: beq Lround$to$nearest
814: cmpw IMM (ROUND_TO_PLUS),d6
815: bhi Lround$to$minus
816: blt Lround$to$zero
817: bra Lround$to$plus
818: Lsubdf$1:
819: | Put back the exponent and sign (we don't have overflow). '
820: bclr IMM (DBL_MANT_DIG-1),d0
821: lslw IMM (4),d4 | put exponent back into position
822: swap d0 |
823: orw d4,d0 |
824: swap d0 |
825: bra Ladddf$ret
826:
827: | If one of the numbers was too small (difference of exponents >=
828: | DBL_MANT_DIG+1) we return the other (and now we don't have to '
829: | check for finiteness or zero).
830: Ladddf$a$small:
831: moveml sp@+,a2-a3
832: movel a6@(16),d0
833: movel a6@(20),d1
834: lea SYM (_fpCCR),a0
835: movew IMM (0),a0@
836: moveml sp@+,d2-d7 | restore data registers
837: unlk a6 | and return
838: rts
839:
840: Ladddf$b$small:
841: moveml sp@+,a2-a3
842: movel a6@(8),d0
843: movel a6@(12),d1
844: lea SYM (_fpCCR),a0
845: movew IMM (0),a0@
846: moveml sp@+,d2-d7 | restore data registers
847: unlk a6 | and return
848: rts
849:
850: Ladddf$a$den:
851: movel d7,d4 | d7 contains 0x00200000
852: bra Ladddf$1
853:
854: Ladddf$b$den:
855: movel d7,d5 | d7 contains 0x00200000
856: notl d6
857: bra Ladddf$2
858:
859: Ladddf$b:
860: | Return b (if a is zero)
861: movel d2,d0
862: movel d3,d1
863: bra 1f
864: Ladddf$a:
865: movel a6@(8),d0
866: movel a6@(12),d1
867: 1:
868: movew IMM (ADD),d5
869: | Check for NaN and +/-INFINITY.
870: movel d0,d7 |
871: andl IMM (0x80000000),d7 |
872: bclr IMM (31),d0 |
873: cmpl IMM (0x7ff00000),d0 |
874: bge 2f |
875: movel d0,d0 | check for zero, since we don't '
876: bne Ladddf$ret | want to return -0 by mistake
877: bclr IMM (31),d7 |
878: bra Ladddf$ret |
879: 2:
880: andl IMM (0x000fffff),d0 | check for NaN (nonzero fraction)
881: orl d1,d0 |
882: bne Ld$inop |
883: bra Ld$infty |
884:
885: Ladddf$ret$1:
886: moveml sp@+,a2-a3 | restore regs and exit
887:
888: Ladddf$ret:
889: | Normal exit.
890: lea SYM (_fpCCR),a0
891: movew IMM (0),a0@
892: orl d7,d0 | put sign bit back
893: moveml sp@+,d2-d7
894: unlk a6
895: rts
896:
897: Ladddf$ret$den:
898: | Return a denormalized number.
899: lsrl IMM (1),d0 | shift right once more
900: roxrl IMM (1),d1 |
901: bra Ladddf$ret
902:
903: Ladddf$nf:
904: movew IMM (ADD),d5
905: | This could be faster but it is not worth the effort, since it is not
906: | executed very often. We sacrifice speed for clarity here.
907: movel a6@(8),d0 | get the numbers back (remember that we
908: movel a6@(12),d1 | did some processing already)
909: movel a6@(16),d2 |
910: movel a6@(20),d3 |
911: movel IMM (0x7ff00000),d4 | useful constant (INFINITY)
912: movel d0,d7 | save sign bits
913: movel d2,d6 |
914: bclr IMM (31),d0 | clear sign bits
915: bclr IMM (31),d2 |
916: | We know that one of them is either NaN of +/-INFINITY
917: | Check for NaN (if either one is NaN return NaN)
918: cmpl d4,d0 | check first a (d0)
919: bhi Ld$inop | if d0 > 0x7ff00000 or equal and
920: bne 2f
921: tstl d1 | d1 > 0, a is NaN
922: bne Ld$inop |
923: 2: cmpl d4,d2 | check now b (d1)
924: bhi Ld$inop |
925: bne 3f
926: tstl d3 |
927: bne Ld$inop |
928: 3:
929: | Now comes the check for +/-INFINITY. We know that both are (maybe not
930: | finite) numbers, but we have to check if both are infinite whether we
931: | are adding or subtracting them.
932: eorl d7,d6 | to check sign bits
933: bmi 1f
934: andl IMM (0x80000000),d7 | get (common) sign bit
935: bra Ld$infty
936: 1:
937: | We know one (or both) are infinite, so we test for equality between the
938: | two numbers (if they are equal they have to be infinite both, so we
939: | return NaN).
940: cmpl d2,d0 | are both infinite?
941: bne 1f | if d0 <> d2 they are not equal
942: cmpl d3,d1 | if d0 == d2 test d3 and d1
943: beq Ld$inop | if equal return NaN
944: 1:
945: andl IMM (0x80000000),d7 | get a's sign bit '
946: cmpl d4,d0 | test now for infinity
947: beq Ld$infty | if a is INFINITY return with this sign
948: bchg IMM (31),d7 | else we know b is INFINITY and has
949: bra Ld$infty | the opposite sign
950:
951: |=============================================================================
952: | __muldf3
953: |=============================================================================
954:
955: | double __muldf3(double, double);
956: SYM (__muldf3):
957: link a6,IMM (0)
958: moveml d2-d7,sp@-
959: movel a6@(8),d0 | get a into d0-d1
960: movel a6@(12),d1 |
961: movel a6@(16),d2 | and b into d2-d3
962: movel a6@(20),d3 |
963: movel d0,d7 | d7 will hold the sign of the product
964: eorl d2,d7 |
965: andl IMM (0x80000000),d7 |
966: movel d7,a0 | save sign bit into a0
967: movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
968: movel d7,d6 | another (mask for fraction)
969: notl d6 |
970: bclr IMM (31),d0 | get rid of a's sign bit '
971: movel d0,d4 |
972: orl d1,d4 |
973: beq Lmuldf$a$0 | branch if a is zero
974: movel d0,d4 |
975: bclr IMM (31),d2 | get rid of b's sign bit '
976: movel d2,d5 |
977: orl d3,d5 |
978: beq Lmuldf$b$0 | branch if b is zero
979: movel d2,d5 |
980: cmpl d7,d0 | is a big?
981: bhi Lmuldf$inop | if a is NaN return NaN
982: beq Lmuldf$a$nf | we still have to check d1 and b ...
983: cmpl d7,d2 | now compare b with INFINITY
984: bhi Lmuldf$inop | is b NaN?
985: beq Lmuldf$b$nf | we still have to check d3 ...
986: | Here we have both numbers finite and nonzero (and with no sign bit).
987: | Now we get the exponents into d4 and d5.
988: andl d7,d4 | isolate exponent in d4
989: beq Lmuldf$a$den | if exponent zero, have denormalized
990: andl d6,d0 | isolate fraction
991: orl IMM (0x00100000),d0 | and put hidden bit back
992: swap d4 | I like exponents in the first byte
993: lsrw IMM (4),d4 |
994: Lmuldf$1:
995: andl d7,d5 |
996: beq Lmuldf$b$den |
997: andl d6,d2 |
998: orl IMM (0x00100000),d2 | and put hidden bit back
999: swap d5 |
1000: lsrw IMM (4),d5 |
1001: Lmuldf$2: |
1002: addw d5,d4 | add exponents
1003: subw IMM (D_BIAS+1),d4 | and subtract bias (plus one)
1004:
1005: | We are now ready to do the multiplication. The situation is as follows:
1006: | both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were
1007: | denormalized to start with!), which means that in the product bit 104
1008: | (which will correspond to bit 8 of the fourth long) is set.
1009:
1010: | Here we have to do the product.
1011: | To do it we have to juggle the registers back and forth, as there are not
1012: | enough to keep everything in them. So we use the address registers to keep
1013: | some intermediate data.
1014:
1015: moveml a2-a3,sp@- | save a2 and a3 for temporary use
1016: movel IMM (0),a2 | a2 is a null register
1017: movel d4,a3 | and a3 will preserve the exponent
1018:
1019: | First, shift d2-d3 so bit 20 becomes bit 31:
1020: rorl IMM (5),d2 | rotate d2 5 places right
1021: swap d2 | and swap it
1022: rorl IMM (5),d3 | do the same thing with d3
1023: swap d3 |
1024: movew d3,d6 | get the rightmost 11 bits of d3
1025: andw IMM (0x07ff),d6 |
1026: orw d6,d2 | and put them into d2
1027: andw IMM (0xf800),d3 | clear those bits in d3
1028:
1029: movel d2,d6 | move b into d6-d7
1030: movel d3,d7 | move a into d4-d5
1031: movel d0,d4 | and clear d0-d1-d2-d3 (to put result)
1032: movel d1,d5 |
1033: movel IMM (0),d3 |
1034: movel d3,d2 |
1035: movel d3,d1 |
1036: movel d3,d0 |
1037:
1038: | We use a1 as counter:
1039: movel IMM (DBL_MANT_DIG-1),a1
1040: exg d7,a1
1041:
1042: 1: exg d7,a1 | put counter back in a1
1043: addl d3,d3 | shift sum once left
1044: addxl d2,d2 |
1045: addxl d1,d1 |
1046: addxl d0,d0 |
1047: addl d7,d7 |
1048: addxl d6,d6 |
1049: bcc 2f | if bit clear skip the following
1050: exg d7,a2 |
1051: addl d5,d3 | else add a to the sum
1052: addxl d4,d2 |
1053: addxl d7,d1 |
1054: addxl d7,d0 |
1055: exg d7,a2 |
1056: 2: exg d7,a1 | put counter in d7
1057: dbf d7,1b | decrement and branch
1058:
1059: movel a3,d4 | restore exponent
1060: moveml sp@+,a2-a3
1061:
1062: | Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The
1063: | first thing to do now is to normalize it so bit 8 becomes bit
1064: | DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
1065: swap d0
1066: swap d1
1067: movew d1,d0
1068: swap d2
1069: movew d2,d1
1070: swap d3
1071: movew d3,d2
1072: movew IMM (0),d3
1073: lsrl IMM (1),d0
1074: roxrl IMM (1),d1
1075: roxrl IMM (1),d2
1076: roxrl IMM (1),d3
1077: lsrl IMM (1),d0
1078: roxrl IMM (1),d1
1079: roxrl IMM (1),d2
1080: roxrl IMM (1),d3
1081: lsrl IMM (1),d0
1082: roxrl IMM (1),d1
1083: roxrl IMM (1),d2
1084: roxrl IMM (1),d3
1085:
1086: | Now round, check for over- and underflow, and exit.
1087: movel a0,d7 | get sign bit back into d7
1088: movew IMM (MULTIPLY),d5
1089:
1090: btst IMM (DBL_MANT_DIG+1-32),d0
1091: beq Lround$exit
1092: lsrl IMM (1),d0
1093: roxrl IMM (1),d1
1094: addw IMM (1),d4
1095: bra Lround$exit
1096:
1097: Lmuldf$inop:
1098: movew IMM (MULTIPLY),d5
1099: bra Ld$inop
1100:
1101: Lmuldf$b$nf:
1102: movew IMM (MULTIPLY),d5
1103: movel a0,d7 | get sign bit back into d7
1104: tstl d3 | we know d2 == 0x7ff00000, so check d3
1105: bne Ld$inop | if d3 <> 0 b is NaN
1106: bra Ld$overflow | else we have overflow (since a is finite)
1107:
1108: Lmuldf$a$nf:
1109: movew IMM (MULTIPLY),d5
1110: movel a0,d7 | get sign bit back into d7
1111: tstl d1 | we know d0 == 0x7ff00000, so check d1
1112: bne Ld$inop | if d1 <> 0 a is NaN
1113: bra Ld$overflow | else signal overflow
1114:
1115: | If either number is zero return zero, unless the other is +/-INFINITY or
1116: | NaN, in which case we return NaN.
1117: Lmuldf$b$0:
1118: movew IMM (MULTIPLY),d5
1119: exg d2,d0 | put b (==0) into d0-d1
1120: exg d3,d1 | and a (with sign bit cleared) into d2-d3
1121: bra 1f
1122: Lmuldf$a$0:
1123: movel a6@(16),d2 | put b into d2-d3 again
1124: movel a6@(20),d3 |
1125: bclr IMM (31),d2 | clear sign bit
1126: 1: cmpl IMM (0x7ff00000),d2 | check for non-finiteness
1127: bge Ld$inop | in case NaN or +/-INFINITY return NaN
1128: lea SYM (_fpCCR),a0
1129: movew IMM (0),a0@
1130: moveml sp@+,d2-d7
1131: unlk a6
1132: rts
1133:
1134: | If a number is denormalized we put an exponent of 1 but do not put the
1135: | hidden bit back into the fraction; instead we shift left until bit 21
1136: | (the hidden bit) is set, adjusting the exponent accordingly. We do this
1137: | to ensure that the product of the fractions is close to 1.
1138: Lmuldf$a$den:
1139: movel IMM (1),d4
1140: andl d6,d0
1141: 1: addl d1,d1 | shift a left until bit 20 is set
1142: addxl d0,d0 |
1143: subw IMM (1),d4 | and adjust exponent
1144: btst IMM (20),d0 |
1145: bne Lmuldf$1 |
1146: bra 1b
1147:
1148: Lmuldf$b$den:
1149: movel IMM (1),d5
1150: andl d6,d2
1151: 1: addl d3,d3 | shift b left until bit 20 is set
1152: addxl d2,d2 |
1153: subw IMM (1),d5 | and adjust exponent
1154: btst IMM (20),d2 |
1155: bne Lmuldf$2 |
1156: bra 1b
1157:
1158:
1159: |=============================================================================
1160: | __divdf3
1161: |=============================================================================
1162:
1163: | double __divdf3(double, double);
1164: SYM (__divdf3):
1165: link a6,IMM (0)
1166: moveml d2-d7,sp@-
1167: movel a6@(8),d0 | get a into d0-d1
1168: movel a6@(12),d1 |
1169: movel a6@(16),d2 | and b into d2-d3
1170: movel a6@(20),d3 |
1171: movel d0,d7 | d7 will hold the sign of the result
1172: eorl d2,d7 |
1173: andl IMM (0x80000000),d7
1174: movel d7,a0 | save sign into a0
1175: movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
1176: movel d7,d6 | another (mask for fraction)
1177: notl d6 |
1178: bclr IMM (31),d0 | get rid of a's sign bit '
1179: movel d0,d4 |
1180: orl d1,d4 |
1181: beq Ldivdf$a$0 | branch if a is zero
1182: movel d0,d4 |
1183: bclr IMM (31),d2 | get rid of b's sign bit '
1184: movel d2,d5 |
1185: orl d3,d5 |
1186: beq Ldivdf$b$0 | branch if b is zero
1187: movel d2,d5
1188: cmpl d7,d0 | is a big?
1189: bhi Ldivdf$inop | if a is NaN return NaN
1190: beq Ldivdf$a$nf | if d0 == 0x7ff00000 we check d1
1191: cmpl d7,d2 | now compare b with INFINITY
1192: bhi Ldivdf$inop | if b is NaN return NaN
1193: beq Ldivdf$b$nf | if d2 == 0x7ff00000 we check d3
1194: | Here we have both numbers finite and nonzero (and with no sign bit).
1195: | Now we get the exponents into d4 and d5 and normalize the numbers to
1196: | ensure that the ratio of the fractions is around 1. We do this by
1197: | making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
1198: | set, even if they were denormalized to start with.
1199: | Thus, the result will satisfy: 2 > result > 1/2.
1200: andl d7,d4 | and isolate exponent in d4
1201: beq Ldivdf$a$den | if exponent is zero we have a denormalized
1202: andl d6,d0 | and isolate fraction
1203: orl IMM (0x00100000),d0 | and put hidden bit back
1204: swap d4 | I like exponents in the first byte
1205: lsrw IMM (4),d4 |
1206: Ldivdf$1: |
1207: andl d7,d5 |
1208: beq Ldivdf$b$den |
1209: andl d6,d2 |
1210: orl IMM (0x00100000),d2
1211: swap d5 |
1212: lsrw IMM (4),d5 |
1213: Ldivdf$2: |
1214: subw d5,d4 | subtract exponents
1215: addw IMM (D_BIAS),d4 | and add bias
1216:
1217: | We are now ready to do the division. We have prepared things in such a way
1218: | that the ratio of the fractions will be less than 2 but greater than 1/2.
1219: | At this point the registers in use are:
1220: | d0-d1 hold a (first operand, bit DBL_MANT_DIG-32=0, bit
1221: | DBL_MANT_DIG-1-32=1)
1222: | d2-d3 hold b (second operand, bit DBL_MANT_DIG-32=1)
1223: | d4 holds the difference of the exponents, corrected by the bias
1224: | a0 holds the sign of the ratio
1225:
1226: | To do the rounding correctly we need to keep information about the
1227: | nonsignificant bits. One way to do this would be to do the division
1228: | using four registers; another is to use two registers (as originally
1229: | I did), but use a sticky bit to preserve information about the
1230: | fractional part. Note that we can keep that info in a1, which is not
1231: | used.
1232: movel IMM (0),d6 | d6-d7 will hold the result
1233: movel d6,d7 |
1234: movel IMM (0),a1 | and a1 will hold the sticky bit
1235:
1236: movel IMM (DBL_MANT_DIG-32+1),d5
1237:
1238: 1: cmpl d0,d2 | is a < b?
1239: bhi 3f | if b > a skip the following
1240: beq 4f | if d0==d2 check d1 and d3
1241: 2: subl d3,d1 |
1242: subxl d2,d0 | a <-- a - b
1243: bset d5,d6 | set the corresponding bit in d6
1244: 3: addl d1,d1 | shift a by 1
1245: addxl d0,d0 |
1246: dbra d5,1b | and branch back
1247: bra 5f
1248: 4: cmpl d1,d3 | here d0==d2, so check d1 and d3
1249: bhi 3b | if d1 > d2 skip the subtraction
1250: bra 2b | else go do it
1251: 5:
1252: | Here we have to start setting the bits in the second long.
1253: movel IMM (31),d5 | again d5 is counter
1254:
1255: 1: cmpl d0,d2 | is a < b?
1256: bhi 3f | if b > a skip the following
1257: beq 4f | if d0==d2 check d1 and d3
1258: 2: subl d3,d1 |
1259: subxl d2,d0 | a <-- a - b
1260: bset d5,d7 | set the corresponding bit in d7
1261: 3: addl d1,d1 | shift a by 1
1262: addxl d0,d0 |
1263: dbra d5,1b | and branch back
1264: bra 5f
1265: 4: cmpl d1,d3 | here d0==d2, so check d1 and d3
1266: bhi 3b | if d1 > d2 skip the subtraction
1267: bra 2b | else go do it
1268: 5:
1269: | Now go ahead checking until we hit a one, which we store in d2.
1270: movel IMM (DBL_MANT_DIG),d5
1271: 1: cmpl d2,d0 | is a < b?
1272: bhi 4f | if b < a, exit
1273: beq 3f | if d0==d2 check d1 and d3
1274: 2: addl d1,d1 | shift a by 1
1275: addxl d0,d0 |
1276: dbra d5,1b | and branch back
1277: movel IMM (0),d2 | here no sticky bit was found
1278: movel d2,d3
1279: bra 5f
1280: 3: cmpl d1,d3 | here d0==d2, so check d1 and d3
1281: bhi 2b | if d1 > d2 go back
1282: 4:
1283: | Here put the sticky bit in d2-d3 (in the position which actually corresponds
1284: | to it; if you don't do this the algorithm loses in some cases). '
1285: movel IMM (0),d2
1286: movel d2,d3
1287: subw IMM (DBL_MANT_DIG),d5
1288: addw IMM (63),d5
1289: cmpw IMM (31),d5
1290: bhi 2f
1291: 1: bset d5,d3
1292: bra 5f
1293: subw IMM (32),d5
1294: 2: bset d5,d2
1295: 5:
1296: | Finally we are finished! Move the longs in the address registers to
1297: | their final destination:
1298: movel d6,d0
1299: movel d7,d1
1300: movel IMM (0),d3
1301:
1302: | Here we have finished the division, with the result in d0-d1-d2-d3, with
1303: | 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
1304: | If it is not, then definitely bit 21 is set. Normalize so bit 22 is
1305: | not set:
1306: btst IMM (DBL_MANT_DIG-32+1),d0
1307: beq 1f
1308: lsrl IMM (1),d0
1309: roxrl IMM (1),d1
1310: roxrl IMM (1),d2
1311: roxrl IMM (1),d3
1312: addw IMM (1),d4
1313: 1:
1314: | Now round, check for over- and underflow, and exit.
1315: movel a0,d7 | restore sign bit to d7
1316: movew IMM (DIVIDE),d5
1317: bra Lround$exit
1318:
1319: Ldivdf$inop:
1320: movew IMM (DIVIDE),d5
1321: bra Ld$inop
1322:
1323: Ldivdf$a$0:
1324: | If a is zero check to see whether b is zero also. In that case return
1325: | NaN; then check if b is NaN, and return NaN also in that case. Else
1326: | return zero.
1327: movew IMM (DIVIDE),d5
1328: bclr IMM (31),d2 |
1329: movel d2,d4 |
1330: orl d3,d4 |
1331: beq Ld$inop | if b is also zero return NaN
1332: cmpl IMM (0x7ff00000),d2 | check for NaN
1333: bhi Ld$inop |
1334: blt 1f |
1335: tstl d3 |
1336: bne Ld$inop |
1337: 1: movel IMM (0),d0 | else return zero
1338: movel d0,d1 |
1339: lea SYM (_fpCCR),a0 | clear exception flags
1340: movew IMM (0),a0@ |
1341: moveml sp@+,d2-d7 |
1342: unlk a6 |
1343: rts |
1344:
1345: Ldivdf$b$0:
1346: movew IMM (DIVIDE),d5
1347: | If we got here a is not zero. Check if a is NaN; in that case return NaN,
1348: | else return +/-INFINITY. Remember that a is in d0 with the sign bit
1349: | cleared already.
1350: movel a0,d7 | put a's sign bit back in d7 '
1351: cmpl IMM (0x7ff00000),d0 | compare d0 with INFINITY
1352: bhi Ld$inop | if larger it is NaN
1353: tstl d1 |
1354: bne Ld$inop |
1355: bra Ld$div$0 | else signal DIVIDE_BY_ZERO
1356:
1357: Ldivdf$b$nf:
1358: movew IMM (DIVIDE),d5
1359: | If d2 == 0x7ff00000 we have to check d3.
1360: tstl d3 |
1361: bne Ld$inop | if d3 <> 0, b is NaN
1362: bra Ld$underflow | else b is +/-INFINITY, so signal underflow
1363:
1364: Ldivdf$a$nf:
1365: movew IMM (DIVIDE),d5
1366: | If d0 == 0x7ff00000 we have to check d1.
1367: tstl d1 |
1368: bne Ld$inop | if d1 <> 0, a is NaN
1369: | If a is INFINITY we have to check b
1370: cmpl d7,d2 | compare b with INFINITY
1371: bge Ld$inop | if b is NaN or INFINITY return NaN
1372: tstl d3 |
1373: bne Ld$inop |
1374: bra Ld$overflow | else return overflow
1375:
1376: | If a number is denormalized we put an exponent of 1 but do not put the
1377: | bit back into the fraction.
1378: Ldivdf$a$den:
1379: movel IMM (1),d4
1380: andl d6,d0
1381: 1: addl d1,d1 | shift a left until bit 20 is set
1382: addxl d0,d0
1383: subw IMM (1),d4 | and adjust exponent
1384: btst IMM (DBL_MANT_DIG-32-1),d0
1385: bne Ldivdf$1
1386: bra 1b
1387:
1388: Ldivdf$b$den:
1389: movel IMM (1),d5
1390: andl d6,d2
1391: 1: addl d3,d3 | shift b left until bit 20 is set
1392: addxl d2,d2
1393: subw IMM (1),d5 | and adjust exponent
1394: btst IMM (DBL_MANT_DIG-32-1),d2
1395: bne Ldivdf$2
1396: bra 1b
1397:
1398: Lround$exit:
1399: | This is a common exit point for __muldf3 and __divdf3. When they enter
1400: | this point the sign of the result is in d7, the result in d0-d1, normalized
1401: | so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.
1402:
1403: | First check for underlow in the exponent:
1404: cmpw IMM (-DBL_MANT_DIG-1),d4
1405: blt Ld$underflow
1406: | It could happen that the exponent is less than 1, in which case the
1407: | number is denormalized. In this case we shift right and adjust the
1408: | exponent until it becomes 1 or the fraction is zero (in the latter case
1409: | we signal underflow and return zero).
1410: movel d7,a0 |
1411: movel IMM (0),d6 | use d6-d7 to collect bits flushed right
1412: movel d6,d7 | use d6-d7 to collect bits flushed right
1413: cmpw IMM (1),d4 | if the exponent is less than 1 we
1414: bge 2f | have to shift right (denormalize)
1415: 1: addw IMM (1),d4 | adjust the exponent
1416: lsrl IMM (1),d0 | shift right once
1417: roxrl IMM (1),d1 |
1418: roxrl IMM (1),d2 |
1419: roxrl IMM (1),d3 |
1420: roxrl IMM (1),d6 |
1421: roxrl IMM (1),d7 |
1422: cmpw IMM (1),d4 | is the exponent 1 already?
1423: beq 2f | if not loop back
1424: bra 1b |
1425: bra Ld$underflow | safety check, shouldn't execute '
1426: 2: orl d6,d2 | this is a trick so we don't lose '
1427: orl d7,d3 | the bits which were flushed right
1428: movel a0,d7 | get back sign bit into d7
1429: | Now call the rounding routine (which takes care of denormalized numbers):
1430: lea Lround$0,a0 | to return from rounding routine
1431: lea SYM (_fpCCR),a1 | check the rounding mode
1432: movew a1@(6),d6 | rounding mode in d6
1433: beq Lround$to$nearest
1434: cmpw IMM (ROUND_TO_PLUS),d6
1435: bhi Lround$to$minus
1436: blt Lround$to$zero
1437: bra Lround$to$plus
1438: Lround$0:
1439: | Here we have a correctly rounded result (either normalized or denormalized).
1440:
1441: | Here we should have either a normalized number or a denormalized one, and
1442: | the exponent is necessarily larger or equal to 1 (so we don't have to '
1443: | check again for underflow!). We have to check for overflow or for a
1444: | denormalized number (which also signals underflow).
1445: | Check for overflow (i.e., exponent >= 0x7ff).
1446: cmpw IMM (0x07ff),d4
1447: bge Ld$overflow
1448: | Now check for a denormalized number (exponent==0):
1449: movew d4,d4
1450: beq Ld$den
1451: 1:
1452: | Put back the exponents and sign and return.
1453: lslw IMM (4),d4 | exponent back to fourth byte
1454: bclr IMM (DBL_MANT_DIG-32-1),d0
1455: swap d0 | and put back exponent
1456: orw d4,d0 |
1457: swap d0 |
1458: orl d7,d0 | and sign also
1459:
1460: lea SYM (_fpCCR),a0
1461: movew IMM (0),a0@
1462: moveml sp@+,d2-d7
1463: unlk a6
1464: rts
1465:
1466: |=============================================================================
1467: | __negdf2
1468: |=============================================================================
1469:
1470: | double __negdf2(double, double);
1471: SYM (__negdf2):
1472: link a6,IMM (0)
1473: moveml d2-d7,sp@-
1474: movew IMM (NEGATE),d5
1475: movel a6@(8),d0 | get number to negate in d0-d1
1476: movel a6@(12),d1 |
1477: bchg IMM (31),d0 | negate
1478: movel d0,d2 | make a positive copy (for the tests)
1479: bclr IMM (31),d2 |
1480: movel d2,d4 | check for zero
1481: orl d1,d4 |
1482: beq 2f | if zero (either sign) return +zero
1483: cmpl IMM (0x7ff00000),d2 | compare to +INFINITY
1484: blt 1f | if finite, return
1485: bhi Ld$inop | if larger (fraction not zero) is NaN
1486: tstl d1 | if d2 == 0x7ff00000 check d1
1487: bne Ld$inop |
1488: movel d0,d7 | else get sign and return INFINITY
1489: andl IMM (0x80000000),d7
1490: bra Ld$infty
1491: 1: lea SYM (_fpCCR),a0
1492: movew IMM (0),a0@
1493: moveml sp@+,d2-d7
1494: unlk a6
1495: rts
1496: 2: bclr IMM (31),d0
1497: bra 1b
1498:
1499: |=============================================================================
1500: | __cmpdf2
1501: |=============================================================================
1502:
1503: GREATER = 1
1504: LESS = -1
1505: EQUAL = 0
1506:
1507: | int __cmpdf2(double, double);
1508: SYM (__cmpdf2):
1509: link a6,IMM (0)
1510: moveml d2-d7,sp@- | save registers
1511: movew IMM (COMPARE),d5
1512: movel a6@(8),d0 | get first operand
1513: movel a6@(12),d1 |
1514: movel a6@(16),d2 | get second operand
1515: movel a6@(20),d3 |
1516: | First check if a and/or b are (+/-) zero and in that case clear
1517: | the sign bit.
1518: movel d0,d6 | copy signs into d6 (a) and d7(b)
1519: bclr IMM (31),d0 | and clear signs in d0 and d2
1520: movel d2,d7 |
1521: bclr IMM (31),d2 |
1522: cmpl IMM (0x7fff0000),d0 | check for a == NaN
1523: bhi Ld$inop | if d0 > 0x7ff00000, a is NaN
1524: beq Lcmpdf$a$nf | if equal can be INFINITY, so check d1
1525: movel d0,d4 | copy into d4 to test for zero
1526: orl d1,d4 |
1527: beq Lcmpdf$a$0 |
1528: Lcmpdf$0:
1529: cmpl IMM (0x7fff0000),d2 | check for b == NaN
1530: bhi Ld$inop | if d2 > 0x7ff00000, b is NaN
1531: beq Lcmpdf$b$nf | if equal can be INFINITY, so check d3
1532: movel d2,d4 |
1533: orl d3,d4 |
1534: beq Lcmpdf$b$0 |
1535: Lcmpdf$1:
1536: | Check the signs
1537: eorl d6,d7
1538: bpl 1f
1539: | If the signs are not equal check if a >= 0
1540: tstl d6
1541: bpl Lcmpdf$a$gt$b | if (a >= 0 && b < 0) => a > b
1542: bmi Lcmpdf$b$gt$a | if (a < 0 && b >= 0) => a < b
1543: 1:
1544: | If the signs are equal check for < 0
1545: tstl d6
1546: bpl 1f
1547: | If both are negative exchange them
1548: exg d0,d2
1549: exg d1,d3
1550: 1:
1551: | Now that they are positive we just compare them as longs (does this also
1552: | work for denormalized numbers?).
1553: cmpl d0,d2
1554: bhi Lcmpdf$b$gt$a | |b| > |a|
1555: bne Lcmpdf$a$gt$b | |b| < |a|
1556: | If we got here d0 == d2, so we compare d1 and d3.
1557: cmpl d1,d3
1558: bhi Lcmpdf$b$gt$a | |b| > |a|
1559: bne Lcmpdf$a$gt$b | |b| < |a|
1560: | If we got here a == b.
1561: movel IMM (EQUAL),d0
1562: moveml sp@+,d2-d7 | put back the registers
1563: unlk a6
1564: rts
1565: Lcmpdf$a$gt$b:
1566: movel IMM (GREATER),d0
1567: moveml sp@+,d2-d7 | put back the registers
1568: unlk a6
1569: rts
1570: Lcmpdf$b$gt$a:
1571: movel IMM (LESS),d0
1572: moveml sp@+,d2-d7 | put back the registers
1573: unlk a6
1574: rts
1575:
1576: Lcmpdf$a$0:
1577: bclr IMM (31),d6
1578: bra Lcmpdf$0
1579: Lcmpdf$b$0:
1580: bclr IMM (31),d7
1581: bra Lcmpdf$1
1582:
1583: Lcmpdf$a$nf:
1584: tstl d1
1585: bne Ld$inop
1586: bra Lcmpdf$0
1587:
1588: Lcmpdf$b$nf:
1589: tstl d3
1590: bne Ld$inop
1591: bra Lcmpdf$1
1592:
1593: |=============================================================================
1594: | rounding routines
1595: |=============================================================================
1596:
1597: | The rounding routines expect the number to be normalized in registers
1598: | d0-d1-d2-d3, with the exponent in register d4. They assume that the
1599: | exponent is larger or equal to 1. They return a properly normalized number
1600: | if possible, and a denormalized number otherwise. The exponent is returned
1601: | in d4.
1602:
1603: Lround$to$nearest:
1604: | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
1605: | Here we assume that the exponent is not too small (this should be checked
1606: | before entering the rounding routine), but the number could be denormalized.
1607:
1608: | Check for denormalized numbers:
1609: 1: btst IMM (DBL_MANT_DIG-32),d0
1610: bne 2f | if set the number is normalized
1611: | Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent
1612: | is one (remember that a denormalized number corresponds to an
1613: | exponent of -D_BIAS+1).
1614: cmpw IMM (1),d4 | remember that the exponent is at least one
1615: beq 2f | an exponent of one means denormalized
1616: addl d3,d3 | else shift and adjust the exponent
1617: addxl d2,d2 |
1618: addxl d1,d1 |
1619: addxl d0,d0 |
1620: dbra d4,1b |
1621: 2:
1622: | Now round: we do it as follows: after the shifting we can write the
1623: | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
1624: | If delta < 1, do nothing. If delta > 1, add 1 to f.
1625: | If delta == 1, we make sure the rounded number will be even (odd?)
1626: | (after shifting).
1627: btst IMM (0),d1 | is delta < 1?
1628: beq 2f | if so, do not do anything
1629: orl d2,d3 | is delta == 1?
1630: bne 1f | if so round to even
1631: movel d1,d3 |
1632: andl IMM (2),d3 | bit 1 is the last significant bit
1633: movel IMM (0),d2 |
1634: addl d3,d1 |
1635: addxl d2,d0 |
1636: bra 2f |
1637: 1: movel IMM (1),d3 | else add 1
1638: movel IMM (0),d2 |
1639: addl d3,d1 |
1640: addxl d2,d0
1641: | Shift right once (because we used bit #DBL_MANT_DIG-32!).
1642: 2: lsrl IMM (1),d0
1643: roxrl IMM (1),d1
1644:
1645: | Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
1646: | 'fraction overflow' ...).
1647: btst IMM (DBL_MANT_DIG-32),d0
1648: beq 1f
1649: lsrl IMM (1),d0
1650: roxrl IMM (1),d1
1651: addw IMM (1),d4
1652: 1:
1653: | If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we
1654: | have to put the exponent to zero and return a denormalized number.
1655: btst IMM (DBL_MANT_DIG-32-1),d0
1656: beq 1f
1657: jmp a0@
1658: 1: movel IMM (0),d4
1659: jmp a0@
1660:
1661: Lround$to$zero:
1662: Lround$to$plus:
1663: Lround$to$minus:
1664: jmp a0@
1665: #endif /* L_double */
1666:
1667: #ifdef L_float
1668:
1669: .globl SYM (_fpCCR)
1670: .globl $_exception_handler
1671:
1672: QUIET_NaN = 0xffffffff
1673: SIGNL_NaN = 0x7f800001
1674: INFINITY = 0x7f800000
1675:
1676: F_MAX_EXP = 0xff
1677: F_BIAS = 126
1678: FLT_MAX_EXP = F_MAX_EXP - F_BIAS
1679: FLT_MIN_EXP = 1 - F_BIAS
1680: FLT_MANT_DIG = 24
1681:
1682: INEXACT_RESULT = 0x0001
1683: UNDERFLOW = 0x0002
1684: OVERFLOW = 0x0004
1685: DIVIDE_BY_ZERO = 0x0008
1686: INVALID_OPERATION = 0x0010
1687:
1688: SINGLE_FLOAT = 1
1689:
1690: NOOP = 0
1691: ADD = 1
1692: MULTIPLY = 2
1693: DIVIDE = 3
1694: NEGATE = 4
1695: COMPARE = 5
1696: EXTENDSFDF = 6
1697: TRUNCDFSF = 7
1698:
1699: UNKNOWN = -1
1700: ROUND_TO_NEAREST = 0 | round result to nearest representable value
1701: ROUND_TO_ZERO = 1 | round result towards zero
1702: ROUND_TO_PLUS = 2 | round result towards plus infinity
1703: ROUND_TO_MINUS = 3 | round result towards minus infinity
1704:
1705: | Entry points:
1706:
1707: .globl SYM (__addsf3)
1708: .globl SYM (__subsf3)
1709: .globl SYM (__mulsf3)
1710: .globl SYM (__divsf3)
1711: .globl SYM (__negsf2)
1712: .globl SYM (__cmpsf2)
1713:
1714: | These are common routines to return and signal exceptions.
1715:
1716: .text
1717: .even
1718:
1719: Lf$den:
1720: | Return and signal a denormalized number
1721: orl d7,d0
1722: movew IMM (UNDERFLOW),d7
1723: orw IMM (INEXACT_RESULT),d7
1724: movew IMM (SINGLE_FLOAT),d6
1725: jmp $_exception_handler
1726:
1727: Lf$infty:
1728: Lf$overflow:
1729: | Return a properly signed INFINITY and set the exception flags
1730: movel IMM (INFINITY),d0
1731: orl d7,d0
1732: movew IMM (OVERFLOW),d7
1733: orw IMM (INEXACT_RESULT),d7
1734: movew IMM (SINGLE_FLOAT),d6
1735: jmp $_exception_handler
1736:
1737: Lf$underflow:
1738: | Return 0 and set the exception flags
1739: movel IMM (0),d0
1740: movew IMM (UNDERFLOW),d7
1741: orw IMM (INEXACT_RESULT),d7
1742: movew IMM (SINGLE_FLOAT),d6
1743: jmp $_exception_handler
1744:
1745: Lf$inop:
1746: | Return a quiet NaN and set the exception flags
1747: movel IMM (QUIET_NaN),d0
1748: movew IMM (INVALID_OPERATION),d7
1749: orw IMM (INEXACT_RESULT),d7
1750: movew IMM (SINGLE_FLOAT),d6
1751: jmp $_exception_handler
1752:
1753: Lf$div$0:
1754: | Return a properly signed INFINITY and set the exception flags
1755: movel IMM (INFINITY),d0
1756: orl d7,d0
1757: movew IMM (DIVIDE_BY_ZERO),d7
1758: orw IMM (INEXACT_RESULT),d7
1759: movew IMM (SINGLE_FLOAT),d6
1760: jmp $_exception_handler
1761:
1762: |=============================================================================
1763: |=============================================================================
1764: | single precision routines
1765: |=============================================================================
1766: |=============================================================================
1767:
1768: | A single precision floating point number (float) has the format:
1769: |
1770: | struct _float {
1771: | unsigned int sign : 1; /* sign bit */
1772: | unsigned int exponent : 8; /* exponent, shifted by 126 */
1773: | unsigned int fraction : 23; /* fraction */
1774: | } float;
1775: |
1776: | Thus sizeof(float) = 4 (32 bits).
1777: |
1778: | All the routines are callable from C programs, and return the result
1779: | in the single register d0. They also preserve all registers except
1780: | d0-d1 and a0-a1.
1781:
1782: |=============================================================================
1783: | __subsf3
1784: |=============================================================================
1785:
1786: | float __subsf3(float, float);
1787: SYM (__subsf3):
1788: bchg IMM (31),sp@(8) | change sign of second operand
1789: | and fall through
1790: |=============================================================================
1791: | __addsf3
1792: |=============================================================================
1793:
1794: | float __addsf3(float, float);
1795: SYM (__addsf3):
1796: link a6,IMM (0) | everything will be done in registers
1797: moveml d2-d7,sp@- | save all data registers but d0-d1
1798: movel a6@(8),d0 | get first operand
1799: movel a6@(12),d1 | get second operand
1800: movel d0,d6 | get d0's sign bit '
1801: addl d0,d0 | check and clear sign bit of a
1802: beq Laddsf$b | if zero return second operand
1803: movel d1,d7 | save b's sign bit '
1804: addl d1,d1 | get rid of sign bit
1805: beq Laddsf$a | if zero return first operand
1806:
1807: movel d6,a0 | save signs in address registers
1808: movel d7,a1 | so we can use d6 and d7
1809:
1810: | Get the exponents and check for denormalized and/or infinity.
1811:
1812: movel IMM (0x00ffffff),d4 | mask to get fraction
1813: movel IMM (0x01000000),d5 | mask to put hidden bit back
1814:
1815: movel d0,d6 | save a to get exponent
1816: andl d4,d0 | get fraction in d0
1817: notl d4 | make d4 into a mask for the exponent
1818: andl d4,d6 | get exponent in d6
1819: beq Laddsf$a$den | branch if a is denormalized
1820: cmpl d4,d6 | check for INFINITY or NaN
1821: beq Laddsf$nf
1822: swap d6 | put exponent into first word
1823: orl d5,d0 | and put hidden bit back
1824: Laddsf$1:
1825: | Now we have a's exponent in d6 (second byte) and the mantissa in d0. '
1826: movel d1,d7 | get exponent in d7
1827: andl d4,d7 |
1828: beq Laddsf$b$den | branch if b is denormalized
1829: cmpl d4,d7 | check for INFINITY or NaN
1830: beq Laddsf$nf
1831: swap d7 | put exponent into first word
1832: notl d4 | make d4 into a mask for the fraction
1833: andl d4,d1 | get fraction in d1
1834: orl d5,d1 | and put hidden bit back
1835: Laddsf$2:
1836: | Now we have b's exponent in d7 (second byte) and the mantissa in d1. '
1837:
1838: | Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we
1839: | shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
1840: | bit).
1841:
1842: movel d1,d2 | move b to d2, since we want to use
1843: | two registers to do the sum
1844: movel IMM (0),d1 | and clear the new ones
1845: movel d1,d3 |
1846:
1847: | Here we shift the numbers in registers d0 and d1 so the exponents are the
1848: | same, and put the largest exponent in d6. Note that we are using two
1849: | registers for each number (see the discussion by D. Knuth in "Seminumerical
1850: | Algorithms").
1851: cmpw d6,d7 | compare exponents
1852: beq Laddsf$3 | if equal don't shift '
1853: bhi 5f | branch if second exponent largest
1854: 1:
1855: subl d6,d7 | keep the largest exponent
1856: negl d7
1857: lsrw IMM (8),d7 | put difference in lower byte
1858: | if difference is too large we don't shift (actually, we can just exit) '
1859: cmpw IMM (FLT_MANT_DIG+2),d7
1860: bge Laddsf$b$small
1861: cmpw IMM (16),d7 | if difference >= 16 swap
1862: bge 4f
1863: 2:
1864: subw IMM (1),d7
1865: 3: lsrl IMM (1),d2 | shift right second operand
1866: roxrl IMM (1),d3
1867: dbra d7,3b
1868: bra Laddsf$3
1869: 4:
1870: movew d2,d3
1871: swap d3
1872: movew d3,d2
1873: swap d2
1874: subw IMM (16),d7
1875: bne 2b | if still more bits, go back to normal case
1876: bra Laddsf$3
1877: 5:
1878: exg d6,d7 | exchange the exponents
1879: subl d6,d7 | keep the largest exponent
1880: negl d7 |
1881: lsrw IMM (8),d7 | put difference in lower byte
1882: | if difference is too large we don't shift (and exit!) '
1883: cmpw IMM (FLT_MANT_DIG+2),d7
1884: bge Laddsf$a$small
1885: cmpw IMM (16),d7 | if difference >= 16 swap
1886: bge 8f
1887: 6:
1888: subw IMM (1),d7
1889: 7: lsrl IMM (1),d0 | shift right first operand
1890: roxrl IMM (1),d1
1891: dbra d7,7b
1892: bra Laddsf$3
1893: 8:
1894: movew d0,d1
1895: swap d1
1896: movew d1,d0
1897: swap d0
1898: subw IMM (16),d7
1899: bne 6b | if still more bits, go back to normal case
1900: | otherwise we fall through
1901:
1902: | Now we have a in d0-d1, b in d2-d3, and the largest exponent in d6 (the
1903: | signs are stored in a0 and a1).
1904:
1905: Laddsf$3:
1906: | Here we have to decide whether to add or subtract the numbers
1907: exg d6,a0 | get signs back
1908: exg d7,a1 | and save the exponents
1909: eorl d6,d7 | combine sign bits
1910: bmi Lsubsf$0 | if negative a and b have opposite
1911: | sign so we actually subtract the
1912: | numbers
1913:
1914: | Here we have both positive or both negative
1915: exg d6,a0 | now we have the exponent in d6
1916: movel a0,d7 | and sign in d7
1917: andl IMM (0x80000000),d7
1918: | Here we do the addition.
1919: addl d3,d1
1920: addxl d2,d0
1921: | Note: now we have d2, d3, d4 and d5 to play with!
1922:
1923: | Put the exponent, in the first byte, in d2, to use the "standard" rounding
1924: | routines:
1925: movel d6,d2
1926: lsrw IMM (8),d2
1927:
1928: | Before rounding normalize so bit #FLT_MANT_DIG is set (we will consider
1929: | the case of denormalized numbers in the rounding routine itself).
1930: | As in the addition (not in the subtraction!) we could have set
1931: | one more bit we check this:
1932: btst IMM (FLT_MANT_DIG+1),d0
1933: beq 1f
1934: lsrl IMM (1),d0
1935: roxrl IMM (1),d1
1936: addl IMM (1),d2
1937: 1:
1938: lea Laddsf$4,a0 | to return from rounding routine
1939: lea SYM (_fpCCR),a1 | check the rounding mode
1940: movew a1@(6),d6 | rounding mode in d6
1941: beq Lround$to$nearest
1942: cmpw IMM (ROUND_TO_PLUS),d6
1943: bhi Lround$to$minus
1944: blt Lround$to$zero
1945: bra Lround$to$plus
1946: Laddsf$4:
1947: | Put back the exponent, but check for overflow.
1948: cmpw IMM (0xff),d2
1949: bhi 1f
1950: bclr IMM (FLT_MANT_DIG-1),d0
1951: lslw IMM (7),d2
1952: swap d2
1953: orl d2,d0
1954: bra Laddsf$ret
1955: 1:
1956: movew IMM (ADD),d5
1957: bra Lf$overflow
1958:
1959: Lsubsf$0:
1960: | We are here if a > 0 and b < 0 (sign bits cleared).
1961: | Here we do the subtraction.
1962: movel d6,d7 | put sign in d7
1963: andl IMM (0x80000000),d7
1964:
1965: subl d3,d1 | result in d0-d1
1966: subxl d2,d0 |
1967: beq Laddsf$ret | if zero just exit
1968: bpl 1f | if positive skip the following
1969: bchg IMM (31),d7 | change sign bit in d7
1970: negl d1
1971: negxl d0
1972: 1:
1973: exg d2,a0 | now we have the exponent in d2
1974: lsrw IMM (8),d2 | put it in the first byte
1975:
1976: | Now d0-d1 is positive and the sign bit is in d7.
1977:
1978: | Note that we do not have to normalize, since in the subtraction bit
1979: | #FLT_MANT_DIG+1 is never set, and denormalized numbers are handled by
1980: | the rounding routines themselves.
1981: lea Lsubsf$1,a0 | to return from rounding routine
1982: lea SYM (_fpCCR),a1 | check the rounding mode
1983: movew a1@(6),d6 | rounding mode in d6
1984: beq Lround$to$nearest
1985: cmpw IMM (ROUND_TO_PLUS),d6
1986: bhi Lround$to$minus
1987: blt Lround$to$zero
1988: bra Lround$to$plus
1989: Lsubsf$1:
1990: | Put back the exponent (we can't have overflow!). '
1991: bclr IMM (FLT_MANT_DIG-1),d0
1992: lslw IMM (7),d2
1993: swap d2
1994: orl d2,d0
1995: bra Laddsf$ret
1996:
1997: | If one of the numbers was too small (difference of exponents >=
1998: | FLT_MANT_DIG+2) we return the other (and now we don't have to '
1999: | check for finiteness or zero).
2000: Laddsf$a$small:
2001: movel a6@(12),d0
2002: lea SYM (_fpCCR),a0
2003: movew IMM (0),a0@
2004: moveml sp@+,d2-d7 | restore data registers
2005: unlk a6 | and return
2006: rts
2007:
2008: Laddsf$b$small:
2009: movel a6@(8),d0
2010: lea SYM (_fpCCR),a0
2011: movew IMM (0),a0@
2012: moveml sp@+,d2-d7 | restore data registers
2013: unlk a6 | and return
2014: rts
2015:
2016: | If the numbers are denormalized remember to put exponent equal to 1.
2017:
2018: Laddsf$a$den:
2019: movel d5,d6 | d5 contains 0x01000000
2020: swap d6
2021: bra Laddsf$1
2022:
2023: Laddsf$b$den:
2024: movel d5,d7
2025: swap d7
2026: notl d4 | make d4 into a mask for the fraction
2027: | (this was not executed after the jump)
2028: bra Laddsf$2
2029:
2030: | The rest is mainly code for the different results which can be
2031: | returned (checking always for +/-INFINITY and NaN).
2032:
2033: Laddsf$b:
2034: | Return b (if a is zero).
2035: movel a6@(12),d0
2036: bra 1f
2037: Laddsf$a:
2038: | Return a (if b is zero).
2039: movel a6@(8),d0
2040: 1:
2041: movew IMM (ADD),d5
2042: | We have to check for NaN and +/-infty.
2043: movel d0,d7
2044: andl IMM (0x80000000),d7 | put sign in d7
2045: bclr IMM (31),d0 | clear sign
2046: cmpl IMM (INFINITY),d0 | check for infty or NaN
2047: bge 2f
2048: movel d0,d0 | check for zero (we do this because we don't '
2049: bne Laddsf$ret | want to return -0 by mistake
2050: bclr IMM (31),d7 | if zero be sure to clear sign
2051: bra Laddsf$ret | if everything OK just return
2052: 2:
2053: | The value to be returned is either +/-infty or NaN
2054: andl IMM (0x007fffff),d0 | check for NaN
2055: bne Lf$inop | if mantissa not zero is NaN
2056: bra Lf$infty
2057:
2058: Laddsf$ret:
2059: | Normal exit (a and b nonzero, result is not NaN nor +/-infty).
2060: | We have to clear the exception flags (just the exception type).
2061: lea SYM (_fpCCR),a0
2062: movew IMM (0),a0@
2063: orl d7,d0 | put sign bit
2064: moveml sp@+,d2-d7 | restore data registers
2065: unlk a6 | and return
2066: rts
2067:
2068: Laddsf$ret$den:
2069: | Return a denormalized number (for addition we don't signal underflow) '
2070: lsrl IMM (1),d0 | remember to shift right back once
2071: bra Laddsf$ret | and return
2072:
2073: | Note: when adding two floats of the same sign if either one is
2074: | NaN we return NaN without regard to whether the other is finite or
2075: | not. When subtracting them (i.e., when adding two numbers of
2076: | opposite signs) things are more complicated: if both are INFINITY
2077: | we return NaN, if only one is INFINITY and the other is NaN we return
2078: | NaN, but if it is finite we return INFINITY with the corresponding sign.
2079:
2080: Laddsf$nf:
2081: movew IMM (ADD),d5
2082: | This could be faster but it is not worth the effort, since it is not
2083: | executed very often. We sacrifice speed for clarity here.
2084: movel a6@(8),d0 | get the numbers back (remember that we
2085: movel a6@(12),d1 | did some processing already)
2086: movel IMM (INFINITY),d4 | useful constant (INFINITY)
2087: movel d0,d2 | save sign bits
2088: movel d1,d3
2089: bclr IMM (31),d0 | clear sign bits
2090: bclr IMM (31),d1
2091: | We know that one of them is either NaN of +/-INFINITY
2092: | Check for NaN (if either one is NaN return NaN)
2093: cmpl d4,d0 | check first a (d0)
2094: bhi Lf$inop
2095: cmpl d4,d1 | check now b (d1)
2096: bhi Lf$inop
2097: | Now comes the check for +/-INFINITY. We know that both are (maybe not
2098: | finite) numbers, but we have to check if both are infinite whether we
2099: | are adding or subtracting them.
2100: eorl d3,d2 | to check sign bits
2101: bmi 1f
2102: movel d0,d7
2103: andl IMM (0x80000000),d7 | get (common) sign bit
2104: bra Lf$infty
2105: 1:
2106: | We know one (or both) are infinite, so we test for equality between the
2107: | two numbers (if they are equal they have to be infinite both, so we
2108: | return NaN).
2109: cmpl d1,d0 | are both infinite?
2110: beq Lf$inop | if so return NaN
2111:
2112: movel d0,d7
2113: andl IMM (0x80000000),d7 | get a's sign bit '
2114: cmpl d4,d0 | test now for infinity
2115: beq Lf$infty | if a is INFINITY return with this sign
2116: bchg IMM (31),d7 | else we know b is INFINITY and has
2117: bra Lf$infty | the opposite sign
2118:
2119: |=============================================================================
2120: | __mulsf3
2121: |=============================================================================
2122:
2123: | float __mulsf3(float, float);
2124: SYM (__mulsf3):
2125: link a6,IMM (0)
2126: moveml d2-d7,sp@-
2127: movel a6@(8),d0 | get a into d0
2128: movel a6@(12),d1 | and b into d1
2129: movel d0,d7 | d7 will hold the sign of the product
2130: eorl d1,d7 |
2131: andl IMM (0x80000000),d7
2132: movel IMM (INFINITY),d6 | useful constant (+INFINITY)
2133: movel d6,d5 | another (mask for fraction)
2134: notl d5 |
2135: movel IMM (0x00800000),d4 | this is to put hidden bit back
2136: bclr IMM (31),d0 | get rid of a's sign bit '
2137: movel d0,d2 |
2138: beq Lmulsf$a$0 | branch if a is zero
2139: bclr IMM (31),d1 | get rid of b's sign bit '
2140: movel d1,d3 |
2141: beq Lmulsf$b$0 | branch if b is zero
2142: cmpl d6,d0 | is a big?
2143: bhi Lmulsf$inop | if a is NaN return NaN
2144: beq Lmulsf$inf | if a is INFINITY we have to check b
2145: cmpl d6,d1 | now compare b with INFINITY
2146: bhi Lmulsf$inop | is b NaN?
2147: beq Lmulsf$overflow | is b INFINITY?
2148: | Here we have both numbers finite and nonzero (and with no sign bit).
2149: | Now we get the exponents into d2 and d3.
2150: andl d6,d2 | and isolate exponent in d2
2151: beq Lmulsf$a$den | if exponent is zero we have a denormalized
2152: andl d5,d0 | and isolate fraction
2153: orl d4,d0 | and put hidden bit back
2154: swap d2 | I like exponents in the first byte
2155: lsrw IMM (7),d2 |
2156: Lmulsf$1: | number
2157: andl d6,d3 |
2158: beq Lmulsf$b$den |
2159: andl d5,d1 |
2160: orl d4,d1 |
2161: swap d3 |
2162: lsrw IMM (7),d3 |
2163: Lmulsf$2: |
2164: addw d3,d2 | add exponents
2165: subw IMM (F_BIAS+1),d2 | and subtract bias (plus one)
2166:
2167: | We are now ready to do the multiplication. The situation is as follows:
2168: | both a and b have bit FLT_MANT_DIG-1 set (even if they were
2169: | denormalized to start with!), which means that in the product
2170: | bit 2*(FLT_MANT_DIG-1) (that is, bit 2*FLT_MANT_DIG-2-32 of the
2171: | high long) is set.
2172:
2173: | To do the multiplication let us move the number a little bit around ...
2174: movel d1,d6 | second operand in d6
2175: movel d0,d5 | first operand in d4-d5
2176: movel IMM (0),d4
2177: movel d4,d1 | the sums will go in d0-d1
2178: movel d4,d0
2179:
2180: | now bit FLT_MANT_DIG-1 becomes bit 31:
2181: lsll IMM (31-FLT_MANT_DIG+1),d6
2182:
2183: | Start the loop (we loop #FLT_MANT_DIG times):
2184: movew IMM (FLT_MANT_DIG-1),d3
2185: 1: addl d1,d1 | shift sum
2186: addxl d0,d0
2187: lsll IMM (1),d6 | get bit bn
2188: bcc 2f | if not set skip sum
2189: addl d5,d1 | add a
2190: addxl d4,d0
2191: 2: dbf d3,1b | loop back
2192:
2193: | Now we have the product in d0-d1, with bit (FLT_MANT_DIG - 1) + FLT_MANT_DIG
2194: | (mod 32) of d0 set. The first thing to do now is to normalize it so bit
2195: | FLT_MANT_DIG is set (to do the rounding).
2196: rorl IMM (6),d1
2197: swap d1
2198: movew d1,d3
2199: andw IMM (0x03ff),d3
2200: andw IMM (0xfd00),d1
2201: lsll IMM (8),d0
2202: addl d0,d0
2203: addl d0,d0
2204: orw d3,d0
2205:
2206: movew IMM (MULTIPLY),d5
2207:
2208: btst IMM (FLT_MANT_DIG+1),d0
2209: beq Lround$exit
2210: lsrl IMM (1),d0
2211: roxrl IMM (1),d1
2212: addw IMM (1),d2
2213: bra Lround$exit
2214:
2215: Lmulsf$inop:
2216: movew IMM (MULTIPLY),d5
2217: bra Lf$inop
2218:
2219: Lmulsf$overflow:
2220: movew IMM (MULTIPLY),d5
2221: bra Lf$overflow
2222:
2223: Lmulsf$inf:
2224: movew IMM (MULTIPLY),d5
2225: | If either is NaN return NaN; else both are (maybe infinite) numbers, so
2226: | return INFINITY with the correct sign (which is in d7).
2227: cmpl d6,d1 | is b NaN?
2228: bhi Lf$inop | if so return NaN
2229: bra Lf$overflow | else return +/-INFINITY
2230:
2231: | If either number is zero return zero, unless the other is +/-INFINITY,
2232: | or NaN, in which case we return NaN.
2233: Lmulsf$b$0:
2234: | Here d1 (==b) is zero.
2235: movel d1,d0 | put b into d0 (just a zero)
2236: movel a6@(8),d1 | get a again to check for non-finiteness
2237: bra 1f
2238: Lmulsf$a$0:
2239: movel a6@(12),d1 | get b again to check for non-finiteness
2240: 1: bclr IMM (31),d1 | clear sign bit
2241: cmpl IMM (INFINITY),d1 | and check for a large exponent
2242: bge Lf$inop | if b is +/-INFINITY or NaN return NaN
2243: lea SYM (_fpCCR),a0 | else return zero
2244: movew IMM (0),a0@ |
2245: moveml sp@+,d2-d7 |
2246: unlk a6 |
2247: rts |
2248:
2249: | If a number is denormalized we put an exponent of 1 but do not put the
2250: | hidden bit back into the fraction; instead we shift left until bit 23
2251: | (the hidden bit) is set, adjusting the exponent accordingly. We do this
2252: | to ensure that the product of the fractions is close to 1.
2253: Lmulsf$a$den:
2254: movel IMM (1),d2
2255: andl d5,d0
2256: 1: addl d0,d0 | shift a left (until bit 23 is set)
2257: subw IMM (1),d2 | and adjust exponent
2258: btst IMM (FLT_MANT_DIG-1),d0
2259: bne Lmulsf$1 |
2260: bra 1b | else loop back
2261:
2262: Lmulsf$b$den:
2263: movel IMM (1),d3
2264: andl d5,d1
2265: 1: addl d1,d1 | shift b left until bit 23 is set
2266: subw IMM (1),d3 | and adjust exponent
2267: btst IMM (FLT_MANT_DIG-1),d1
2268: bne Lmulsf$2 |
2269: bra 1b | else loop back
2270:
2271: |=============================================================================
2272: | __divsf3
2273: |=============================================================================
2274:
2275: | float __divsf3(float, float);
2276: SYM (__divsf3):
2277: link a6,IMM (0)
2278: moveml d2-d7,sp@-
2279: movel a6@(8),d0 | get a into d0
2280: movel a6@(12),d1 | and b into d1
2281: movel d0,d7 | d7 will hold the sign of the result
2282: eorl d1,d7 |
2283: andl IMM (0x80000000),d7 |
2284: movel IMM (INFINITY),d6 | useful constant (+INFINITY)
2285: movel d6,d5 | another (mask for fraction)
2286: notl d5 |
2287: movel IMM (0x00800000),d4 | this is to put hidden bit back
2288: bclr IMM (31),d0 | get rid of a's sign bit '
2289: movel d0,d2 |
2290: beq Ldivsf$a$0 | branch if a is zero
2291: bclr IMM (31),d1 | get rid of b's sign bit '
2292: movel d1,d3 |
2293: beq Ldivsf$b$0 | branch if b is zero
2294: cmpl d6,d0 | is a big?
2295: bhi Ldivsf$inop | if a is NaN return NaN
2296: beq Ldivsf$inf | if a is INFINITY we have to check b
2297: cmpl d6,d1 | now compare b with INFINITY
2298: bhi Ldivsf$inop | if b is NaN return NaN
2299: beq Ldivsf$underflow
2300: | Here we have both numbers finite and nonzero (and with no sign bit).
2301: | Now we get the exponents into d2 and d3 and normalize the numbers to
2302: | ensure that the ratio of the fractions is close to 1. We do this by
2303: | making sure that bit #FLT_MANT_DIG-1 (hidden bit) is set.
2304: andl d6,d2 | and isolate exponent in d2
2305: beq Ldivsf$a$den | if exponent is zero we have a denormalized
2306: andl d5,d0 | and isolate fraction
2307: orl d4,d0 | and put hidden bit back
2308: swap d2 | I like exponents in the first byte
2309: lsrw IMM (7),d2 |
2310: Ldivsf$1: |
2311: andl d6,d3 |
2312: beq Ldivsf$b$den |
2313: andl d5,d1 |
2314: orl d4,d1 |
2315: swap d3 |
2316: lsrw IMM (7),d3 |
2317: Ldivsf$2: |
2318: subw d3,d2 | subtract exponents
2319: addw IMM (F_BIAS),d2 | and add bias
2320:
2321: | We are now ready to do the division. We have prepared things in such a way
2322: | that the ratio of the fractions will be less than 2 but greater than 1/2.
2323: | At this point the registers in use are:
2324: | d0 holds a (first operand, bit FLT_MANT_DIG=0, bit FLT_MANT_DIG-1=1)
2325: | d1 holds b (second operand, bit FLT_MANT_DIG=1)
2326: | d2 holds the difference of the exponents, corrected by the bias
2327: | d7 holds the sign of the ratio
2328: | d4, d5, d6 hold some constants
2329: movel d7,a0 | d6-d7 will hold the ratio of the fractions
2330: movel IMM (0),d6 |
2331: movel d6,d7
2332:
2333: movew IMM (FLT_MANT_DIG+1),d3
2334: 1: cmpl d0,d1 | is a < b?
2335: bhi 2f |
2336: bset d3,d6 | set a bit in d6
2337: subl d1,d0 | if a >= b a <-- a-b
2338: beq 3f | if a is zero, exit
2339: 2: addl d0,d0 | multiply a by 2
2340: dbra d3,1b
2341:
2342: | Now we keep going to set the sticky bit ...
2343: movew IMM (FLT_MANT_DIG),d3
2344: 1: cmpl d0,d1
2345: ble 2f
2346: addl d0,d0
2347: dbra d3,1b
2348: movel IMM (0),d1
2349: bra 3f
2350: 2: movel IMM (0),d1
2351: subw IMM (FLT_MANT_DIG),d3
2352: addw IMM (31),d3
2353: bset d3,d1
2354: 3:
2355: movel d6,d0 | put the ratio in d0-d1
2356: movel a0,d7 | get sign back
2357:
2358: | Because of the normalization we did before we are guaranteed that
2359: | d0 is smaller than 2^26 but larger than 2^24. Thus bit 26 is not set,
2360: | bit 25 could be set, and if it is not set then bit 24 is necessarily set.
2361: btst IMM (FLT_MANT_DIG+1),d0
2362: beq 1f | if it is not set, then bit 24 is set
2363: lsrl IMM (1),d0 |
2364: addw IMM (1),d2 |
2365: 1:
2366: | Now round, check for over- and underflow, and exit.
2367: movew IMM (DIVIDE),d5
2368: bra Lround$exit
2369:
2370: Ldivsf$inop:
2371: movew IMM (DIVIDE),d5
2372: bra Lf$inop
2373:
2374: Ldivsf$overflow:
2375: movew IMM (DIVIDE),d5
2376: bra Lf$overflow
2377:
2378: Ldivsf$underflow:
2379: movew IMM (DIVIDE),d5
2380: bra Lf$underflow
2381:
2382: Ldivsf$a$0:
2383: movew IMM (DIVIDE),d5
2384: | If a is zero check to see whether b is zero also. In that case return
2385: | NaN; then check if b is NaN, and return NaN also in that case. Else
2386: | return zero.
2387: andl IMM (0x7fffffff),d1 | clear sign bit and test b
2388: beq Lf$inop | if b is also zero return NaN
2389: cmpl IMM (INFINITY),d1 | check for NaN
2390: bhi Lf$inop |
2391: movel IMM (0),d0 | else return zero
2392: lea SYM (_fpCCR),a0 |
2393: movew IMM (0),a0@ |
2394: moveml sp@+,d2-d7 |
2395: unlk a6 |
2396: rts |
2397:
2398: Ldivsf$b$0:
2399: movew IMM (DIVIDE),d5
2400: | If we got here a is not zero. Check if a is NaN; in that case return NaN,
2401: | else return +/-INFINITY. Remember that a is in d0 with the sign bit
2402: | cleared already.
2403: cmpl IMM (INFINITY),d0 | compare d0 with INFINITY
2404: bhi Lf$inop | if larger it is NaN
2405: bra Lf$div$0 | else signal DIVIDE_BY_ZERO
2406:
2407: Ldivsf$inf:
2408: movew IMM (DIVIDE),d5
2409: | If a is INFINITY we have to check b
2410: cmpl IMM (INFINITY),d1 | compare b with INFINITY
2411: bge Lf$inop | if b is NaN or INFINITY return NaN
2412: bra Lf$overflow | else return overflow
2413:
2414: | If a number is denormalized we put an exponent of 1 but do not put the
2415: | bit back into the fraction.
2416: Ldivsf$a$den:
2417: movel IMM (1),d2
2418: andl d5,d0
2419: 1: addl d0,d0 | shift a left until bit FLT_MANT_DIG-1 is set
2420: subw IMM (1),d2 | and adjust exponent
2421: btst IMM (FLT_MANT_DIG-1),d0
2422: bne Ldivsf$1
2423: bra 1b
2424:
2425: Ldivsf$b$den:
2426: movel IMM (1),d3
2427: andl d5,d1
2428: 1: addl d1,d1 | shift b left until bit FLT_MANT_DIG is set
2429: subw IMM (1),d3 | and adjust exponent
2430: btst IMM (FLT_MANT_DIG-1),d1
2431: bne Ldivsf$2
2432: bra 1b
2433:
2434: Lround$exit:
2435: | This is a common exit point for __mulsf3 and __divsf3.
2436:
2437: | First check for underlow in the exponent:
2438: cmpw IMM (-FLT_MANT_DIG-1),d2
2439: blt Lf$underflow
2440: | It could happen that the exponent is less than 1, in which case the
2441: | number is denormalized. In this case we shift right and adjust the
2442: | exponent until it becomes 1 or the fraction is zero (in the latter case
2443: | we signal underflow and return zero).
2444: movel IMM (0),d6 | d6 is used temporarily
2445: cmpw IMM (1),d2 | if the exponent is less than 1 we
2446: bge 2f | have to shift right (denormalize)
2447: 1: addw IMM (1),d2 | adjust the exponent
2448: lsrl IMM (1),d0 | shift right once
2449: roxrl IMM (1),d1 |
2450: roxrl IMM (1),d6 | d6 collect bits we would lose otherwise
2451: cmpw IMM (1),d2 | is the exponent 1 already?
2452: beq 2f | if not loop back
2453: bra 1b |
2454: bra Lf$underflow | safety check, shouldn't execute '
2455: 2: orl d6,d1 | this is a trick so we don't lose '
2456: | the extra bits which were flushed right
2457: | Now call the rounding routine (which takes care of denormalized numbers):
2458: lea Lround$0,a0 | to return from rounding routine
2459: lea SYM (_fpCCR),a1 | check the rounding mode
2460: movew a1@(6),d6 | rounding mode in d6
2461: beq Lround$to$nearest
2462: cmpw IMM (ROUND_TO_PLUS),d6
2463: bhi Lround$to$minus
2464: blt Lround$to$zero
2465: bra Lround$to$plus
2466: Lround$0:
2467: | Here we have a correctly rounded result (either normalized or denormalized).
2468:
2469: | Here we should have either a normalized number or a denormalized one, and
2470: | the exponent is necessarily larger or equal to 1 (so we don't have to '
2471: | check again for underflow!). We have to check for overflow or for a
2472: | denormalized number (which also signals underflow).
2473: | Check for overflow (i.e., exponent >= 255).
2474: cmpw IMM (0x00ff),d2
2475: bge Lf$overflow
2476: | Now check for a denormalized number (exponent==0).
2477: movew d2,d2
2478: beq Lf$den
2479: 1:
2480: | Put back the exponents and sign and return.
2481: lslw IMM (7),d2 | exponent back to fourth byte
2482: bclr IMM (FLT_MANT_DIG-1),d0
2483: swap d0 | and put back exponent
2484: orw d2,d0 |
2485: swap d0 |
2486: orl d7,d0 | and sign also
2487:
2488: lea SYM (_fpCCR),a0
2489: movew IMM (0),a0@
2490: moveml sp@+,d2-d7
2491: unlk a6
2492: rts
2493:
2494: |=============================================================================
2495: | __negsf2
2496: |=============================================================================
2497:
2498: | This is trivial and could be shorter if we didn't bother checking for NaN '
2499: | and +/-INFINITY.
2500:
2501: | float __negsf2(float);
2502: SYM (__negsf2):
2503: link a6,IMM (0)
2504: moveml d2-d7,sp@-
2505: movew IMM (NEGATE),d5
2506: movel a6@(8),d0 | get number to negate in d0
2507: bchg IMM (31),d0 | negate
2508: movel d0,d1 | make a positive copy
2509: bclr IMM (31),d1 |
2510: tstl d1 | check for zero
2511: beq 2f | if zero (either sign) return +zero
2512: cmpl IMM (INFINITY),d1 | compare to +INFINITY
2513: blt 1f |
2514: bhi Lf$inop | if larger (fraction not zero) is NaN
2515: movel d0,d7 | else get sign and return INFINITY
2516: andl IMM (0x80000000),d7
2517: bra Lf$infty
2518: 1: lea SYM (_fpCCR),a0
2519: movew IMM (0),a0@
2520: moveml sp@+,d2-d7
2521: unlk a6
2522: rts
2523: 2: bclr IMM (31),d0
2524: bra 1b
2525:
2526: |=============================================================================
2527: | __cmpsf2
2528: |=============================================================================
2529:
2530: GREATER = 1
2531: LESS = -1
2532: EQUAL = 0
2533:
2534: | int __cmpsf2(float, float);
2535: SYM (__cmpsf2):
2536: link a6,IMM (0)
2537: moveml d2-d7,sp@- | save registers
2538: movew IMM (COMPARE),d5
2539: movel a6@(8),d0 | get first operand
2540: movel a6@(12),d1 | get second operand
2541: | Check if either is NaN, and in that case return garbage and signal
2542: | INVALID_OPERATION. Check also if either is zero, and clear the signs
2543: | if necessary.
2544: movel d0,d6
2545: andl IMM (0x7fffffff),d0
2546: beq Lcmpsf$a$0
2547: cmpl IMM (0x7f800000),d0
2548: bhi Lf$inop
2549: Lcmpsf$1:
2550: movel d1,d7
2551: andl IMM (0x7fffffff),d1
2552: beq Lcmpsf$b$0
2553: cmpl IMM (0x7f800000),d1
2554: bhi Lf$inop
2555: Lcmpsf$2:
2556: | Check the signs
2557: eorl d6,d7
2558: bpl 1f
2559: | If the signs are not equal check if a >= 0
2560: tstl d6
2561: bpl Lcmpsf$a$gt$b | if (a >= 0 && b < 0) => a > b
2562: bmi Lcmpsf$b$gt$a | if (a < 0 && b >= 0) => a < b
2563: 1:
2564: | If the signs are equal check for < 0
2565: tstl d6
2566: bpl 1f
2567: | If both are negative exchange them
2568: exg d0,d1
2569: 1:
2570: | Now that they are positive we just compare them as longs (does this also
2571: | work for denormalized numbers?).
2572: cmpl d0,d1
2573: bhi Lcmpsf$b$gt$a | |b| > |a|
2574: bne Lcmpsf$a$gt$b | |b| < |a|
2575: | If we got here a == b.
2576: movel IMM (EQUAL),d0
2577: moveml sp@+,d2-d7 | put back the registers
2578: unlk a6
2579: rts
2580: Lcmpsf$a$gt$b:
2581: movel IMM (GREATER),d0
2582: moveml sp@+,d2-d7 | put back the registers
2583: unlk a6
2584: rts
2585: Lcmpsf$b$gt$a:
2586: movel IMM (LESS),d0
2587: moveml sp@+,d2-d7 | put back the registers
2588: unlk a6
2589: rts
2590:
2591: Lcmpsf$a$0:
2592: bclr IMM (31),d6
2593: bra Lcmpsf$1
2594: Lcmpsf$b$0:
2595: bclr IMM (31),d7
2596: bra Lcmpsf$2
2597:
2598: |=============================================================================
2599: | rounding routines
2600: |=============================================================================
2601:
2602: | The rounding routines expect the number to be normalized in registers
2603: | d0-d1, with the exponent in register d2. They assume that the
2604: | exponent is larger or equal to 1. They return a properly normalized number
2605: | if possible, and a denormalized number otherwise. The exponent is returned
2606: | in d2.
2607:
2608: Lround$to$nearest:
2609: | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
2610: | Here we assume that the exponent is not too small (this should be checked
2611: | before entering the rounding routine), but the number could be denormalized.
2612:
2613: | Check for denormalized numbers:
2614: 1: btst IMM (FLT_MANT_DIG),d0
2615: bne 2f | if set the number is normalized
2616: | Normalize shifting left until bit #FLT_MANT_DIG is set or the exponent
2617: | is one (remember that a denormalized number corresponds to an
2618: | exponent of -F_BIAS+1).
2619: cmpw IMM (1),d2 | remember that the exponent is at least one
2620: beq 2f | an exponent of one means denormalized
2621: addl d1,d1 | else shift and adjust the exponent
2622: addxl d0,d0 |
2623: dbra d2,1b |
2624: 2:
2625: | Now round: we do it as follows: after the shifting we can write the
2626: | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
2627: | If delta < 1, do nothing. If delta > 1, add 1 to f.
2628: | If delta == 1, we make sure the rounded number will be even (odd?)
2629: | (after shifting).
2630: btst IMM (0),d0 | is delta < 1?
2631: beq 2f | if so, do not do anything
2632: tstl d1 | is delta == 1?
2633: bne 1f | if so round to even
2634: movel d0,d1 |
2635: andl IMM (2),d1 | bit 1 is the last significant bit
2636: addl d1,d0 |
2637: bra 2f |
2638: 1: movel IMM (1),d1 | else add 1
2639: addl d1,d0 |
2640: | Shift right once (because we used bit #FLT_MANT_DIG!).
2641: 2: lsrl IMM (1),d0
2642: | Now check again bit #FLT_MANT_DIG (rounding could have produced a
2643: | 'fraction overflow' ...).
2644: btst IMM (FLT_MANT_DIG),d0
2645: beq 1f
2646: lsrl IMM (1),d0
2647: addw IMM (1),d2
2648: 1:
2649: | If bit #FLT_MANT_DIG-1 is clear we have a denormalized number, so we
2650: | have to put the exponent to zero and return a denormalized number.
2651: btst IMM (FLT_MANT_DIG-1),d0
2652: beq 1f
2653: jmp a0@
2654: 1: movel IMM (0),d2
2655: jmp a0@
2656:
2657: Lround$to$zero:
2658: Lround$to$plus:
2659: Lround$to$minus:
2660: jmp a0@
2661: #endif /* L_float */
2662:
2663: | gcc expects the routines __eqdf2, __nedf2, __gtdf2, __gedf2,
2664: | __ledf2, __ltdf2 to all return the same value as a direct call to
2665: | __cmpdf2 would. In this implementation, each of these routines
2666: | simply calls __cmpdf2. It would be more efficient to give the
2667: | __cmpdf2 routine several names, but separating them out will make it
2668: | easier to write efficient versions of these routines someday.
2669:
2670: #ifdef L_eqdf2
2671: LL0:
2672: .text
2673: .proc
2674: |#PROC# 04
2675: LF18 = 4
2676: LS18 = 128
2677: LFF18 = 0
2678: LSS18 = 0
2679: LV18 = 0
2680: .text
2681: .globl SYM (__eqdf2)
2682: SYM (__eqdf2):
2683: |#PROLOGUE# 0
2684: link a6,IMM (0)
2685: |#PROLOGUE# 1
2686: movl a6@(20),sp@-
2687: movl a6@(16),sp@-
2688: movl a6@(12),sp@-
2689: movl a6@(8),sp@-
2690: jbsr SYM (__cmpdf2)
2691: |#PROLOGUE# 2
2692: unlk a6
2693: |#PROLOGUE# 3
2694: rts
2695: #endif /* L_eqdf2 */
2696:
2697: #ifdef L_nedf2
2698: LL0:
2699: .text
2700: .proc
2701: |#PROC# 04
2702: LF18 = 8
2703: LS18 = 132
2704: LFF18 = 0
2705: LSS18 = 0
2706: LV18 = 0
2707: .text
2708: .globl SYM (__nedf2)
2709: SYM (__nedf2):
2710: |#PROLOGUE# 0
2711: link a6,IMM (0)
2712: |#PROLOGUE# 1
2713: movl a6@(20),sp@-
2714: movl a6@(16),sp@-
2715: movl a6@(12),sp@-
2716: movl a6@(8),sp@-
2717: jbsr SYM (__cmpdf2)
2718: |#PROLOGUE# 2
2719: unlk a6
2720: |#PROLOGUE# 3
2721: rts
2722: #endif /* L_nedf2 */
2723:
2724: #ifdef L_gtdf2
2725: .text
2726: .proc
2727: |#PROC# 04
2728: LF18 = 8
2729: LS18 = 132
2730: LFF18 = 0
2731: LSS18 = 0
2732: LV18 = 0
2733: .text
2734: .globl SYM (__gtdf2)
2735: SYM (__gtdf2):
2736: |#PROLOGUE# 0
2737: link a6,IMM (0)
2738: |#PROLOGUE# 1
2739: movl a6@(20),sp@-
2740: movl a6@(16),sp@-
2741: movl a6@(12),sp@-
2742: movl a6@(8),sp@-
2743: jbsr SYM (__cmpdf2)
2744: |#PROLOGUE# 2
2745: unlk a6
2746: |#PROLOGUE# 3
2747: rts
2748: #endif /* L_gtdf2 */
2749:
2750: #ifdef L_gedf2
2751: LL0:
2752: .text
2753: .proc
2754: |#PROC# 04
2755: LF18 = 8
2756: LS18 = 132
2757: LFF18 = 0
2758: LSS18 = 0
2759: LV18 = 0
2760: .text
2761: .globl SYM (__gedf2)
2762: SYM (__gedf2):
2763: |#PROLOGUE# 0
2764: link a6,IMM (0)
2765: |#PROLOGUE# 1
2766: movl a6@(20),sp@-
2767: movl a6@(16),sp@-
2768: movl a6@(12),sp@-
2769: movl a6@(8),sp@-
2770: jbsr SYM (__cmpdf2)
2771: |#PROLOGUE# 2
2772: unlk a6
2773: |#PROLOGUE# 3
2774: rts
2775: #endif /* L_gedf2 */
2776:
2777: #ifdef L_ltdf2
2778: LL0:
2779: .text
2780: .proc
2781: |#PROC# 04
2782: LF18 = 8
2783: LS18 = 132
2784: LFF18 = 0
2785: LSS18 = 0
2786: LV18 = 0
2787: .text
2788: .globl SYM (__ltdf2)
2789: SYM (__ltdf2):
2790: |#PROLOGUE# 0
2791: link a6,IMM (0)
2792: |#PROLOGUE# 1
2793: movl a6@(20),sp@-
2794: movl a6@(16),sp@-
2795: movl a6@(12),sp@-
2796: movl a6@(8),sp@-
2797: jbsr SYM (__cmpdf2)
2798: |#PROLOGUE# 2
2799: unlk a6
2800: |#PROLOGUE# 3
2801: rts
2802: #endif /* L_ltdf2 */
2803:
2804: #ifdef L_ledf2
2805: .text
2806: .proc
2807: |#PROC# 04
2808: LF18 = 8
2809: LS18 = 132
2810: LFF18 = 0
2811: LSS18 = 0
2812: LV18 = 0
2813: .text
2814: .globl SYM (__ledf2)
2815: SYM (__ledf2):
2816: |#PROLOGUE# 0
2817: link a6,IMM (0)
2818: |#PROLOGUE# 1
2819: movl a6@(20),sp@-
2820: movl a6@(16),sp@-
2821: movl a6@(12),sp@-
2822: movl a6@(8),sp@-
2823: jbsr SYM (__cmpdf2)
2824: |#PROLOGUE# 2
2825: unlk a6
2826: |#PROLOGUE# 3
2827: rts
2828: #endif /* L_ledf2 */
2829:
2830: | The comments above about __eqdf2, et. al., also apply to __eqsf2,
2831: | et. al., except that the latter call __cmpsf2 rather than __cmpdf2.
2832:
2833: #ifdef L_eqsf2
2834: .text
2835: .proc
2836: |#PROC# 04
2837: LF18 = 4
2838: LS18 = 128
2839: LFF18 = 0
2840: LSS18 = 0
2841: LV18 = 0
2842: .text
2843: .globl SYM (__eqsf2)
2844: SYM (__eqsf2):
2845: |#PROLOGUE# 0
2846: link a6,IMM (0)
2847: |#PROLOGUE# 1
2848: movl a6@(12),sp@-
2849: movl a6@(8),sp@-
2850: jbsr SYM (__cmpsf2)
2851: |#PROLOGUE# 2
2852: unlk a6
2853: |#PROLOGUE# 3
2854: rts
2855: #endif /* L_eqsf2 */
2856:
2857: #ifdef L_nesf2
2858: .text
2859: .proc
2860: |#PROC# 04
2861: LF18 = 8
2862: LS18 = 132
2863: LFF18 = 0
2864: LSS18 = 0
2865: LV18 = 0
2866: .text
2867: .globl SYM (__nesf2)
2868: SYM (__nesf2):
2869: |#PROLOGUE# 0
2870: link a6,IMM (0)
2871: |#PROLOGUE# 1
2872: movl a6@(12),sp@-
2873: movl a6@(8),sp@-
2874: jbsr SYM (__cmpsf2)
2875: |#PROLOGUE# 2
2876: unlk a6
2877: |#PROLOGUE# 3
2878: rts
2879: #endif /* L_nesf2 */
2880:
2881: #ifdef L_gtsf2
2882: .text
2883: .proc
2884: |#PROC# 04
2885: LF18 = 8
2886: LS18 = 132
2887: LFF18 = 0
2888: LSS18 = 0
2889: LV18 = 0
2890: .text
2891: .globl SYM (__gtsf2)
2892: SYM (__gtsf2):
2893: |#PROLOGUE# 0
2894: link a6,IMM (0)
2895: |#PROLOGUE# 1
2896: movl a6@(12),sp@-
2897: movl a6@(8),sp@-
2898: jbsr SYM (__cmpsf2)
2899: |#PROLOGUE# 2
2900: unlk a6
2901: |#PROLOGUE# 3
2902: rts
2903: #endif /* L_gtsf2 */
2904:
2905: #ifdef L_gesf2
2906: .text
2907: .proc
2908: |#PROC# 04
2909: LF18 = 8
2910: LS18 = 132
2911: LFF18 = 0
2912: LSS18 = 0
2913: LV18 = 0
2914: .text
2915: .globl SYM (__gesf2)
2916: SYM (__gesf2):
2917: |#PROLOGUE# 0
2918: link a6,IMM (0)
2919: |#PROLOGUE# 1
2920: movl a6@(12),sp@-
2921: movl a6@(8),sp@-
2922: jbsr SYM (__cmpsf2)
2923: |#PROLOGUE# 2
2924: unlk a6
2925: |#PROLOGUE# 3
2926: rts
2927: #endif /* L_gesf2 */
2928:
2929: #ifdef L_ltsf2
2930: .text
2931: .proc
2932: |#PROC# 04
2933: LF18 = 8
2934: LS18 = 132
2935: LFF18 = 0
2936: LSS18 = 0
2937: LV18 = 0
2938: .text
2939: .globl SYM (__ltsf2)
2940: SYM (__ltsf2):
2941: |#PROLOGUE# 0
2942: link a6,IMM (0)
2943: |#PROLOGUE# 1
2944: movl a6@(12),sp@-
2945: movl a6@(8),sp@-
2946: jbsr SYM (__cmpsf2)
2947: |#PROLOGUE# 2
2948: unlk a6
2949: |#PROLOGUE# 3
2950: rts
2951: #endif /* L_ltsf2 */
2952:
2953: #ifdef L_lesf2
2954: .text
2955: .proc
2956: |#PROC# 04
2957: LF18 = 8
2958: LS18 = 132
2959: LFF18 = 0
2960: LSS18 = 0
2961: LV18 = 0
2962: .text
2963: .globl SYM (__lesf2)
2964: SYM (__lesf2):
2965: |#PROLOGUE# 0
2966: link a6,IMM (0)
2967: |#PROLOGUE# 1
2968: movl a6@(12),sp@-
2969: movl a6@(8),sp@-
2970: jbsr SYM (__cmpsf2)
2971: |#PROLOGUE# 2
2972: unlk a6
2973: |#PROLOGUE# 3
2974: rts
2975: #endif /* L_lesf2 */
2976:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.