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