This is the mail archive of the
gcc-bugs@gcc.gnu.org
mailing list for the GCC project.
gcc 3.4 uclinux Floating Point Patch
- From: Felix Daners <f dot daners at swissworld dot com>
- To: gcc-bugs at gcc dot gnu dot org
- Cc: bernie at develer dot com
- Date: Tue, 05 Dec 2006 23:29:20 +0100
- Subject: gcc 3.4 uclinux Floating Point Patch
Hello
Recently I have tested floating point behaviour of gcc 3.4 uclinux
target built by the script you contribute on
http://www.uclinux.org/pub/uClinux/uclinux-elf-tools/gcc-3/
The result was not very satisfactory. Results are ok but wrong
in overflow / underflow conditions. Long doubles are supported but
in fact degenerate to doubles in all the elementary math functions.
The patch does the following:
- define long doubles as doubles, ie. the size of a long double
is 64 resulting in math precision to meet the number of bits used
to represent it
- removes the bogus floating point routines in fpgnulib.c
- removes most of the asm functions in lb1sf68.asm
- uses the gcc fp-bit.c softfloat library
- removes a prototype problem in strsignal.c
For testing softfloat precision and IEE conformance I have used
http://www.netlib.org/paranoia/paranoia.c
Feel free to add this patch to your repository and if you like
to the gcc community.
Felix
--
+-------------------------------------
| Felix Daners Engineering
| F.Daners
| Platz 3
| CH-8200 Schaffhausen
| Tel: +41 52 624 92 32
| Fax: +41 52 624 92 31
| mailto:f.daners@swissworld.com
| http://f.daners.swissworld.com
+-------------------------------------
diff -urNa gcc-3.4.0/gcc/config/m68k/lb1sf68.asm gcc-3.4.0.mod/gcc/config/m68k/lb1sf68.asm
--- gcc-3.4.0/gcc/config/m68k/lb1sf68.asm 2003-10-11 16:22:54.000000000 +0200
+++ gcc-3.4.0.mod/gcc/config/m68k/lb1sf68.asm 2005-12-16 23:33:35.000000000 +0100
@@ -157,213 +157,6 @@
#endif /* __PIC__ */
-#ifdef L_floatex
-
-| This is an attempt at a decent floating point (single, double and
-| extended double) code for the GNU C compiler. It should be easy to
-| adapt to other compilers (but beware of the local labels!).
-
-| Starting date: 21 October, 1990
-
-| It is convenient to introduce the notation (s,e,f) for a floating point
-| number, where s=sign, e=exponent, f=fraction. We will call a floating
-| point number fpn to abbreviate, independently of the precision.
-| Let MAX_EXP be in each case the maximum exponent (255 for floats, 1023
-| for doubles and 16383 for long doubles). We then have the following
-| different cases:
-| 1. Normalized fpns have 0 < e < MAX_EXP. They correspond to
-| (-1)^s x 1.f x 2^(e-bias-1).
-| 2. Denormalized fpns have e=0. They correspond to numbers of the form
-| (-1)^s x 0.f x 2^(-bias).
-| 3. +/-INFINITY have e=MAX_EXP, f=0.
-| 4. Quiet NaN (Not a Number) have all bits set.
-| 5. Signaling NaN (Not a Number) have s=0, e=MAX_EXP, f=1.
-
-|=============================================================================
-| exceptions
-|=============================================================================
-
-| This is the floating point condition code register (_fpCCR):
-|
-| struct {
-| short _exception_bits;
-| short _trap_enable_bits;
-| short _sticky_bits;
-| short _rounding_mode;
-| short _format;
-| short _last_operation;
-| union {
-| float sf;
-| double df;
-| } _operand1;
-| union {
-| float sf;
-| double df;
-| } _operand2;
-| } _fpCCR;
-
- .data
- .even
-
- .globl SYM (_fpCCR)
-
-SYM (_fpCCR):
-__exception_bits:
- .word 0
-__trap_enable_bits:
- .word 0
-__sticky_bits:
- .word 0
-__rounding_mode:
- .word ROUND_TO_NEAREST
-__format:
- .word NIL
-__last_operation:
- .word NOOP
-__operand1:
- .long 0
- .long 0
-__operand2:
- .long 0
- .long 0
-
-| Offsets:
-EBITS = __exception_bits - SYM (_fpCCR)
-TRAPE = __trap_enable_bits - SYM (_fpCCR)
-STICK = __sticky_bits - SYM (_fpCCR)
-ROUND = __rounding_mode - SYM (_fpCCR)
-FORMT = __format - SYM (_fpCCR)
-LASTO = __last_operation - SYM (_fpCCR)
-OPER1 = __operand1 - SYM (_fpCCR)
-OPER2 = __operand2 - SYM (_fpCCR)
-
-| The following exception types are supported:
-INEXACT_RESULT = 0x0001
-UNDERFLOW = 0x0002
-OVERFLOW = 0x0004
-DIVIDE_BY_ZERO = 0x0008
-INVALID_OPERATION = 0x0010
-
-| The allowed rounding modes are:
-UNKNOWN = -1
-ROUND_TO_NEAREST = 0 | round result to nearest representable value
-ROUND_TO_ZERO = 1 | round result towards zero
-ROUND_TO_PLUS = 2 | round result towards plus infinity
-ROUND_TO_MINUS = 3 | round result towards minus infinity
-
-| The allowed values of format are:
-NIL = 0
-SINGLE_FLOAT = 1
-DOUBLE_FLOAT = 2
-LONG_FLOAT = 3
-
-| The allowed values for the last operation are:
-NOOP = 0
-ADD = 1
-MULTIPLY = 2
-DIVIDE = 3
-NEGATE = 4
-COMPARE = 5
-EXTENDSFDF = 6
-TRUNCDFSF = 7
-
-|=============================================================================
-| __clear_sticky_bits
-|=============================================================================
-
-| The sticky bits are normally not cleared (thus the name), whereas the
-| exception type and exception value reflect the last computation.
-| This routine is provided to clear them (you can also write to _fpCCR,
-| since it is globally visible).
-
- .globl SYM (__clear_sticky_bit)
-
- .text
- .even
-
-| void __clear_sticky_bits(void);
-SYM (__clear_sticky_bit):
- PICLEA SYM (_fpCCR),a0
-#ifndef __mcoldfire__
- movew IMM (0),a0@(STICK)
-#else
- clr.w a0@(STICK)
-#endif
- rts
-
-|=============================================================================
-| $_exception_handler
-|=============================================================================
-
- .globl $_exception_handler
-
- .text
- .even
-
-| This is the common exit point if an exception occurs.
-| NOTE: it is NOT callable from C!
-| It expects the exception type in d7, the format (SINGLE_FLOAT,
-| DOUBLE_FLOAT or LONG_FLOAT) in d6, and the last operation code in d5.
-| It sets the corresponding exception and sticky bits, and the format.
-| Depending on the format if fills the corresponding slots for the
-| operands which produced the exception (all this information is provided
-| so if you write your own exception handlers you have enough information
-| to deal with the problem).
-| Then checks to see if the corresponding exception is trap-enabled,
-| in which case it pushes the address of _fpCCR and traps through
-| trap FPTRAP (15 for the moment).
-
-FPTRAP = 15
-
-$_exception_handler:
- PICLEA SYM (_fpCCR),a0
- movew d7,a0@(EBITS) | set __exception_bits
-#ifndef __mcoldfire__
- orw d7,a0@(STICK) | and __sticky_bits
-#else
- movew a0@(STICK),d4
- orl d7,d4
- movew d4,a0@(STICK)
-#endif
- movew d6,a0@(FORMT) | and __format
- movew d5,a0@(LASTO) | and __last_operation
-
-| Now put the operands in place:
-#ifndef __mcoldfire__
- cmpw IMM (SINGLE_FLOAT),d6
-#else
- cmpl IMM (SINGLE_FLOAT),d6
-#endif
- beq 1f
- movel a6@(8),a0@(OPER1)
- movel a6@(12),a0@(OPER1+4)
- movel a6@(16),a0@(OPER2)
- movel a6@(20),a0@(OPER2+4)
- bra 2f
-1: movel a6@(8),a0@(OPER1)
- movel a6@(12),a0@(OPER2)
-2:
-| And check whether the exception is trap-enabled:
-#ifndef __mcoldfire__
- andw a0@(TRAPE),d7 | is exception trap-enabled?
-#else
- clrl d6
- movew a0@(TRAPE),d6
- andl d6,d7
-#endif
- beq 1f | no, exit
- PICPEA SYM (_fpCCR),a1 | yes, push address of _fpCCR
- trap IMM (FPTRAP) | and trap
-#ifndef __mcoldfire__
-1: moveml sp@+,d2-d7 | restore data registers
-#else
-1: moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 | and return
- rts
-#endif /* L_floatex */
#ifdef L_mulsi3
.text
@@ -555,3420 +348,3 @@
#endif /* L_modsi3 */
-#ifdef L_double
-
- .globl SYM (_fpCCR)
- .globl $_exception_handler
-
-QUIET_NaN = 0xffffffff
-
-D_MAX_EXP = 0x07ff
-D_BIAS = 1022
-DBL_MAX_EXP = D_MAX_EXP - D_BIAS
-DBL_MIN_EXP = 1 - D_BIAS
-DBL_MANT_DIG = 53
-
-INEXACT_RESULT = 0x0001
-UNDERFLOW = 0x0002
-OVERFLOW = 0x0004
-DIVIDE_BY_ZERO = 0x0008
-INVALID_OPERATION = 0x0010
-
-DOUBLE_FLOAT = 2
-
-NOOP = 0
-ADD = 1
-MULTIPLY = 2
-DIVIDE = 3
-NEGATE = 4
-COMPARE = 5
-EXTENDSFDF = 6
-TRUNCDFSF = 7
-
-UNKNOWN = -1
-ROUND_TO_NEAREST = 0 | round result to nearest representable value
-ROUND_TO_ZERO = 1 | round result towards zero
-ROUND_TO_PLUS = 2 | round result towards plus infinity
-ROUND_TO_MINUS = 3 | round result towards minus infinity
-
-| Entry points:
-
- .globl SYM (__adddf3)
- .globl SYM (__subdf3)
- .globl SYM (__muldf3)
- .globl SYM (__divdf3)
- .globl SYM (__negdf2)
- .globl SYM (__cmpdf2)
-
- .text
- .even
-
-| These are common routines to return and signal exceptions.
-
-Ld$den:
-| Return and signal a denormalized number
- orl d7,d0
- movew IMM (INEXACT_RESULT+UNDERFLOW),d7
- moveq IMM (DOUBLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Ld$infty:
-Ld$overflow:
-| Return a properly signed INFINITY and set the exception flags
- movel IMM (0x7ff00000),d0
- movel IMM (0),d1
- orl d7,d0
- movew IMM (INEXACT_RESULT+OVERFLOW),d7
- moveq IMM (DOUBLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Ld$underflow:
-| Return 0 and set the exception flags
- movel IMM (0),d0
- movel d0,d1
- movew IMM (INEXACT_RESULT+UNDERFLOW),d7
- moveq IMM (DOUBLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Ld$inop:
-| Return a quiet NaN and set the exception flags
- movel IMM (QUIET_NaN),d0
- movel d0,d1
- movew IMM (INEXACT_RESULT+INVALID_OPERATION),d7
- moveq IMM (DOUBLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Ld$div$0:
-| Return a properly signed INFINITY and set the exception flags
- movel IMM (0x7ff00000),d0
- movel IMM (0),d1
- orl d7,d0
- movew IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
- moveq IMM (DOUBLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-|=============================================================================
-|=============================================================================
-| double precision routines
-|=============================================================================
-|=============================================================================
-
-| A double precision floating point number (double) has the format:
-|
-| struct _double {
-| unsigned int sign : 1; /* sign bit */
-| unsigned int exponent : 11; /* exponent, shifted by 126 */
-| unsigned int fraction : 52; /* fraction */
-| } double;
-|
-| Thus sizeof(double) = 8 (64 bits).
-|
-| All the routines are callable from C programs, and return the result
-| in the register pair d0-d1. They also preserve all registers except
-| d0-d1 and a0-a1.
-
-|=============================================================================
-| __subdf3
-|=============================================================================
-
-| double __subdf3(double, double);
-SYM (__subdf3):
- bchg IMM (31),sp@(12) | change sign of second operand
- | and fall through, so we always add
-|=============================================================================
-| __adddf3
-|=============================================================================
-
-| double __adddf3(double, double);
-SYM (__adddf3):
-#ifndef __mcoldfire__
- link a6,IMM (0) | everything will be done in registers
- moveml d2-d7,sp@- | save all data registers and a2 (but d0-d1)
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movel a6@(8),d0 | get first operand
- movel a6@(12),d1 |
- movel a6@(16),d2 | get second operand
- movel a6@(20),d3 |
-
- movel d0,d7 | get d0's sign bit in d7 '
- addl d1,d1 | check and clear sign bit of a, and gain one
- addxl d0,d0 | bit of extra precision
- beq Ladddf$b | if zero return second operand
-
- movel d2,d6 | save sign in d6
- addl d3,d3 | get rid of sign bit and gain one bit of
- addxl d2,d2 | extra precision
- beq Ladddf$a | if zero return first operand
-
- andl IMM (0x80000000),d7 | isolate a's sign bit '
- swap d6 | and also b's sign bit '
-#ifndef __mcoldfire__
- andw IMM (0x8000),d6 |
- orw d6,d7 | and combine them into d7, so that a's sign '
- | bit is in the high word and b's is in the '
- | low word, so d6 is free to be used
-#else
- andl IMM (0x8000),d6
- orl d6,d7
-#endif
- movel d7,a0 | now save d7 into a0, so d7 is free to
- | be used also
-
-| Get the exponents and check for denormalized and/or infinity.
-
- movel IMM (0x001fffff),d6 | mask for the fraction
- movel IMM (0x00200000),d7 | mask to put hidden bit back
-
- movel d0,d4 |
- andl d6,d0 | get fraction in d0
- notl d6 | make d6 into mask for the exponent
- andl d6,d4 | get exponent in d4
- beq Ladddf$a$den | branch if a is denormalized
- cmpl d6,d4 | check for INFINITY or NaN
- beq Ladddf$nf |
- orl d7,d0 | and put hidden bit back
-Ladddf$1:
- swap d4 | shift right exponent so that it starts
-#ifndef __mcoldfire__
- lsrw IMM (5),d4 | in bit 0 and not bit 20
-#else
- lsrl IMM (5),d4 | in bit 0 and not bit 20
-#endif
-| Now we have a's exponent in d4 and fraction in d0-d1 '
- movel d2,d5 | save b to get exponent
- andl d6,d5 | get exponent in d5
- beq Ladddf$b$den | branch if b is denormalized
- cmpl d6,d5 | check for INFINITY or NaN
- beq Ladddf$nf
- notl d6 | make d6 into mask for the fraction again
- andl d6,d2 | and get fraction in d2
- orl d7,d2 | and put hidden bit back
-Ladddf$2:
- swap d5 | shift right exponent so that it starts
-#ifndef __mcoldfire__
- lsrw IMM (5),d5 | in bit 0 and not bit 20
-#else
- lsrl IMM (5),d5 | in bit 0 and not bit 20
-#endif
-
-| Now we have b's exponent in d5 and fraction in d2-d3. '
-
-| The situation now is as follows: the signs are combined in a0, the
-| numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
-| and d5 (b). To do the rounding correctly we need to keep all the
-| bits until the end, so we need to use d0-d1-d2-d3 for the first number
-| and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
-| exponents in a2-a3.
-
-#ifndef __mcoldfire__
- moveml a2-a3,sp@- | save the address registers
-#else
- movel a2,sp@-
- movel a3,sp@-
- movel a4,sp@-
-#endif
-
- movel d4,a2 | save the exponents
- movel d5,a3 |
-
- movel IMM (0),d7 | and move the numbers around
- movel d7,d6 |
- movel d3,d5 |
- movel d2,d4 |
- movel d7,d3 |
- movel d7,d2 |
-
-| Here we shift the numbers until the exponents are the same, and put
-| the largest exponent in a2.
-#ifndef __mcoldfire__
- exg d4,a2 | get exponents back
- exg d5,a3 |
- cmpw d4,d5 | compare the exponents
-#else
- movel d4,a4 | get exponents back
- movel a2,d4
- movel a4,a2
- movel d5,a4
- movel a3,d5
- movel a4,a3
- cmpl d4,d5 | compare the exponents
-#endif
- beq Ladddf$3 | if equal don't shift '
- bhi 9f | branch if second exponent is higher
-
-| Here we have a's exponent larger than b's, so we have to shift b. We do
-| this by using as counter d2:
-1: movew d4,d2 | move largest exponent to d2
-#ifndef __mcoldfire__
- subw d5,d2 | and subtract second exponent
- exg d4,a2 | get back the longs we saved
- exg d5,a3 |
-#else
- subl d5,d2 | and subtract second exponent
- movel d4,a4 | get back the longs we saved
- movel a2,d4
- movel a4,a2
- movel d5,a4
- movel a3,d5
- movel a4,a3
-#endif
-| if difference is too large we don't shift (actually, we can just exit) '
-#ifndef __mcoldfire__
- cmpw IMM (DBL_MANT_DIG+2),d2
-#else
- cmpl IMM (DBL_MANT_DIG+2),d2
-#endif
- bge Ladddf$b$small
-#ifndef __mcoldfire__
- cmpw IMM (32),d2 | if difference >= 32, shift by longs
-#else
- cmpl IMM (32),d2 | if difference >= 32, shift by longs
-#endif
- bge 5f
-2:
-#ifndef __mcoldfire__
- cmpw IMM (16),d2 | if difference >= 16, shift by words
-#else
- cmpl IMM (16),d2 | if difference >= 16, shift by words
-#endif
- bge 6f
- bra 3f | enter dbra loop
-
-4:
-#ifndef __mcoldfire__
- lsrl IMM (1),d4
- roxrl IMM (1),d5
- roxrl IMM (1),d6
- roxrl IMM (1),d7
-#else
- lsrl IMM (1),d7
- btst IMM (0),d6
- beq 10f
- bset IMM (31),d7
-10: lsrl IMM (1),d6
- btst IMM (0),d5
- beq 11f
- bset IMM (31),d6
-11: lsrl IMM (1),d5
- btst IMM (0),d4
- beq 12f
- bset IMM (31),d5
-12: lsrl IMM (1),d4
-#endif
-3:
-#ifndef __mcoldfire__
- dbra d2,4b
-#else
- subql IMM (1),d2
- bpl 4b
-#endif
- movel IMM (0),d2
- movel d2,d3
- bra Ladddf$4
-5:
- movel d6,d7
- movel d5,d6
- movel d4,d5
- movel IMM (0),d4
-#ifndef __mcoldfire__
- subw IMM (32),d2
-#else
- subl IMM (32),d2
-#endif
- bra 2b
-6:
- movew d6,d7
- swap d7
- movew d5,d6
- swap d6
- movew d4,d5
- swap d5
- movew IMM (0),d4
- swap d4
-#ifndef __mcoldfire__
- subw IMM (16),d2
-#else
- subl IMM (16),d2
-#endif
- bra 3b
-
-9:
-#ifndef __mcoldfire__
- exg d4,d5
- movew d4,d6
- subw d5,d6 | keep d5 (largest exponent) in d4
- exg d4,a2
- exg d5,a3
-#else
- movel d5,d6
- movel d4,d5
- movel d6,d4
- subl d5,d6
- movel d4,a4
- movel a2,d4
- movel a4,a2
- movel d5,a4
- movel a3,d5
- movel a4,a3
-#endif
-| if difference is too large we don't shift (actually, we can just exit) '
-#ifndef __mcoldfire__
- cmpw IMM (DBL_MANT_DIG+2),d6
-#else
- cmpl IMM (DBL_MANT_DIG+2),d6
-#endif
- bge Ladddf$a$small
-#ifndef __mcoldfire__
- cmpw IMM (32),d6 | if difference >= 32, shift by longs
-#else
- cmpl IMM (32),d6 | if difference >= 32, shift by longs
-#endif
- bge 5f
-2:
-#ifndef __mcoldfire__
- cmpw IMM (16),d6 | if difference >= 16, shift by words
-#else
- cmpl IMM (16),d6 | if difference >= 16, shift by words
-#endif
- bge 6f
- bra 3f | enter dbra loop
-
-4:
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- roxrl IMM (1),d2
- roxrl IMM (1),d3
-#else
- lsrl IMM (1),d3
- btst IMM (0),d2
- beq 10f
- bset IMM (31),d3
-10: lsrl IMM (1),d2
- btst IMM (0),d1
- beq 11f
- bset IMM (31),d2
-11: lsrl IMM (1),d1
- btst IMM (0),d0
- beq 12f
- bset IMM (31),d1
-12: lsrl IMM (1),d0
-#endif
-3:
-#ifndef __mcoldfire__
- dbra d6,4b
-#else
- subql IMM (1),d6
- bpl 4b
-#endif
- movel IMM (0),d7
- movel d7,d6
- bra Ladddf$4
-5:
- movel d2,d3
- movel d1,d2
- movel d0,d1
- movel IMM (0),d0
-#ifndef __mcoldfire__
- subw IMM (32),d6
-#else
- subl IMM (32),d6
-#endif
- bra 2b
-6:
- movew d2,d3
- swap d3
- movew d1,d2
- swap d2
- movew d0,d1
- swap d1
- movew IMM (0),d0
- swap d0
-#ifndef __mcoldfire__
- subw IMM (16),d6
-#else
- subl IMM (16),d6
-#endif
- bra 3b
-Ladddf$3:
-#ifndef __mcoldfire__
- exg d4,a2
- exg d5,a3
-#else
- movel d4,a4
- movel a2,d4
- movel a4,a2
- movel d5,a4
- movel a3,d5
- movel a4,a3
-#endif
-Ladddf$4:
-| Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
-| the signs in a4.
-
-| Here we have to decide whether to add or subtract the numbers:
-#ifndef __mcoldfire__
- exg d7,a0 | get the signs
- exg d6,a3 | a3 is free to be used
-#else
- movel d7,a4
- movel a0,d7
- movel a4,a0
- movel d6,a4
- movel a3,d6
- movel a4,a3
-#endif
- movel d7,d6 |
- movew IMM (0),d7 | get a's sign in d7 '
- swap d6 |
- movew IMM (0),d6 | and b's sign in d6 '
- eorl d7,d6 | compare the signs
- bmi Lsubdf$0 | if the signs are different we have
- | to subtract
-#ifndef __mcoldfire__
- exg d7,a0 | else we add the numbers
- exg d6,a3 |
-#else
- movel d7,a4
- movel a0,d7
- movel a4,a0
- movel d6,a4
- movel a3,d6
- movel a4,a3
-#endif
- addl d7,d3 |
- addxl d6,d2 |
- addxl d5,d1 |
- addxl d4,d0 |
-
- movel a2,d4 | return exponent to d4
- movel a0,d7 |
- andl IMM (0x80000000),d7 | d7 now has the sign
-
-#ifndef __mcoldfire__
- moveml sp@+,a2-a3
-#else
- movel sp@+,a4
- movel sp@+,a3
- movel sp@+,a2
-#endif
-
-| Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
-| the case of denormalized numbers in the rounding routine itself).
-| As in the addition (not in the subtraction!) we could have set
-| one more bit we check this:
- btst IMM (DBL_MANT_DIG+1),d0
- beq 1f
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- roxrl IMM (1),d2
- roxrl IMM (1),d3
- addw IMM (1),d4
-#else
- lsrl IMM (1),d3
- btst IMM (0),d2
- beq 10f
- bset IMM (31),d3
-10: lsrl IMM (1),d2
- btst IMM (0),d1
- beq 11f
- bset IMM (31),d2
-11: lsrl IMM (1),d1
- btst IMM (0),d0
- beq 12f
- bset IMM (31),d1
-12: lsrl IMM (1),d0
- addl IMM (1),d4
-#endif
-1:
- lea pc@(Ladddf$5),a0 | to return from rounding routine
- PICLEA SYM (_fpCCR),a1 | check the rounding mode
-#ifdef __mcoldfire__
- clrl d6
-#endif
- movew a1@(6),d6 | rounding mode in d6
- beq Lround$to$nearest
-#ifndef __mcoldfire__
- cmpw IMM (ROUND_TO_PLUS),d6
-#else
- cmpl IMM (ROUND_TO_PLUS),d6
-#endif
- bhi Lround$to$minus
- blt Lround$to$zero
- bra Lround$to$plus
-Ladddf$5:
-| Put back the exponent and check for overflow
-#ifndef __mcoldfire__
- cmpw IMM (0x7ff),d4 | is the exponent big?
-#else
- cmpl IMM (0x7ff),d4 | is the exponent big?
-#endif
- bge 1f
- bclr IMM (DBL_MANT_DIG-1),d0
-#ifndef __mcoldfire__
- lslw IMM (4),d4 | put exponent back into position
-#else
- lsll IMM (4),d4 | put exponent back into position
-#endif
- swap d0 |
-#ifndef __mcoldfire__
- orw d4,d0 |
-#else
- orl d4,d0 |
-#endif
- swap d0 |
- bra Ladddf$ret
-1:
- movew IMM (ADD),d5
- bra Ld$overflow
-
-Lsubdf$0:
-| Here we do the subtraction.
-#ifndef __mcoldfire__
- exg d7,a0 | put sign back in a0
- exg d6,a3 |
-#else
- movel d7,a4
- movel a0,d7
- movel a4,a0
- movel d6,a4
- movel a3,d6
- movel a4,a3
-#endif
- subl d7,d3 |
- subxl d6,d2 |
- subxl d5,d1 |
- subxl d4,d0 |
- beq Ladddf$ret$1 | if zero just exit
- bpl 1f | if positive skip the following
- movel a0,d7 |
- bchg IMM (31),d7 | change sign bit in d7
- movel d7,a0 |
- negl d3 |
- negxl d2 |
- negxl d1 | and negate result
- negxl d0 |
-1:
- movel a2,d4 | return exponent to d4
- movel a0,d7
- andl IMM (0x80000000),d7 | isolate sign bit
-#ifndef __mcoldfire__
- moveml sp@+,a2-a3 |
-#else
- movel sp@+,a4
- movel sp@+,a3
- movel sp@+,a2
-#endif
-
-| Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
-| the case of denormalized numbers in the rounding routine itself).
-| As in the addition (not in the subtraction!) we could have set
-| one more bit we check this:
- btst IMM (DBL_MANT_DIG+1),d0
- beq 1f
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- roxrl IMM (1),d2
- roxrl IMM (1),d3
- addw IMM (1),d4
-#else
- lsrl IMM (1),d3
- btst IMM (0),d2
- beq 10f
- bset IMM (31),d3
-10: lsrl IMM (1),d2
- btst IMM (0),d1
- beq 11f
- bset IMM (31),d2
-11: lsrl IMM (1),d1
- btst IMM (0),d0
- beq 12f
- bset IMM (31),d1
-12: lsrl IMM (1),d0
- addl IMM (1),d4
-#endif
-1:
- lea pc@(Lsubdf$1),a0 | to return from rounding routine
- PICLEA SYM (_fpCCR),a1 | check the rounding mode
-#ifdef __mcoldfire__
- clrl d6
-#endif
- movew a1@(6),d6 | rounding mode in d6
- beq Lround$to$nearest
-#ifndef __mcoldfire__
- cmpw IMM (ROUND_TO_PLUS),d6
-#else
- cmpl IMM (ROUND_TO_PLUS),d6
-#endif
- bhi Lround$to$minus
- blt Lround$to$zero
- bra Lround$to$plus
-Lsubdf$1:
-| Put back the exponent and sign (we don't have overflow). '
- bclr IMM (DBL_MANT_DIG-1),d0
-#ifndef __mcoldfire__
- lslw IMM (4),d4 | put exponent back into position
-#else
- lsll IMM (4),d4 | put exponent back into position
-#endif
- swap d0 |
-#ifndef __mcoldfire__
- orw d4,d0 |
-#else
- orl d4,d0 |
-#endif
- swap d0 |
- bra Ladddf$ret
-
-| If one of the numbers was too small (difference of exponents >=
-| DBL_MANT_DIG+1) we return the other (and now we don't have to '
-| check for finiteness or zero).
-Ladddf$a$small:
-#ifndef __mcoldfire__
- moveml sp@+,a2-a3
-#else
- movel sp@+,a4
- movel sp@+,a3
- movel sp@+,a2
-#endif
- movel a6@(16),d0
- movel a6@(20),d1
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | restore data registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 | and return
- rts
-
-Ladddf$b$small:
-#ifndef __mcoldfire__
- moveml sp@+,a2-a3
-#else
- movel sp@+,a4
- movel sp@+,a3
- movel sp@+,a2
-#endif
- movel a6@(8),d0
- movel a6@(12),d1
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | restore data registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 | and return
- rts
-
-Ladddf$a$den:
- movel d7,d4 | d7 contains 0x00200000
- bra Ladddf$1
-
-Ladddf$b$den:
- movel d7,d5 | d7 contains 0x00200000
- notl d6
- bra Ladddf$2
-
-Ladddf$b:
-| Return b (if a is zero)
- movel d2,d0
- movel d3,d1
- bra 1f
-Ladddf$a:
- movel a6@(8),d0
- movel a6@(12),d1
-1:
- movew IMM (ADD),d5
-| Check for NaN and +/-INFINITY.
- movel d0,d7 |
- andl IMM (0x80000000),d7 |
- bclr IMM (31),d0 |
- cmpl IMM (0x7ff00000),d0 |
- bge 2f |
- movel d0,d0 | check for zero, since we don't '
- bne Ladddf$ret | want to return -0 by mistake
- bclr IMM (31),d7 |
- bra Ladddf$ret |
-2:
- andl IMM (0x000fffff),d0 | check for NaN (nonzero fraction)
- orl d1,d0 |
- bne Ld$inop |
- bra Ld$infty |
-
-Ladddf$ret$1:
-#ifndef __mcoldfire__
- moveml sp@+,a2-a3 | restore regs and exit
-#else
- movel sp@+,a4
- movel sp@+,a3
- movel sp@+,a2
-#endif
-
-Ladddf$ret:
-| Normal exit.
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
- orl d7,d0 | put sign bit back
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-
-Ladddf$ret$den:
-| Return a denormalized number.
-#ifndef __mcoldfire__
- lsrl IMM (1),d0 | shift right once more
- roxrl IMM (1),d1 |
-#else
- lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
-#endif
- bra Ladddf$ret
-
-Ladddf$nf:
- movew IMM (ADD),d5
-| This could be faster but it is not worth the effort, since it is not
-| executed very often. We sacrifice speed for clarity here.
- movel a6@(8),d0 | get the numbers back (remember that we
- movel a6@(12),d1 | did some processing already)
- movel a6@(16),d2 |
- movel a6@(20),d3 |
- movel IMM (0x7ff00000),d4 | useful constant (INFINITY)
- movel d0,d7 | save sign bits
- movel d2,d6 |
- bclr IMM (31),d0 | clear sign bits
- bclr IMM (31),d2 |
-| We know that one of them is either NaN of +/-INFINITY
-| Check for NaN (if either one is NaN return NaN)
- cmpl d4,d0 | check first a (d0)
- bhi Ld$inop | if d0 > 0x7ff00000 or equal and
- bne 2f
- tstl d1 | d1 > 0, a is NaN
- bne Ld$inop |
-2: cmpl d4,d2 | check now b (d1)
- bhi Ld$inop |
- bne 3f
- tstl d3 |
- bne Ld$inop |
-3:
-| Now comes the check for +/-INFINITY. We know that both are (maybe not
-| finite) numbers, but we have to check if both are infinite whether we
-| are adding or subtracting them.
- eorl d7,d6 | to check sign bits
- bmi 1f
- andl IMM (0x80000000),d7 | get (common) sign bit
- bra Ld$infty
-1:
-| We know one (or both) are infinite, so we test for equality between the
-| two numbers (if they are equal they have to be infinite both, so we
-| return NaN).
- cmpl d2,d0 | are both infinite?
- bne 1f | if d0 <> d2 they are not equal
- cmpl d3,d1 | if d0 == d2 test d3 and d1
- beq Ld$inop | if equal return NaN
-1:
- andl IMM (0x80000000),d7 | get a's sign bit '
- cmpl d4,d0 | test now for infinity
- beq Ld$infty | if a is INFINITY return with this sign
- bchg IMM (31),d7 | else we know b is INFINITY and has
- bra Ld$infty | the opposite sign
-
-|=============================================================================
-| __muldf3
-|=============================================================================
-
-| double __muldf3(double, double);
-SYM (__muldf3):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@-
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movel a6@(8),d0 | get a into d0-d1
- movel a6@(12),d1 |
- movel a6@(16),d2 | and b into d2-d3
- movel a6@(20),d3 |
- movel d0,d7 | d7 will hold the sign of the product
- eorl d2,d7 |
- andl IMM (0x80000000),d7 |
- movel d7,a0 | save sign bit into a0
- movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
- movel d7,d6 | another (mask for fraction)
- notl d6 |
- bclr IMM (31),d0 | get rid of a's sign bit '
- movel d0,d4 |
- orl d1,d4 |
- beq Lmuldf$a$0 | branch if a is zero
- movel d0,d4 |
- bclr IMM (31),d2 | get rid of b's sign bit '
- movel d2,d5 |
- orl d3,d5 |
- beq Lmuldf$b$0 | branch if b is zero
- movel d2,d5 |
- cmpl d7,d0 | is a big?
- bhi Lmuldf$inop | if a is NaN return NaN
- beq Lmuldf$a$nf | we still have to check d1 and b ...
- cmpl d7,d2 | now compare b with INFINITY
- bhi Lmuldf$inop | is b NaN?
- beq Lmuldf$b$nf | we still have to check d3 ...
-| Here we have both numbers finite and nonzero (and with no sign bit).
-| Now we get the exponents into d4 and d5.
- andl d7,d4 | isolate exponent in d4
- beq Lmuldf$a$den | if exponent zero, have denormalized
- andl d6,d0 | isolate fraction
- orl IMM (0x00100000),d0 | and put hidden bit back
- swap d4 | I like exponents in the first byte
-#ifndef __mcoldfire__
- lsrw IMM (4),d4 |
-#else
- lsrl IMM (4),d4 |
-#endif
-Lmuldf$1:
- andl d7,d5 |
- beq Lmuldf$b$den |
- andl d6,d2 |
- orl IMM (0x00100000),d2 | and put hidden bit back
- swap d5 |
-#ifndef __mcoldfire__
- lsrw IMM (4),d5 |
-#else
- lsrl IMM (4),d5 |
-#endif
-Lmuldf$2: |
-#ifndef __mcoldfire__
- addw d5,d4 | add exponents
- subw IMM (D_BIAS+1),d4 | and subtract bias (plus one)
-#else
- addl d5,d4 | add exponents
- subl IMM (D_BIAS+1),d4 | and subtract bias (plus one)
-#endif
-
-| We are now ready to do the multiplication. The situation is as follows:
-| both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were
-| denormalized to start with!), which means that in the product bit 104
-| (which will correspond to bit 8 of the fourth long) is set.
-
-| Here we have to do the product.
-| To do it we have to juggle the registers back and forth, as there are not
-| enough to keep everything in them. So we use the address registers to keep
-| some intermediate data.
-
-#ifndef __mcoldfire__
- moveml a2-a3,sp@- | save a2 and a3 for temporary use
-#else
- movel a2,sp@-
- movel a3,sp@-
- movel a4,sp@-
-#endif
- movel IMM (0),a2 | a2 is a null register
- movel d4,a3 | and a3 will preserve the exponent
-
-| First, shift d2-d3 so bit 20 becomes bit 31:
-#ifndef __mcoldfire__
- rorl IMM (5),d2 | rotate d2 5 places right
- swap d2 | and swap it
- rorl IMM (5),d3 | do the same thing with d3
- swap d3 |
- movew d3,d6 | get the rightmost 11 bits of d3
- andw IMM (0x07ff),d6 |
- orw d6,d2 | and put them into d2
- andw IMM (0xf800),d3 | clear those bits in d3
-#else
- moveq IMM (11),d7 | left shift d2 11 bits
- lsll d7,d2
- movel d3,d6 | get a copy of d3
- lsll d7,d3 | left shift d3 11 bits
- andl IMM (0xffe00000),d6 | get the top 11 bits of d3
- moveq IMM (21),d7 | right shift them 21 bits
- lsrl d7,d6
- orl d6,d2 | stick them at the end of d2
-#endif
-
- movel d2,d6 | move b into d6-d7
- movel d3,d7 | move a into d4-d5
- movel d0,d4 | and clear d0-d1-d2-d3 (to put result)
- movel d1,d5 |
- movel IMM (0),d3 |
- movel d3,d2 |
- movel d3,d1 |
- movel d3,d0 |
-
-| We use a1 as counter:
- movel IMM (DBL_MANT_DIG-1),a1
-#ifndef __mcoldfire__
- exg d7,a1
-#else
- movel d7,a4
- movel a1,d7
- movel a4,a1
-#endif
-
-1:
-#ifndef __mcoldfire__
- exg d7,a1 | put counter back in a1
-#else
- movel d7,a4
- movel a1,d7
- movel a4,a1
-#endif
- addl d3,d3 | shift sum once left
- addxl d2,d2 |
- addxl d1,d1 |
- addxl d0,d0 |
- addl d7,d7 |
- addxl d6,d6 |
- bcc 2f | if bit clear skip the following
-#ifndef __mcoldfire__
- exg d7,a2 |
-#else
- movel d7,a4
- movel a2,d7
- movel a4,a2
-#endif
- addl d5,d3 | else add a to the sum
- addxl d4,d2 |
- addxl d7,d1 |
- addxl d7,d0 |
-#ifndef __mcoldfire__
- exg d7,a2 |
-#else
- movel d7,a4
- movel a2,d7
- movel a4,a2
-#endif
-2:
-#ifndef __mcoldfire__
- exg d7,a1 | put counter in d7
- dbf d7,1b | decrement and branch
-#else
- movel d7,a4
- movel a1,d7
- movel a4,a1
- subql IMM (1),d7
- bpl 1b
-#endif
-
- movel a3,d4 | restore exponent
-#ifndef __mcoldfire__
- moveml sp@+,a2-a3
-#else
- movel sp@+,a4
- movel sp@+,a3
- movel sp@+,a2
-#endif
-
-| Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The
-| first thing to do now is to normalize it so bit 8 becomes bit
-| DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
- swap d0
- swap d1
- movew d1,d0
- swap d2
- movew d2,d1
- swap d3
- movew d3,d2
- movew IMM (0),d3
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- roxrl IMM (1),d2
- roxrl IMM (1),d3
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- roxrl IMM (1),d2
- roxrl IMM (1),d3
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- roxrl IMM (1),d2
- roxrl IMM (1),d3
-#else
- moveq IMM (29),d6
- lsrl IMM (3),d3
- movel d2,d7
- lsll d6,d7
- orl d7,d3
- lsrl IMM (3),d2
- movel d1,d7
- lsll d6,d7
- orl d7,d2
- lsrl IMM (3),d1
- movel d0,d7
- lsll d6,d7
- orl d7,d1
- lsrl IMM (3),d0
-#endif
-
-| Now round, check for over- and underflow, and exit.
- movel a0,d7 | get sign bit back into d7
- movew IMM (MULTIPLY),d5
-
- btst IMM (DBL_MANT_DIG+1-32),d0
- beq Lround$exit
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- addw IMM (1),d4
-#else
- lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
- addl IMM (1),d4
-#endif
- bra Lround$exit
-
-Lmuldf$inop:
- movew IMM (MULTIPLY),d5
- bra Ld$inop
-
-Lmuldf$b$nf:
- movew IMM (MULTIPLY),d5
- movel a0,d7 | get sign bit back into d7
- tstl d3 | we know d2 == 0x7ff00000, so check d3
- bne Ld$inop | if d3 <> 0 b is NaN
- bra Ld$overflow | else we have overflow (since a is finite)
-
-Lmuldf$a$nf:
- movew IMM (MULTIPLY),d5
- movel a0,d7 | get sign bit back into d7
- tstl d1 | we know d0 == 0x7ff00000, so check d1
- bne Ld$inop | if d1 <> 0 a is NaN
- bra Ld$overflow | else signal overflow
-
-| If either number is zero return zero, unless the other is +/-INFINITY or
-| NaN, in which case we return NaN.
-Lmuldf$b$0:
- movew IMM (MULTIPLY),d5
-#ifndef __mcoldfire__
- exg d2,d0 | put b (==0) into d0-d1
- exg d3,d1 | and a (with sign bit cleared) into d2-d3
-#else
- movel d2,d7
- movel d0,d2
- movel d7,d0
- movel d3,d7
- movel d1,d3
- movel d7,d1
-#endif
- bra 1f
-Lmuldf$a$0:
- movel a6@(16),d2 | put b into d2-d3 again
- movel a6@(20),d3 |
- bclr IMM (31),d2 | clear sign bit
-1: cmpl IMM (0x7ff00000),d2 | check for non-finiteness
- bge Ld$inop | in case NaN or +/-INFINITY return NaN
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-
-| If a number is denormalized we put an exponent of 1 but do not put the
-| hidden bit back into the fraction; instead we shift left until bit 21
-| (the hidden bit) is set, adjusting the exponent accordingly. We do this
-| to ensure that the product of the fractions is close to 1.
-Lmuldf$a$den:
- movel IMM (1),d4
- andl d6,d0
-1: addl d1,d1 | shift a left until bit 20 is set
- addxl d0,d0 |
-#ifndef __mcoldfire__
- subw IMM (1),d4 | and adjust exponent
-#else
- subl IMM (1),d4 | and adjust exponent
-#endif
- btst IMM (20),d0 |
- bne Lmuldf$1 |
- bra 1b
-
-Lmuldf$b$den:
- movel IMM (1),d5
- andl d6,d2
-1: addl d3,d3 | shift b left until bit 20 is set
- addxl d2,d2 |
-#ifndef __mcoldfire__
- subw IMM (1),d5 | and adjust exponent
-#else
- subql IMM (1),d5 | and adjust exponent
-#endif
- btst IMM (20),d2 |
- bne Lmuldf$2 |
- bra 1b
-
-
-|=============================================================================
-| __divdf3
-|=============================================================================
-
-| double __divdf3(double, double);
-SYM (__divdf3):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@-
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movel a6@(8),d0 | get a into d0-d1
- movel a6@(12),d1 |
- movel a6@(16),d2 | and b into d2-d3
- movel a6@(20),d3 |
- movel d0,d7 | d7 will hold the sign of the result
- eorl d2,d7 |
- andl IMM (0x80000000),d7
- movel d7,a0 | save sign into a0
- movel IMM (0x7ff00000),d7 | useful constant (+INFINITY)
- movel d7,d6 | another (mask for fraction)
- notl d6 |
- bclr IMM (31),d0 | get rid of a's sign bit '
- movel d0,d4 |
- orl d1,d4 |
- beq Ldivdf$a$0 | branch if a is zero
- movel d0,d4 |
- bclr IMM (31),d2 | get rid of b's sign bit '
- movel d2,d5 |
- orl d3,d5 |
- beq Ldivdf$b$0 | branch if b is zero
- movel d2,d5
- cmpl d7,d0 | is a big?
- bhi Ldivdf$inop | if a is NaN return NaN
- beq Ldivdf$a$nf | if d0 == 0x7ff00000 we check d1
- cmpl d7,d2 | now compare b with INFINITY
- bhi Ldivdf$inop | if b is NaN return NaN
- beq Ldivdf$b$nf | if d2 == 0x7ff00000 we check d3
-| Here we have both numbers finite and nonzero (and with no sign bit).
-| Now we get the exponents into d4 and d5 and normalize the numbers to
-| ensure that the ratio of the fractions is around 1. We do this by
-| making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
-| set, even if they were denormalized to start with.
-| Thus, the result will satisfy: 2 > result > 1/2.
- andl d7,d4 | and isolate exponent in d4
- beq Ldivdf$a$den | if exponent is zero we have a denormalized
- andl d6,d0 | and isolate fraction
- orl IMM (0x00100000),d0 | and put hidden bit back
- swap d4 | I like exponents in the first byte
-#ifndef __mcoldfire__
- lsrw IMM (4),d4 |
-#else
- lsrl IMM (4),d4 |
-#endif
-Ldivdf$1: |
- andl d7,d5 |
- beq Ldivdf$b$den |
- andl d6,d2 |
- orl IMM (0x00100000),d2
- swap d5 |
-#ifndef __mcoldfire__
- lsrw IMM (4),d5 |
-#else
- lsrl IMM (4),d5 |
-#endif
-Ldivdf$2: |
-#ifndef __mcoldfire__
- subw d5,d4 | subtract exponents
- addw IMM (D_BIAS),d4 | and add bias
-#else
- subl d5,d4 | subtract exponents
- addl IMM (D_BIAS),d4 | and add bias
-#endif
-
-| We are now ready to do the division. We have prepared things in such a way
-| that the ratio of the fractions will be less than 2 but greater than 1/2.
-| At this point the registers in use are:
-| d0-d1 hold a (first operand, bit DBL_MANT_DIG-32=0, bit
-| DBL_MANT_DIG-1-32=1)
-| d2-d3 hold b (second operand, bit DBL_MANT_DIG-32=1)
-| d4 holds the difference of the exponents, corrected by the bias
-| a0 holds the sign of the ratio
-
-| To do the rounding correctly we need to keep information about the
-| nonsignificant bits. One way to do this would be to do the division
-| using four registers; another is to use two registers (as originally
-| I did), but use a sticky bit to preserve information about the
-| fractional part. Note that we can keep that info in a1, which is not
-| used.
- movel IMM (0),d6 | d6-d7 will hold the result
- movel d6,d7 |
- movel IMM (0),a1 | and a1 will hold the sticky bit
-
- movel IMM (DBL_MANT_DIG-32+1),d5
-
-1: cmpl d0,d2 | is a < b?
- bhi 3f | if b > a skip the following
- beq 4f | if d0==d2 check d1 and d3
-2: subl d3,d1 |
- subxl d2,d0 | a <-- a - b
- bset d5,d6 | set the corresponding bit in d6
-3: addl d1,d1 | shift a by 1
- addxl d0,d0 |
-#ifndef __mcoldfire__
- dbra d5,1b | and branch back
-#else
- subql IMM (1), d5
- bpl 1b
-#endif
- bra 5f
-4: cmpl d1,d3 | here d0==d2, so check d1 and d3
- bhi 3b | if d1 > d2 skip the subtraction
- bra 2b | else go do it
-5:
-| Here we have to start setting the bits in the second long.
- movel IMM (31),d5 | again d5 is counter
-
-1: cmpl d0,d2 | is a < b?
- bhi 3f | if b > a skip the following
- beq 4f | if d0==d2 check d1 and d3
-2: subl d3,d1 |
- subxl d2,d0 | a <-- a - b
- bset d5,d7 | set the corresponding bit in d7
-3: addl d1,d1 | shift a by 1
- addxl d0,d0 |
-#ifndef __mcoldfire__
- dbra d5,1b | and branch back
-#else
- subql IMM (1), d5
- bpl 1b
-#endif
- bra 5f
-4: cmpl d1,d3 | here d0==d2, so check d1 and d3
- bhi 3b | if d1 > d2 skip the subtraction
- bra 2b | else go do it
-5:
-| Now go ahead checking until we hit a one, which we store in d2.
- movel IMM (DBL_MANT_DIG),d5
-1: cmpl d2,d0 | is a < b?
- bhi 4f | if b < a, exit
- beq 3f | if d0==d2 check d1 and d3
-2: addl d1,d1 | shift a by 1
- addxl d0,d0 |
-#ifndef __mcoldfire__
- dbra d5,1b | and branch back
-#else
- subql IMM (1), d5
- bpl 1b
-#endif
- movel IMM (0),d2 | here no sticky bit was found
- movel d2,d3
- bra 5f
-3: cmpl d1,d3 | here d0==d2, so check d1 and d3
- bhi 2b | if d1 > d2 go back
-4:
-| Here put the sticky bit in d2-d3 (in the position which actually corresponds
-| to it; if you don't do this the algorithm loses in some cases). '
- movel IMM (0),d2
- movel d2,d3
-#ifndef __mcoldfire__
- subw IMM (DBL_MANT_DIG),d5
- addw IMM (63),d5
- cmpw IMM (31),d5
-#else
- subl IMM (DBL_MANT_DIG),d5
- addl IMM (63),d5
- cmpl IMM (31),d5
-#endif
- bhi 2f
-1: bset d5,d3
- bra 5f
-#ifndef __mcoldfire__
- subw IMM (32),d5
-#else
- subl IMM (32),d5
-#endif
-2: bset d5,d2
-5:
-| Finally we are finished! Move the longs in the address registers to
-| their final destination:
- movel d6,d0
- movel d7,d1
- movel IMM (0),d3
-
-| Here we have finished the division, with the result in d0-d1-d2-d3, with
-| 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
-| If it is not, then definitely bit 21 is set. Normalize so bit 22 is
-| not set:
- btst IMM (DBL_MANT_DIG-32+1),d0
- beq 1f
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- roxrl IMM (1),d2
- roxrl IMM (1),d3
- addw IMM (1),d4
-#else
- lsrl IMM (1),d3
- btst IMM (0),d2
- beq 10f
- bset IMM (31),d3
-10: lsrl IMM (1),d2
- btst IMM (0),d1
- beq 11f
- bset IMM (31),d2
-11: lsrl IMM (1),d1
- btst IMM (0),d0
- beq 12f
- bset IMM (31),d1
-12: lsrl IMM (1),d0
- addl IMM (1),d4
-#endif
-1:
-| Now round, check for over- and underflow, and exit.
- movel a0,d7 | restore sign bit to d7
- movew IMM (DIVIDE),d5
- bra Lround$exit
-
-Ldivdf$inop:
- movew IMM (DIVIDE),d5
- bra Ld$inop
-
-Ldivdf$a$0:
-| If a is zero check to see whether b is zero also. In that case return
-| NaN; then check if b is NaN, and return NaN also in that case. Else
-| return zero.
- movew IMM (DIVIDE),d5
- bclr IMM (31),d2 |
- movel d2,d4 |
- orl d3,d4 |
- beq Ld$inop | if b is also zero return NaN
- cmpl IMM (0x7ff00000),d2 | check for NaN
- bhi Ld$inop |
- blt 1f |
- tstl d3 |
- bne Ld$inop |
-1: movel IMM (0),d0 | else return zero
- movel d0,d1 |
- PICLEA SYM (_fpCCR),a0 | clear exception flags
- movew IMM (0),a0@ |
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 |
-#else
- moveml sp@,d2-d7 |
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 |
- rts |
-
-Ldivdf$b$0:
- movew IMM (DIVIDE),d5
-| If we got here a is not zero. Check if a is NaN; in that case return NaN,
-| else return +/-INFINITY. Remember that a is in d0 with the sign bit
-| cleared already.
- movel a0,d7 | put a's sign bit back in d7 '
- cmpl IMM (0x7ff00000),d0 | compare d0 with INFINITY
- bhi Ld$inop | if larger it is NaN
- tstl d1 |
- bne Ld$inop |
- bra Ld$div$0 | else signal DIVIDE_BY_ZERO
-
-Ldivdf$b$nf:
- movew IMM (DIVIDE),d5
-| If d2 == 0x7ff00000 we have to check d3.
- tstl d3 |
- bne Ld$inop | if d3 <> 0, b is NaN
- bra Ld$underflow | else b is +/-INFINITY, so signal underflow
-
-Ldivdf$a$nf:
- movew IMM (DIVIDE),d5
-| If d0 == 0x7ff00000 we have to check d1.
- tstl d1 |
- bne Ld$inop | if d1 <> 0, a is NaN
-| If a is INFINITY we have to check b
- cmpl d7,d2 | compare b with INFINITY
- bge Ld$inop | if b is NaN or INFINITY return NaN
- tstl d3 |
- bne Ld$inop |
- bra Ld$overflow | else return overflow
-
-| If a number is denormalized we put an exponent of 1 but do not put the
-| bit back into the fraction.
-Ldivdf$a$den:
- movel IMM (1),d4
- andl d6,d0
-1: addl d1,d1 | shift a left until bit 20 is set
- addxl d0,d0
-#ifndef __mcoldfire__
- subw IMM (1),d4 | and adjust exponent
-#else
- subl IMM (1),d4 | and adjust exponent
-#endif
- btst IMM (DBL_MANT_DIG-32-1),d0
- bne Ldivdf$1
- bra 1b
-
-Ldivdf$b$den:
- movel IMM (1),d5
- andl d6,d2
-1: addl d3,d3 | shift b left until bit 20 is set
- addxl d2,d2
-#ifndef __mcoldfire__
- subw IMM (1),d5 | and adjust exponent
-#else
- subql IMM (1),d5 | and adjust exponent
-#endif
- btst IMM (DBL_MANT_DIG-32-1),d2
- bne Ldivdf$2
- bra 1b
-
-Lround$exit:
-| This is a common exit point for __muldf3 and __divdf3. When they enter
-| this point the sign of the result is in d7, the result in d0-d1, normalized
-| so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.
-
-| First check for underlow in the exponent:
-#ifndef __mcoldfire__
- cmpw IMM (-DBL_MANT_DIG-1),d4
-#else
- cmpl IMM (-DBL_MANT_DIG-1),d4
-#endif
- blt Ld$underflow
-| It could happen that the exponent is less than 1, in which case the
-| number is denormalized. In this case we shift right and adjust the
-| exponent until it becomes 1 or the fraction is zero (in the latter case
-| we signal underflow and return zero).
- movel d7,a0 |
- movel IMM (0),d6 | use d6-d7 to collect bits flushed right
- movel d6,d7 | use d6-d7 to collect bits flushed right
-#ifndef __mcoldfire__
- cmpw IMM (1),d4 | if the exponent is less than 1 we
-#else
- cmpl IMM (1),d4 | if the exponent is less than 1 we
-#endif
- bge 2f | have to shift right (denormalize)
-1:
-#ifndef __mcoldfire__
- addw IMM (1),d4 | adjust the exponent
- lsrl IMM (1),d0 | shift right once
- roxrl IMM (1),d1 |
- roxrl IMM (1),d2 |
- roxrl IMM (1),d3 |
- roxrl IMM (1),d6 |
- roxrl IMM (1),d7 |
- cmpw IMM (1),d4 | is the exponent 1 already?
-#else
- addl IMM (1),d4 | adjust the exponent
- lsrl IMM (1),d7
- btst IMM (0),d6
- beq 13f
- bset IMM (31),d7
-13: lsrl IMM (1),d6
- btst IMM (0),d3
- beq 14f
- bset IMM (31),d6
-14: lsrl IMM (1),d3
- btst IMM (0),d2
- beq 10f
- bset IMM (31),d3
-10: lsrl IMM (1),d2
- btst IMM (0),d1
- beq 11f
- bset IMM (31),d2
-11: lsrl IMM (1),d1
- btst IMM (0),d0
- beq 12f
- bset IMM (31),d1
-12: lsrl IMM (1),d0
- cmpl IMM (1),d4 | is the exponent 1 already?
-#endif
- beq 2f | if not loop back
- bra 1b |
- bra Ld$underflow | safety check, shouldn't execute '
-2: orl d6,d2 | this is a trick so we don't lose '
- orl d7,d3 | the bits which were flushed right
- movel a0,d7 | get back sign bit into d7
-| Now call the rounding routine (which takes care of denormalized numbers):
- lea pc@(Lround$0),a0 | to return from rounding routine
- PICLEA SYM (_fpCCR),a1 | check the rounding mode
-#ifdef __mcoldfire__
- clrl d6
-#endif
- movew a1@(6),d6 | rounding mode in d6
- beq Lround$to$nearest
-#ifndef __mcoldfire__
- cmpw IMM (ROUND_TO_PLUS),d6
-#else
- cmpl IMM (ROUND_TO_PLUS),d6
-#endif
- bhi Lround$to$minus
- blt Lround$to$zero
- bra Lround$to$plus
-Lround$0:
-| Here we have a correctly rounded result (either normalized or denormalized).
-
-| Here we should have either a normalized number or a denormalized one, and
-| the exponent is necessarily larger or equal to 1 (so we don't have to '
-| check again for underflow!). We have to check for overflow or for a
-| denormalized number (which also signals underflow).
-| Check for overflow (i.e., exponent >= 0x7ff).
-#ifndef __mcoldfire__
- cmpw IMM (0x07ff),d4
-#else
- cmpl IMM (0x07ff),d4
-#endif
- bge Ld$overflow
-| Now check for a denormalized number (exponent==0):
- movew d4,d4
- beq Ld$den
-1:
-| Put back the exponents and sign and return.
-#ifndef __mcoldfire__
- lslw IMM (4),d4 | exponent back to fourth byte
-#else
- lsll IMM (4),d4 | exponent back to fourth byte
-#endif
- bclr IMM (DBL_MANT_DIG-32-1),d0
- swap d0 | and put back exponent
-#ifndef __mcoldfire__
- orw d4,d0 |
-#else
- orl d4,d0 |
-#endif
- swap d0 |
- orl d7,d0 | and sign also
-
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-
-|=============================================================================
-| __negdf2
-|=============================================================================
-
-| double __negdf2(double, double);
-SYM (__negdf2):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@-
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movew IMM (NEGATE),d5
- movel a6@(8),d0 | get number to negate in d0-d1
- movel a6@(12),d1 |
- bchg IMM (31),d0 | negate
- movel d0,d2 | make a positive copy (for the tests)
- bclr IMM (31),d2 |
- movel d2,d4 | check for zero
- orl d1,d4 |
- beq 2f | if zero (either sign) return +zero
- cmpl IMM (0x7ff00000),d2 | compare to +INFINITY
- blt 1f | if finite, return
- bhi Ld$inop | if larger (fraction not zero) is NaN
- tstl d1 | if d2 == 0x7ff00000 check d1
- bne Ld$inop |
- movel d0,d7 | else get sign and return INFINITY
- andl IMM (0x80000000),d7
- bra Ld$infty
-1: PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-2: bclr IMM (31),d0
- bra 1b
-
-|=============================================================================
-| __cmpdf2
-|=============================================================================
-
-GREATER = 1
-LESS = -1
-EQUAL = 0
-
-| int __cmpdf2(double, double);
-SYM (__cmpdf2):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@- | save registers
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movew IMM (COMPARE),d5
- movel a6@(8),d0 | get first operand
- movel a6@(12),d1 |
- movel a6@(16),d2 | get second operand
- movel a6@(20),d3 |
-| First check if a and/or b are (+/-) zero and in that case clear
-| the sign bit.
- movel d0,d6 | copy signs into d6 (a) and d7(b)
- bclr IMM (31),d0 | and clear signs in d0 and d2
- movel d2,d7 |
- bclr IMM (31),d2 |
- cmpl IMM (0x7fff0000),d0 | check for a == NaN
- bhi Ld$inop | if d0 > 0x7ff00000, a is NaN
- beq Lcmpdf$a$nf | if equal can be INFINITY, so check d1
- movel d0,d4 | copy into d4 to test for zero
- orl d1,d4 |
- beq Lcmpdf$a$0 |
-Lcmpdf$0:
- cmpl IMM (0x7fff0000),d2 | check for b == NaN
- bhi Ld$inop | if d2 > 0x7ff00000, b is NaN
- beq Lcmpdf$b$nf | if equal can be INFINITY, so check d3
- movel d2,d4 |
- orl d3,d4 |
- beq Lcmpdf$b$0 |
-Lcmpdf$1:
-| Check the signs
- eorl d6,d7
- bpl 1f
-| If the signs are not equal check if a >= 0
- tstl d6
- bpl Lcmpdf$a$gt$b | if (a >= 0 && b < 0) => a > b
- bmi Lcmpdf$b$gt$a | if (a < 0 && b >= 0) => a < b
-1:
-| If the signs are equal check for < 0
- tstl d6
- bpl 1f
-| If both are negative exchange them
-#ifndef __mcoldfire__
- exg d0,d2
- exg d1,d3
-#else
- movel d0,d7
- movel d2,d0
- movel d7,d2
- movel d1,d7
- movel d3,d1
- movel d7,d3
-#endif
-1:
-| Now that they are positive we just compare them as longs (does this also
-| work for denormalized numbers?).
- cmpl d0,d2
- bhi Lcmpdf$b$gt$a | |b| > |a|
- bne Lcmpdf$a$gt$b | |b| < |a|
-| If we got here d0 == d2, so we compare d1 and d3.
- cmpl d1,d3
- bhi Lcmpdf$b$gt$a | |b| > |a|
- bne Lcmpdf$a$gt$b | |b| < |a|
-| If we got here a == b.
- movel IMM (EQUAL),d0
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | put back the registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-Lcmpdf$a$gt$b:
- movel IMM (GREATER),d0
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | put back the registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-Lcmpdf$b$gt$a:
- movel IMM (LESS),d0
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | put back the registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-
-Lcmpdf$a$0:
- bclr IMM (31),d6
- bra Lcmpdf$0
-Lcmpdf$b$0:
- bclr IMM (31),d7
- bra Lcmpdf$1
-
-Lcmpdf$a$nf:
- tstl d1
- bne Ld$inop
- bra Lcmpdf$0
-
-Lcmpdf$b$nf:
- tstl d3
- bne Ld$inop
- bra Lcmpdf$1
-
-|=============================================================================
-| rounding routines
-|=============================================================================
-
-| The rounding routines expect the number to be normalized in registers
-| d0-d1-d2-d3, with the exponent in register d4. They assume that the
-| exponent is larger or equal to 1. They return a properly normalized number
-| if possible, and a denormalized number otherwise. The exponent is returned
-| in d4.
-
-Lround$to$nearest:
-| We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
-| Here we assume that the exponent is not too small (this should be checked
-| before entering the rounding routine), but the number could be denormalized.
-
-| Check for denormalized numbers:
-1: btst IMM (DBL_MANT_DIG-32),d0
- bne 2f | if set the number is normalized
-| Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent
-| is one (remember that a denormalized number corresponds to an
-| exponent of -D_BIAS+1).
-#ifndef __mcoldfire__
- cmpw IMM (1),d4 | remember that the exponent is at least one
-#else
- cmpl IMM (1),d4 | remember that the exponent is at least one
-#endif
- beq 2f | an exponent of one means denormalized
- addl d3,d3 | else shift and adjust the exponent
- addxl d2,d2 |
- addxl d1,d1 |
- addxl d0,d0 |
-#ifndef __mcoldfire__
- dbra d4,1b |
-#else
- subql IMM (1), d4
- bpl 1b
-#endif
-2:
-| Now round: we do it as follows: after the shifting we can write the
-| fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
-| If delta < 1, do nothing. If delta > 1, add 1 to f.
-| If delta == 1, we make sure the rounded number will be even (odd?)
-| (after shifting).
- btst IMM (0),d1 | is delta < 1?
- beq 2f | if so, do not do anything
- orl d2,d3 | is delta == 1?
- bne 1f | if so round to even
- movel d1,d3 |
- andl IMM (2),d3 | bit 1 is the last significant bit
- movel IMM (0),d2 |
- addl d3,d1 |
- addxl d2,d0 |
- bra 2f |
-1: movel IMM (1),d3 | else add 1
- movel IMM (0),d2 |
- addl d3,d1 |
- addxl d2,d0
-| Shift right once (because we used bit #DBL_MANT_DIG-32!).
-2:
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
-#else
- lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
-#endif
-
-| Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
-| 'fraction overflow' ...).
- btst IMM (DBL_MANT_DIG-32),d0
- beq 1f
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- addw IMM (1),d4
-#else
- lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
- addl IMM (1),d4
-#endif
-1:
-| If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we
-| have to put the exponent to zero and return a denormalized number.
- btst IMM (DBL_MANT_DIG-32-1),d0
- beq 1f
- jmp a0@
-1: movel IMM (0),d4
- jmp a0@
-
-Lround$to$zero:
-Lround$to$plus:
-Lround$to$minus:
- jmp a0@
-#endif /* L_double */
-
-#ifdef L_float
-
- .globl SYM (_fpCCR)
- .globl $_exception_handler
-
-QUIET_NaN = 0xffffffff
-SIGNL_NaN = 0x7f800001
-INFINITY = 0x7f800000
-
-F_MAX_EXP = 0xff
-F_BIAS = 126
-FLT_MAX_EXP = F_MAX_EXP - F_BIAS
-FLT_MIN_EXP = 1 - F_BIAS
-FLT_MANT_DIG = 24
-
-INEXACT_RESULT = 0x0001
-UNDERFLOW = 0x0002
-OVERFLOW = 0x0004
-DIVIDE_BY_ZERO = 0x0008
-INVALID_OPERATION = 0x0010
-
-SINGLE_FLOAT = 1
-
-NOOP = 0
-ADD = 1
-MULTIPLY = 2
-DIVIDE = 3
-NEGATE = 4
-COMPARE = 5
-EXTENDSFDF = 6
-TRUNCDFSF = 7
-
-UNKNOWN = -1
-ROUND_TO_NEAREST = 0 | round result to nearest representable value
-ROUND_TO_ZERO = 1 | round result towards zero
-ROUND_TO_PLUS = 2 | round result towards plus infinity
-ROUND_TO_MINUS = 3 | round result towards minus infinity
-
-| Entry points:
-
- .globl SYM (__addsf3)
- .globl SYM (__subsf3)
- .globl SYM (__mulsf3)
- .globl SYM (__divsf3)
- .globl SYM (__negsf2)
- .globl SYM (__cmpsf2)
-
-| These are common routines to return and signal exceptions.
-
- .text
- .even
-
-Lf$den:
-| Return and signal a denormalized number
- orl d7,d0
- movew IMM (INEXACT_RESULT+UNDERFLOW),d7
- moveq IMM (SINGLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Lf$infty:
-Lf$overflow:
-| Return a properly signed INFINITY and set the exception flags
- movel IMM (INFINITY),d0
- orl d7,d0
- movew IMM (INEXACT_RESULT+OVERFLOW),d7
- moveq IMM (SINGLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Lf$underflow:
-| Return 0 and set the exception flags
- movel IMM (0),d0
- movew IMM (INEXACT_RESULT+UNDERFLOW),d7
- moveq IMM (SINGLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Lf$inop:
-| Return a quiet NaN and set the exception flags
- movel IMM (QUIET_NaN),d0
- movew IMM (INEXACT_RESULT+INVALID_OPERATION),d7
- moveq IMM (SINGLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-Lf$div$0:
-| Return a properly signed INFINITY and set the exception flags
- movel IMM (INFINITY),d0
- orl d7,d0
- movew IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
- moveq IMM (SINGLE_FLOAT),d6
- PICJUMP $_exception_handler
-
-|=============================================================================
-|=============================================================================
-| single precision routines
-|=============================================================================
-|=============================================================================
-
-| A single precision floating point number (float) has the format:
-|
-| struct _float {
-| unsigned int sign : 1; /* sign bit */
-| unsigned int exponent : 8; /* exponent, shifted by 126 */
-| unsigned int fraction : 23; /* fraction */
-| } float;
-|
-| Thus sizeof(float) = 4 (32 bits).
-|
-| All the routines are callable from C programs, and return the result
-| in the single register d0. They also preserve all registers except
-| d0-d1 and a0-a1.
-
-|=============================================================================
-| __subsf3
-|=============================================================================
-
-| float __subsf3(float, float);
-SYM (__subsf3):
- bchg IMM (31),sp@(8) | change sign of second operand
- | and fall through
-|=============================================================================
-| __addsf3
-|=============================================================================
-
-| float __addsf3(float, float);
-SYM (__addsf3):
-#ifndef __mcoldfire__
- link a6,IMM (0) | everything will be done in registers
- moveml d2-d7,sp@- | save all data registers but d0-d1
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movel a6@(8),d0 | get first operand
- movel a6@(12),d1 | get second operand
- movel d0,d6 | get d0's sign bit '
- addl d0,d0 | check and clear sign bit of a
- beq Laddsf$b | if zero return second operand
- movel d1,d7 | save b's sign bit '
- addl d1,d1 | get rid of sign bit
- beq Laddsf$a | if zero return first operand
-
- movel d6,a0 | save signs in address registers
- movel d7,a1 | so we can use d6 and d7
-
-| Get the exponents and check for denormalized and/or infinity.
-
- movel IMM (0x00ffffff),d4 | mask to get fraction
- movel IMM (0x01000000),d5 | mask to put hidden bit back
-
- movel d0,d6 | save a to get exponent
- andl d4,d0 | get fraction in d0
- notl d4 | make d4 into a mask for the exponent
- andl d4,d6 | get exponent in d6
- beq Laddsf$a$den | branch if a is denormalized
- cmpl d4,d6 | check for INFINITY or NaN
- beq Laddsf$nf
- swap d6 | put exponent into first word
- orl d5,d0 | and put hidden bit back
-Laddsf$1:
-| Now we have a's exponent in d6 (second byte) and the mantissa in d0. '
- movel d1,d7 | get exponent in d7
- andl d4,d7 |
- beq Laddsf$b$den | branch if b is denormalized
- cmpl d4,d7 | check for INFINITY or NaN
- beq Laddsf$nf
- swap d7 | put exponent into first word
- notl d4 | make d4 into a mask for the fraction
- andl d4,d1 | get fraction in d1
- orl d5,d1 | and put hidden bit back
-Laddsf$2:
-| Now we have b's exponent in d7 (second byte) and the mantissa in d1. '
-
-| Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we
-| shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
-| bit).
-
- movel d1,d2 | move b to d2, since we want to use
- | two registers to do the sum
- movel IMM (0),d1 | and clear the new ones
- movel d1,d3 |
-
-| Here we shift the numbers in registers d0 and d1 so the exponents are the
-| same, and put the largest exponent in d6. Note that we are using two
-| registers for each number (see the discussion by D. Knuth in "Seminumerical
-| Algorithms").
-#ifndef __mcoldfire__
- cmpw d6,d7 | compare exponents
-#else
- cmpl d6,d7 | compare exponents
-#endif
- beq Laddsf$3 | if equal don't shift '
- bhi 5f | branch if second exponent largest
-1:
- subl d6,d7 | keep the largest exponent
- negl d7
-#ifndef __mcoldfire__
- lsrw IMM (8),d7 | put difference in lower byte
-#else
- lsrl IMM (8),d7 | put difference in lower byte
-#endif
-| if difference is too large we don't shift (actually, we can just exit) '
-#ifndef __mcoldfire__
- cmpw IMM (FLT_MANT_DIG+2),d7
-#else
- cmpl IMM (FLT_MANT_DIG+2),d7
-#endif
- bge Laddsf$b$small
-#ifndef __mcoldfire__
- cmpw IMM (16),d7 | if difference >= 16 swap
-#else
- cmpl IMM (16),d7 | if difference >= 16 swap
-#endif
- bge 4f
-2:
-#ifndef __mcoldfire__
- subw IMM (1),d7
-#else
- subql IMM (1), d7
-#endif
-3:
-#ifndef __mcoldfire__
- lsrl IMM (1),d2 | shift right second operand
- roxrl IMM (1),d3
- dbra d7,3b
-#else
- lsrl IMM (1),d3
- btst IMM (0),d2
- beq 10f
- bset IMM (31),d3
-10: lsrl IMM (1),d2
- subql IMM (1), d7
- bpl 3b
-#endif
- bra Laddsf$3
-4:
- movew d2,d3
- swap d3
- movew d3,d2
- swap d2
-#ifndef __mcoldfire__
- subw IMM (16),d7
-#else
- subl IMM (16),d7
-#endif
- bne 2b | if still more bits, go back to normal case
- bra Laddsf$3
-5:
-#ifndef __mcoldfire__
- exg d6,d7 | exchange the exponents
-#else
- eorl d6,d7
- eorl d7,d6
- eorl d6,d7
-#endif
- subl d6,d7 | keep the largest exponent
- negl d7 |
-#ifndef __mcoldfire__
- lsrw IMM (8),d7 | put difference in lower byte
-#else
- lsrl IMM (8),d7 | put difference in lower byte
-#endif
-| if difference is too large we don't shift (and exit!) '
-#ifndef __mcoldfire__
- cmpw IMM (FLT_MANT_DIG+2),d7
-#else
- cmpl IMM (FLT_MANT_DIG+2),d7
-#endif
- bge Laddsf$a$small
-#ifndef __mcoldfire__
- cmpw IMM (16),d7 | if difference >= 16 swap
-#else
- cmpl IMM (16),d7 | if difference >= 16 swap
-#endif
- bge 8f
-6:
-#ifndef __mcoldfire__
- subw IMM (1),d7
-#else
- subl IMM (1),d7
-#endif
-7:
-#ifndef __mcoldfire__
- lsrl IMM (1),d0 | shift right first operand
- roxrl IMM (1),d1
- dbra d7,7b
-#else
- lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
- subql IMM (1),d7
- bpl 7b
-#endif
- bra Laddsf$3
-8:
- movew d0,d1
- swap d1
- movew d1,d0
- swap d0
-#ifndef __mcoldfire__
- subw IMM (16),d7
-#else
- subl IMM (16),d7
-#endif
- bne 6b | if still more bits, go back to normal case
- | otherwise we fall through
-
-| Now we have a in d0-d1, b in d2-d3, and the largest exponent in d6 (the
-| signs are stored in a0 and a1).
-
-Laddsf$3:
-| Here we have to decide whether to add or subtract the numbers
-#ifndef __mcoldfire__
- exg d6,a0 | get signs back
- exg d7,a1 | and save the exponents
-#else
- movel d6,d4
- movel a0,d6
- movel d4,a0
- movel d7,d4
- movel a1,d7
- movel d4,a1
-#endif
- eorl d6,d7 | combine sign bits
- bmi Lsubsf$0 | if negative a and b have opposite
- | sign so we actually subtract the
- | numbers
-
-| Here we have both positive or both negative
-#ifndef __mcoldfire__
- exg d6,a0 | now we have the exponent in d6
-#else
- movel d6,d4
- movel a0,d6
- movel d4,a0
-#endif
- movel a0,d7 | and sign in d7
- andl IMM (0x80000000),d7
-| Here we do the addition.
- addl d3,d1
- addxl d2,d0
-| Note: now we have d2, d3, d4 and d5 to play with!
-
-| Put the exponent, in the first byte, in d2, to use the "standard" rounding
-| routines:
- movel d6,d2
-#ifndef __mcoldfire__
- lsrw IMM (8),d2
-#else
- lsrl IMM (8),d2
-#endif
-
-| Before rounding normalize so bit #FLT_MANT_DIG is set (we will consider
-| the case of denormalized numbers in the rounding routine itself).
-| As in the addition (not in the subtraction!) we could have set
-| one more bit we check this:
- btst IMM (FLT_MANT_DIG+1),d0
- beq 1f
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
-#else
- lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
-#endif
- addl IMM (1),d2
-1:
- lea pc@(Laddsf$4),a0 | to return from rounding routine
- PICLEA SYM (_fpCCR),a1 | check the rounding mode
-#ifdef __mcoldfire__
- clrl d6
-#endif
- movew a1@(6),d6 | rounding mode in d6
- beq Lround$to$nearest
-#ifndef __mcoldfire__
- cmpw IMM (ROUND_TO_PLUS),d6
-#else
- cmpl IMM (ROUND_TO_PLUS),d6
-#endif
- bhi Lround$to$minus
- blt Lround$to$zero
- bra Lround$to$plus
-Laddsf$4:
-| Put back the exponent, but check for overflow.
-#ifndef __mcoldfire__
- cmpw IMM (0xff),d2
-#else
- cmpl IMM (0xff),d2
-#endif
- bhi 1f
- bclr IMM (FLT_MANT_DIG-1),d0
-#ifndef __mcoldfire__
- lslw IMM (7),d2
-#else
- lsll IMM (7),d2
-#endif
- swap d2
- orl d2,d0
- bra Laddsf$ret
-1:
- movew IMM (ADD),d5
- bra Lf$overflow
-
-Lsubsf$0:
-| We are here if a > 0 and b < 0 (sign bits cleared).
-| Here we do the subtraction.
- movel d6,d7 | put sign in d7
- andl IMM (0x80000000),d7
-
- subl d3,d1 | result in d0-d1
- subxl d2,d0 |
- beq Laddsf$ret | if zero just exit
- bpl 1f | if positive skip the following
- bchg IMM (31),d7 | change sign bit in d7
- negl d1
- negxl d0
-1:
-#ifndef __mcoldfire__
- exg d2,a0 | now we have the exponent in d2
- lsrw IMM (8),d2 | put it in the first byte
-#else
- movel d2,d4
- movel a0,d2
- movel d4,a0
- lsrl IMM (8),d2 | put it in the first byte
-#endif
-
-| Now d0-d1 is positive and the sign bit is in d7.
-
-| Note that we do not have to normalize, since in the subtraction bit
-| #FLT_MANT_DIG+1 is never set, and denormalized numbers are handled by
-| the rounding routines themselves.
- lea pc@(Lsubsf$1),a0 | to return from rounding routine
- PICLEA SYM (_fpCCR),a1 | check the rounding mode
-#ifdef __mcoldfire__
- clrl d6
-#endif
- movew a1@(6),d6 | rounding mode in d6
- beq Lround$to$nearest
-#ifndef __mcoldfire__
- cmpw IMM (ROUND_TO_PLUS),d6
-#else
- cmpl IMM (ROUND_TO_PLUS),d6
-#endif
- bhi Lround$to$minus
- blt Lround$to$zero
- bra Lround$to$plus
-Lsubsf$1:
-| Put back the exponent (we can't have overflow!). '
- bclr IMM (FLT_MANT_DIG-1),d0
-#ifndef __mcoldfire__
- lslw IMM (7),d2
-#else
- lsll IMM (7),d2
-#endif
- swap d2
- orl d2,d0
- bra Laddsf$ret
-
-| If one of the numbers was too small (difference of exponents >=
-| FLT_MANT_DIG+2) we return the other (and now we don't have to '
-| check for finiteness or zero).
-Laddsf$a$small:
- movel a6@(12),d0
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | restore data registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 | and return
- rts
-
-Laddsf$b$small:
- movel a6@(8),d0
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | restore data registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 | and return
- rts
-
-| If the numbers are denormalized remember to put exponent equal to 1.
-
-Laddsf$a$den:
- movel d5,d6 | d5 contains 0x01000000
- swap d6
- bra Laddsf$1
-
-Laddsf$b$den:
- movel d5,d7
- swap d7
- notl d4 | make d4 into a mask for the fraction
- | (this was not executed after the jump)
- bra Laddsf$2
-
-| The rest is mainly code for the different results which can be
-| returned (checking always for +/-INFINITY and NaN).
-
-Laddsf$b:
-| Return b (if a is zero).
- movel a6@(12),d0
- bra 1f
-Laddsf$a:
-| Return a (if b is zero).
- movel a6@(8),d0
-1:
- movew IMM (ADD),d5
-| We have to check for NaN and +/-infty.
- movel d0,d7
- andl IMM (0x80000000),d7 | put sign in d7
- bclr IMM (31),d0 | clear sign
- cmpl IMM (INFINITY),d0 | check for infty or NaN
- bge 2f
- movel d0,d0 | check for zero (we do this because we don't '
- bne Laddsf$ret | want to return -0 by mistake
- bclr IMM (31),d7 | if zero be sure to clear sign
- bra Laddsf$ret | if everything OK just return
-2:
-| The value to be returned is either +/-infty or NaN
- andl IMM (0x007fffff),d0 | check for NaN
- bne Lf$inop | if mantissa not zero is NaN
- bra Lf$infty
-
-Laddsf$ret:
-| Normal exit (a and b nonzero, result is not NaN nor +/-infty).
-| We have to clear the exception flags (just the exception type).
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
- orl d7,d0 | put sign bit
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | restore data registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 | and return
- rts
-
-Laddsf$ret$den:
-| Return a denormalized number (for addition we don't signal underflow) '
- lsrl IMM (1),d0 | remember to shift right back once
- bra Laddsf$ret | and return
-
-| Note: when adding two floats of the same sign if either one is
-| NaN we return NaN without regard to whether the other is finite or
-| not. When subtracting them (i.e., when adding two numbers of
-| opposite signs) things are more complicated: if both are INFINITY
-| we return NaN, if only one is INFINITY and the other is NaN we return
-| NaN, but if it is finite we return INFINITY with the corresponding sign.
-
-Laddsf$nf:
- movew IMM (ADD),d5
-| This could be faster but it is not worth the effort, since it is not
-| executed very often. We sacrifice speed for clarity here.
- movel a6@(8),d0 | get the numbers back (remember that we
- movel a6@(12),d1 | did some processing already)
- movel IMM (INFINITY),d4 | useful constant (INFINITY)
- movel d0,d2 | save sign bits
- movel d1,d3
- bclr IMM (31),d0 | clear sign bits
- bclr IMM (31),d1
-| We know that one of them is either NaN of +/-INFINITY
-| Check for NaN (if either one is NaN return NaN)
- cmpl d4,d0 | check first a (d0)
- bhi Lf$inop
- cmpl d4,d1 | check now b (d1)
- bhi Lf$inop
-| Now comes the check for +/-INFINITY. We know that both are (maybe not
-| finite) numbers, but we have to check if both are infinite whether we
-| are adding or subtracting them.
- eorl d3,d2 | to check sign bits
- bmi 1f
- movel d0,d7
- andl IMM (0x80000000),d7 | get (common) sign bit
- bra Lf$infty
-1:
-| We know one (or both) are infinite, so we test for equality between the
-| two numbers (if they are equal they have to be infinite both, so we
-| return NaN).
- cmpl d1,d0 | are both infinite?
- beq Lf$inop | if so return NaN
-
- movel d0,d7
- andl IMM (0x80000000),d7 | get a's sign bit '
- cmpl d4,d0 | test now for infinity
- beq Lf$infty | if a is INFINITY return with this sign
- bchg IMM (31),d7 | else we know b is INFINITY and has
- bra Lf$infty | the opposite sign
-
-|=============================================================================
-| __mulsf3
-|=============================================================================
-
-| float __mulsf3(float, float);
-SYM (__mulsf3):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@-
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movel a6@(8),d0 | get a into d0
- movel a6@(12),d1 | and b into d1
- movel d0,d7 | d7 will hold the sign of the product
- eorl d1,d7 |
- andl IMM (0x80000000),d7
- movel IMM (INFINITY),d6 | useful constant (+INFINITY)
- movel d6,d5 | another (mask for fraction)
- notl d5 |
- movel IMM (0x00800000),d4 | this is to put hidden bit back
- bclr IMM (31),d0 | get rid of a's sign bit '
- movel d0,d2 |
- beq Lmulsf$a$0 | branch if a is zero
- bclr IMM (31),d1 | get rid of b's sign bit '
- movel d1,d3 |
- beq Lmulsf$b$0 | branch if b is zero
- cmpl d6,d0 | is a big?
- bhi Lmulsf$inop | if a is NaN return NaN
- beq Lmulsf$inf | if a is INFINITY we have to check b
- cmpl d6,d1 | now compare b with INFINITY
- bhi Lmulsf$inop | is b NaN?
- beq Lmulsf$overflow | is b INFINITY?
-| Here we have both numbers finite and nonzero (and with no sign bit).
-| Now we get the exponents into d2 and d3.
- andl d6,d2 | and isolate exponent in d2
- beq Lmulsf$a$den | if exponent is zero we have a denormalized
- andl d5,d0 | and isolate fraction
- orl d4,d0 | and put hidden bit back
- swap d2 | I like exponents in the first byte
-#ifndef __mcoldfire__
- lsrw IMM (7),d2 |
-#else
- lsrl IMM (7),d2 |
-#endif
-Lmulsf$1: | number
- andl d6,d3 |
- beq Lmulsf$b$den |
- andl d5,d1 |
- orl d4,d1 |
- swap d3 |
-#ifndef __mcoldfire__
- lsrw IMM (7),d3 |
-#else
- lsrl IMM (7),d3 |
-#endif
-Lmulsf$2: |
-#ifndef __mcoldfire__
- addw d3,d2 | add exponents
- subw IMM (F_BIAS+1),d2 | and subtract bias (plus one)
-#else
- addl d3,d2 | add exponents
- subl IMM (F_BIAS+1),d2 | and subtract bias (plus one)
-#endif
-
-| We are now ready to do the multiplication. The situation is as follows:
-| both a and b have bit FLT_MANT_DIG-1 set (even if they were
-| denormalized to start with!), which means that in the product
-| bit 2*(FLT_MANT_DIG-1) (that is, bit 2*FLT_MANT_DIG-2-32 of the
-| high long) is set.
-
-| To do the multiplication let us move the number a little bit around ...
- movel d1,d6 | second operand in d6
- movel d0,d5 | first operand in d4-d5
- movel IMM (0),d4
- movel d4,d1 | the sums will go in d0-d1
- movel d4,d0
-
-| now bit FLT_MANT_DIG-1 becomes bit 31:
- lsll IMM (31-FLT_MANT_DIG+1),d6
-
-| Start the loop (we loop #FLT_MANT_DIG times):
- movew IMM (FLT_MANT_DIG-1),d3
-1: addl d1,d1 | shift sum
- addxl d0,d0
- lsll IMM (1),d6 | get bit bn
- bcc 2f | if not set skip sum
- addl d5,d1 | add a
- addxl d4,d0
-2:
-#ifndef __mcoldfire__
- dbf d3,1b | loop back
-#else
- subql IMM (1),d3
- bpl 1b
-#endif
-
-| Now we have the product in d0-d1, with bit (FLT_MANT_DIG - 1) + FLT_MANT_DIG
-| (mod 32) of d0 set. The first thing to do now is to normalize it so bit
-| FLT_MANT_DIG is set (to do the rounding).
-#ifndef __mcoldfire__
- rorl IMM (6),d1
- swap d1
- movew d1,d3
- andw IMM (0x03ff),d3
- andw IMM (0xfd00),d1
-#else
- movel d1,d3
- lsll IMM (8),d1
- addl d1,d1
- addl d1,d1
- moveq IMM (22),d5
- lsrl d5,d3
- orl d3,d1
- andl IMM (0xfffffd00),d1
-#endif
- lsll IMM (8),d0
- addl d0,d0
- addl d0,d0
-#ifndef __mcoldfire__
- orw d3,d0
-#else
- orl d3,d0
-#endif
-
- movew IMM (MULTIPLY),d5
-
- btst IMM (FLT_MANT_DIG+1),d0
- beq Lround$exit
-#ifndef __mcoldfire__
- lsrl IMM (1),d0
- roxrl IMM (1),d1
- addw IMM (1),d2
-#else
- lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
- addql IMM (1),d2
-#endif
- bra Lround$exit
-
-Lmulsf$inop:
- movew IMM (MULTIPLY),d5
- bra Lf$inop
-
-Lmulsf$overflow:
- movew IMM (MULTIPLY),d5
- bra Lf$overflow
-
-Lmulsf$inf:
- movew IMM (MULTIPLY),d5
-| If either is NaN return NaN; else both are (maybe infinite) numbers, so
-| return INFINITY with the correct sign (which is in d7).
- cmpl d6,d1 | is b NaN?
- bhi Lf$inop | if so return NaN
- bra Lf$overflow | else return +/-INFINITY
-
-| If either number is zero return zero, unless the other is +/-INFINITY,
-| or NaN, in which case we return NaN.
-Lmulsf$b$0:
-| Here d1 (==b) is zero.
- movel d1,d0 | put b into d0 (just a zero)
- movel a6@(8),d1 | get a again to check for non-finiteness
- bra 1f
-Lmulsf$a$0:
- movel a6@(12),d1 | get b again to check for non-finiteness
-1: bclr IMM (31),d1 | clear sign bit
- cmpl IMM (INFINITY),d1 | and check for a large exponent
- bge Lf$inop | if b is +/-INFINITY or NaN return NaN
- PICLEA SYM (_fpCCR),a0 | else return zero
- movew IMM (0),a0@ |
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 |
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 |
- rts |
-
-| If a number is denormalized we put an exponent of 1 but do not put the
-| hidden bit back into the fraction; instead we shift left until bit 23
-| (the hidden bit) is set, adjusting the exponent accordingly. We do this
-| to ensure that the product of the fractions is close to 1.
-Lmulsf$a$den:
- movel IMM (1),d2
- andl d5,d0
-1: addl d0,d0 | shift a left (until bit 23 is set)
-#ifndef __mcoldfire__
- subw IMM (1),d2 | and adjust exponent
-#else
- subql IMM (1),d2 | and adjust exponent
-#endif
- btst IMM (FLT_MANT_DIG-1),d0
- bne Lmulsf$1 |
- bra 1b | else loop back
-
-Lmulsf$b$den:
- movel IMM (1),d3
- andl d5,d1
-1: addl d1,d1 | shift b left until bit 23 is set
-#ifndef __mcoldfire__
- subw IMM (1),d3 | and adjust exponent
-#else
- subl IMM (1),d3 | and adjust exponent
-#endif
- btst IMM (FLT_MANT_DIG-1),d1
- bne Lmulsf$2 |
- bra 1b | else loop back
-
-|=============================================================================
-| __divsf3
-|=============================================================================
-
-| float __divsf3(float, float);
-SYM (__divsf3):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@-
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movel a6@(8),d0 | get a into d0
- movel a6@(12),d1 | and b into d1
- movel d0,d7 | d7 will hold the sign of the result
- eorl d1,d7 |
- andl IMM (0x80000000),d7 |
- movel IMM (INFINITY),d6 | useful constant (+INFINITY)
- movel d6,d5 | another (mask for fraction)
- notl d5 |
- movel IMM (0x00800000),d4 | this is to put hidden bit back
- bclr IMM (31),d0 | get rid of a's sign bit '
- movel d0,d2 |
- beq Ldivsf$a$0 | branch if a is zero
- bclr IMM (31),d1 | get rid of b's sign bit '
- movel d1,d3 |
- beq Ldivsf$b$0 | branch if b is zero
- cmpl d6,d0 | is a big?
- bhi Ldivsf$inop | if a is NaN return NaN
- beq Ldivsf$inf | if a is INFINITY we have to check b
- cmpl d6,d1 | now compare b with INFINITY
- bhi Ldivsf$inop | if b is NaN return NaN
- beq Ldivsf$underflow
-| Here we have both numbers finite and nonzero (and with no sign bit).
-| Now we get the exponents into d2 and d3 and normalize the numbers to
-| ensure that the ratio of the fractions is close to 1. We do this by
-| making sure that bit #FLT_MANT_DIG-1 (hidden bit) is set.
- andl d6,d2 | and isolate exponent in d2
- beq Ldivsf$a$den | if exponent is zero we have a denormalized
- andl d5,d0 | and isolate fraction
- orl d4,d0 | and put hidden bit back
- swap d2 | I like exponents in the first byte
-#ifndef __mcoldfire__
- lsrw IMM (7),d2 |
-#else
- lsrl IMM (7),d2 |
-#endif
-Ldivsf$1: |
- andl d6,d3 |
- beq Ldivsf$b$den |
- andl d5,d1 |
- orl d4,d1 |
- swap d3 |
-#ifndef __mcoldfire__
- lsrw IMM (7),d3 |
-#else
- lsrl IMM (7),d3 |
-#endif
-Ldivsf$2: |
-#ifndef __mcoldfire__
- subw d3,d2 | subtract exponents
- addw IMM (F_BIAS),d2 | and add bias
-#else
- subl d3,d2 | subtract exponents
- addl IMM (F_BIAS),d2 | and add bias
-#endif
-
-| We are now ready to do the division. We have prepared things in such a way
-| that the ratio of the fractions will be less than 2 but greater than 1/2.
-| At this point the registers in use are:
-| d0 holds a (first operand, bit FLT_MANT_DIG=0, bit FLT_MANT_DIG-1=1)
-| d1 holds b (second operand, bit FLT_MANT_DIG=1)
-| d2 holds the difference of the exponents, corrected by the bias
-| d7 holds the sign of the ratio
-| d4, d5, d6 hold some constants
- movel d7,a0 | d6-d7 will hold the ratio of the fractions
- movel IMM (0),d6 |
- movel d6,d7
-
- movew IMM (FLT_MANT_DIG+1),d3
-1: cmpl d0,d1 | is a < b?
- bhi 2f |
- bset d3,d6 | set a bit in d6
- subl d1,d0 | if a >= b a <-- a-b
- beq 3f | if a is zero, exit
-2: addl d0,d0 | multiply a by 2
-#ifndef __mcoldfire__
- dbra d3,1b
-#else
- subql IMM (1),d3
- bpl 1b
-#endif
-
-| Now we keep going to set the sticky bit ...
- movew IMM (FLT_MANT_DIG),d3
-1: cmpl d0,d1
- ble 2f
- addl d0,d0
-#ifndef __mcoldfire__
- dbra d3,1b
-#else
- subql IMM(1),d3
- bpl 1b
-#endif
- movel IMM (0),d1
- bra 3f
-2: movel IMM (0),d1
-#ifndef __mcoldfire__
- subw IMM (FLT_MANT_DIG),d3
- addw IMM (31),d3
-#else
- subl IMM (FLT_MANT_DIG),d3
- addl IMM (31),d3
-#endif
- bset d3,d1
-3:
- movel d6,d0 | put the ratio in d0-d1
- movel a0,d7 | get sign back
-
-| Because of the normalization we did before we are guaranteed that
-| d0 is smaller than 2^26 but larger than 2^24. Thus bit 26 is not set,
-| bit 25 could be set, and if it is not set then bit 24 is necessarily set.
- btst IMM (FLT_MANT_DIG+1),d0
- beq 1f | if it is not set, then bit 24 is set
- lsrl IMM (1),d0 |
-#ifndef __mcoldfire__
- addw IMM (1),d2 |
-#else
- addl IMM (1),d2 |
-#endif
-1:
-| Now round, check for over- and underflow, and exit.
- movew IMM (DIVIDE),d5
- bra Lround$exit
-
-Ldivsf$inop:
- movew IMM (DIVIDE),d5
- bra Lf$inop
-
-Ldivsf$overflow:
- movew IMM (DIVIDE),d5
- bra Lf$overflow
-
-Ldivsf$underflow:
- movew IMM (DIVIDE),d5
- bra Lf$underflow
-
-Ldivsf$a$0:
- movew IMM (DIVIDE),d5
-| If a is zero check to see whether b is zero also. In that case return
-| NaN; then check if b is NaN, and return NaN also in that case. Else
-| return zero.
- andl IMM (0x7fffffff),d1 | clear sign bit and test b
- beq Lf$inop | if b is also zero return NaN
- cmpl IMM (INFINITY),d1 | check for NaN
- bhi Lf$inop |
- movel IMM (0),d0 | else return zero
- PICLEA SYM (_fpCCR),a0 |
- movew IMM (0),a0@ |
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 |
-#else
- moveml sp@,d2-d7 |
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6 |
- rts |
-
-Ldivsf$b$0:
- movew IMM (DIVIDE),d5
-| If we got here a is not zero. Check if a is NaN; in that case return NaN,
-| else return +/-INFINITY. Remember that a is in d0 with the sign bit
-| cleared already.
- cmpl IMM (INFINITY),d0 | compare d0 with INFINITY
- bhi Lf$inop | if larger it is NaN
- bra Lf$div$0 | else signal DIVIDE_BY_ZERO
-
-Ldivsf$inf:
- movew IMM (DIVIDE),d5
-| If a is INFINITY we have to check b
- cmpl IMM (INFINITY),d1 | compare b with INFINITY
- bge Lf$inop | if b is NaN or INFINITY return NaN
- bra Lf$overflow | else return overflow
-
-| If a number is denormalized we put an exponent of 1 but do not put the
-| bit back into the fraction.
-Ldivsf$a$den:
- movel IMM (1),d2
- andl d5,d0
-1: addl d0,d0 | shift a left until bit FLT_MANT_DIG-1 is set
-#ifndef __mcoldfire__
- subw IMM (1),d2 | and adjust exponent
-#else
- subl IMM (1),d2 | and adjust exponent
-#endif
- btst IMM (FLT_MANT_DIG-1),d0
- bne Ldivsf$1
- bra 1b
-
-Ldivsf$b$den:
- movel IMM (1),d3
- andl d5,d1
-1: addl d1,d1 | shift b left until bit FLT_MANT_DIG is set
-#ifndef __mcoldfire__
- subw IMM (1),d3 | and adjust exponent
-#else
- subl IMM (1),d3 | and adjust exponent
-#endif
- btst IMM (FLT_MANT_DIG-1),d1
- bne Ldivsf$2
- bra 1b
-
-Lround$exit:
-| This is a common exit point for __mulsf3 and __divsf3.
-
-| First check for underlow in the exponent:
-#ifndef __mcoldfire__
- cmpw IMM (-FLT_MANT_DIG-1),d2
-#else
- cmpl IMM (-FLT_MANT_DIG-1),d2
-#endif
- blt Lf$underflow
-| It could happen that the exponent is less than 1, in which case the
-| number is denormalized. In this case we shift right and adjust the
-| exponent until it becomes 1 or the fraction is zero (in the latter case
-| we signal underflow and return zero).
- movel IMM (0),d6 | d6 is used temporarily
-#ifndef __mcoldfire__
- cmpw IMM (1),d2 | if the exponent is less than 1 we
-#else
- cmpl IMM (1),d2 | if the exponent is less than 1 we
-#endif
- bge 2f | have to shift right (denormalize)
-1:
-#ifndef __mcoldfire__
- addw IMM (1),d2 | adjust the exponent
- lsrl IMM (1),d0 | shift right once
- roxrl IMM (1),d1 |
- roxrl IMM (1),d6 | d6 collect bits we would lose otherwise
- cmpw IMM (1),d2 | is the exponent 1 already?
-#else
- addql IMM (1),d2 | adjust the exponent
- lsrl IMM (1),d6
- btst IMM (0),d1
- beq 11f
- bset IMM (31),d6
-11: lsrl IMM (1),d1
- btst IMM (0),d0
- beq 10f
- bset IMM (31),d1
-10: lsrl IMM (1),d0
- cmpl IMM (1),d2 | is the exponent 1 already?
-#endif
- beq 2f | if not loop back
- bra 1b |
- bra Lf$underflow | safety check, shouldn't execute '
-2: orl d6,d1 | this is a trick so we don't lose '
- | the extra bits which were flushed right
-| Now call the rounding routine (which takes care of denormalized numbers):
- lea pc@(Lround$0),a0 | to return from rounding routine
- PICLEA SYM (_fpCCR),a1 | check the rounding mode
-#ifdef __mcoldfire__
- clrl d6
-#endif
- movew a1@(6),d6 | rounding mode in d6
- beq Lround$to$nearest
-#ifndef __mcoldfire__
- cmpw IMM (ROUND_TO_PLUS),d6
-#else
- cmpl IMM (ROUND_TO_PLUS),d6
-#endif
- bhi Lround$to$minus
- blt Lround$to$zero
- bra Lround$to$plus
-Lround$0:
-| Here we have a correctly rounded result (either normalized or denormalized).
-
-| Here we should have either a normalized number or a denormalized one, and
-| the exponent is necessarily larger or equal to 1 (so we don't have to '
-| check again for underflow!). We have to check for overflow or for a
-| denormalized number (which also signals underflow).
-| Check for overflow (i.e., exponent >= 255).
-#ifndef __mcoldfire__
- cmpw IMM (0x00ff),d2
-#else
- cmpl IMM (0x00ff),d2
-#endif
- bge Lf$overflow
-| Now check for a denormalized number (exponent==0).
- movew d2,d2
- beq Lf$den
-1:
-| Put back the exponents and sign and return.
-#ifndef __mcoldfire__
- lslw IMM (7),d2 | exponent back to fourth byte
-#else
- lsll IMM (7),d2 | exponent back to fourth byte
-#endif
- bclr IMM (FLT_MANT_DIG-1),d0
- swap d0 | and put back exponent
-#ifndef __mcoldfire__
- orw d2,d0 |
-#else
- orl d2,d0
-#endif
- swap d0 |
- orl d7,d0 | and sign also
-
- PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-
-|=============================================================================
-| __negsf2
-|=============================================================================
-
-| This is trivial and could be shorter if we didn't bother checking for NaN '
-| and +/-INFINITY.
-
-| float __negsf2(float);
-SYM (__negsf2):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@-
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movew IMM (NEGATE),d5
- movel a6@(8),d0 | get number to negate in d0
- bchg IMM (31),d0 | negate
- movel d0,d1 | make a positive copy
- bclr IMM (31),d1 |
- tstl d1 | check for zero
- beq 2f | if zero (either sign) return +zero
- cmpl IMM (INFINITY),d1 | compare to +INFINITY
- blt 1f |
- bhi Lf$inop | if larger (fraction not zero) is NaN
- movel d0,d7 | else get sign and return INFINITY
- andl IMM (0x80000000),d7
- bra Lf$infty
-1: PICLEA SYM (_fpCCR),a0
- movew IMM (0),a0@
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-2: bclr IMM (31),d0
- bra 1b
-
-|=============================================================================
-| __cmpsf2
-|=============================================================================
-
-GREATER = 1
-LESS = -1
-EQUAL = 0
-
-| int __cmpsf2(float, float);
-SYM (__cmpsf2):
-#ifndef __mcoldfire__
- link a6,IMM (0)
- moveml d2-d7,sp@- | save registers
-#else
- link a6,IMM (-24)
- moveml d2-d7,sp@
-#endif
- movew IMM (COMPARE),d5
- movel a6@(8),d0 | get first operand
- movel a6@(12),d1 | get second operand
-| Check if either is NaN, and in that case return garbage and signal
-| INVALID_OPERATION. Check also if either is zero, and clear the signs
-| if necessary.
- movel d0,d6
- andl IMM (0x7fffffff),d0
- beq Lcmpsf$a$0
- cmpl IMM (0x7f800000),d0
- bhi Lf$inop
-Lcmpsf$1:
- movel d1,d7
- andl IMM (0x7fffffff),d1
- beq Lcmpsf$b$0
- cmpl IMM (0x7f800000),d1
- bhi Lf$inop
-Lcmpsf$2:
-| Check the signs
- eorl d6,d7
- bpl 1f
-| If the signs are not equal check if a >= 0
- tstl d6
- bpl Lcmpsf$a$gt$b | if (a >= 0 && b < 0) => a > b
- bmi Lcmpsf$b$gt$a | if (a < 0 && b >= 0) => a < b
-1:
-| If the signs are equal check for < 0
- tstl d6
- bpl 1f
-| If both are negative exchange them
-#ifndef __mcoldfire__
- exg d0,d1
-#else
- movel d0,d7
- movel d1,d0
- movel d7,d1
-#endif
-1:
-| Now that they are positive we just compare them as longs (does this also
-| work for denormalized numbers?).
- cmpl d0,d1
- bhi Lcmpsf$b$gt$a | |b| > |a|
- bne Lcmpsf$a$gt$b | |b| < |a|
-| If we got here a == b.
- movel IMM (EQUAL),d0
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | put back the registers
-#else
- moveml sp@,d2-d7
-#endif
- unlk a6
- rts
-Lcmpsf$a$gt$b:
- movel IMM (GREATER),d0
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | put back the registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-Lcmpsf$b$gt$a:
- movel IMM (LESS),d0
-#ifndef __mcoldfire__
- moveml sp@+,d2-d7 | put back the registers
-#else
- moveml sp@,d2-d7
- | XXX if frame pointer is ever removed, stack pointer must
- | be adjusted here.
-#endif
- unlk a6
- rts
-
-Lcmpsf$a$0:
- bclr IMM (31),d6
- bra Lcmpsf$1
-Lcmpsf$b$0:
- bclr IMM (31),d7
- bra Lcmpsf$2
-
-|=============================================================================
-| rounding routines
-|=============================================================================
-
-| The rounding routines expect the number to be normalized in registers
-| d0-d1, with the exponent in register d2. They assume that the
-| exponent is larger or equal to 1. They return a properly normalized number
-| if possible, and a denormalized number otherwise. The exponent is returned
-| in d2.
-
-Lround$to$nearest:
-| We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
-| Here we assume that the exponent is not too small (this should be checked
-| before entering the rounding routine), but the number could be denormalized.
-
-| Check for denormalized numbers:
-1: btst IMM (FLT_MANT_DIG),d0
- bne 2f | if set the number is normalized
-| Normalize shifting left until bit #FLT_MANT_DIG is set or the exponent
-| is one (remember that a denormalized number corresponds to an
-| exponent of -F_BIAS+1).
-#ifndef __mcoldfire__
- cmpw IMM (1),d2 | remember that the exponent is at least one
-#else
- cmpl IMM (1),d2 | remember that the exponent is at least one
-#endif
- beq 2f | an exponent of one means denormalized
- addl d1,d1 | else shift and adjust the exponent
- addxl d0,d0 |
-#ifndef __mcoldfire__
- dbra d2,1b |
-#else
- subql IMM (1),d2
- bpl 1b
-#endif
-2:
-| Now round: we do it as follows: after the shifting we can write the
-| fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
-| If delta < 1, do nothing. If delta > 1, add 1 to f.
-| If delta == 1, we make sure the rounded number will be even (odd?)
-| (after shifting).
- btst IMM (0),d0 | is delta < 1?
- beq 2f | if so, do not do anything
- tstl d1 | is delta == 1?
- bne 1f | if so round to even
- movel d0,d1 |
- andl IMM (2),d1 | bit 1 is the last significant bit
- addl d1,d0 |
- bra 2f |
-1: movel IMM (1),d1 | else add 1
- addl d1,d0 |
-| Shift right once (because we used bit #FLT_MANT_DIG!).
-2: lsrl IMM (1),d0
-| Now check again bit #FLT_MANT_DIG (rounding could have produced a
-| 'fraction overflow' ...).
- btst IMM (FLT_MANT_DIG),d0
- beq 1f
- lsrl IMM (1),d0
-#ifndef __mcoldfire__
- addw IMM (1),d2
-#else
- addql IMM (1),d2
-#endif
-1:
-| If bit #FLT_MANT_DIG-1 is clear we have a denormalized number, so we
-| have to put the exponent to zero and return a denormalized number.
- btst IMM (FLT_MANT_DIG-1),d0
- beq 1f
- jmp a0@
-1: movel IMM (0),d2
- jmp a0@
-
-Lround$to$zero:
-Lround$to$plus:
-Lround$to$minus:
- jmp a0@
-#endif /* L_float */
-
-| gcc expects the routines __eqdf2, __nedf2, __gtdf2, __gedf2,
-| __ledf2, __ltdf2 to all return the same value as a direct call to
-| __cmpdf2 would. In this implementation, each of these routines
-| simply calls __cmpdf2. It would be more efficient to give the
-| __cmpdf2 routine several names, but separating them out will make it
-| easier to write efficient versions of these routines someday.
-
-#ifdef L_eqdf2
- .text
- .proc
- .globl SYM (__eqdf2)
-SYM (__eqdf2):
- link a6,IMM (0)
- movl a6@(20),sp@-
- movl a6@(16),sp@-
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpdf2)
- unlk a6
- rts
-#endif /* L_eqdf2 */
-
-#ifdef L_nedf2
- .text
- .proc
- .globl SYM (__nedf2)
-SYM (__nedf2):
- link a6,IMM (0)
- movl a6@(20),sp@-
- movl a6@(16),sp@-
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpdf2)
- unlk a6
- rts
-#endif /* L_nedf2 */
-
-#ifdef L_gtdf2
- .text
- .proc
- .globl SYM (__gtdf2)
-SYM (__gtdf2):
- link a6,IMM (0)
- movl a6@(20),sp@-
- movl a6@(16),sp@-
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpdf2)
- unlk a6
- rts
-#endif /* L_gtdf2 */
-
-#ifdef L_gedf2
- .text
- .proc
- .globl SYM (__gedf2)
-SYM (__gedf2):
- link a6,IMM (0)
- movl a6@(20),sp@-
- movl a6@(16),sp@-
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpdf2)
- unlk a6
- rts
-#endif /* L_gedf2 */
-
-#ifdef L_ltdf2
- .text
- .proc
- .globl SYM (__ltdf2)
-SYM (__ltdf2):
- link a6,IMM (0)
- movl a6@(20),sp@-
- movl a6@(16),sp@-
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpdf2)
- unlk a6
- rts
-#endif /* L_ltdf2 */
-
-#ifdef L_ledf2
- .text
- .proc
- .globl SYM (__ledf2)
-SYM (__ledf2):
- link a6,IMM (0)
- movl a6@(20),sp@-
- movl a6@(16),sp@-
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpdf2)
- unlk a6
- rts
-#endif /* L_ledf2 */
-
-| The comments above about __eqdf2, et. al., also apply to __eqsf2,
-| et. al., except that the latter call __cmpsf2 rather than __cmpdf2.
-
-#ifdef L_eqsf2
- .text
- .proc
- .globl SYM (__eqsf2)
-SYM (__eqsf2):
- link a6,IMM (0)
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpsf2)
- unlk a6
- rts
-#endif /* L_eqsf2 */
-
-#ifdef L_nesf2
- .text
- .proc
- .globl SYM (__nesf2)
-SYM (__nesf2):
- link a6,IMM (0)
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpsf2)
- unlk a6
- rts
-#endif /* L_nesf2 */
-
-#ifdef L_gtsf2
- .text
- .proc
- .globl SYM (__gtsf2)
-SYM (__gtsf2):
- link a6,IMM (0)
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpsf2)
- unlk a6
- rts
-#endif /* L_gtsf2 */
-
-#ifdef L_gesf2
- .text
- .proc
- .globl SYM (__gesf2)
-SYM (__gesf2):
- link a6,IMM (0)
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpsf2)
- unlk a6
- rts
-#endif /* L_gesf2 */
-
-#ifdef L_ltsf2
- .text
- .proc
- .globl SYM (__ltsf2)
-SYM (__ltsf2):
- link a6,IMM (0)
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpsf2)
- unlk a6
- rts
-#endif /* L_ltsf2 */
-
-#ifdef L_lesf2
- .text
- .proc
- .globl SYM (__lesf2)
-SYM (__lesf2):
- link a6,IMM (0)
- movl a6@(12),sp@-
- movl a6@(8),sp@-
- PICCALL SYM (__cmpsf2)
- unlk a6
- rts
-#endif /* L_lesf2 */
diff -urNa gcc-3.4.0/gcc/config/m68k/m68k.h gcc-3.4.0.mod/gcc/config/m68k/m68k.h
--- gcc-3.4.0/gcc/config/m68k/m68k.h 2005-12-18 09:48:13.000000000 +0100
+++ gcc-3.4.0.mod/gcc/config/m68k/m68k.h 2005-12-17 17:22:48.000000000 +0100
@@ -386,7 +386,7 @@
/* target machine storage layout */
/* Define for XFmode extended real floating point support. */
-#define LONG_DOUBLE_TYPE_SIZE 96
+#define LONG_DOUBLE_TYPE_SIZE 64
/* Set the value of FLT_EVAL_METHOD in float.h. When using 68040 fp
instructions, we get proper intermediate rounding, otherwise we
diff -urNa gcc-3.4.0/gcc/config/m68k/t-uclinux gcc-3.4.0.mod/gcc/config/m68k/t-uclinux
--- gcc-3.4.0/gcc/config/m68k/t-uclinux 2003-10-11 06:23:20.000000000 +0200
+++ gcc-3.4.0.mod/gcc/config/m68k/t-uclinux 2005-12-18 09:40:35.000000000 +0100
@@ -1,20 +1,25 @@
LIB1ASMSRC = m68k/lb1sf68.asm
-LIB1ASMFUNCS = _mulsi3 _udivsi3 _divsi3 _umodsi3 _modsi3 \
- _double _float _floatex \
- _eqdf2 _nedf2 _gtdf2 _gedf2 _ltdf2 _ledf2 \
- _eqsf2 _nesf2 _gtsf2 _gesf2 _ltsf2 _lesf2
-
-LIB2FUNCS_EXTRA = fpgnulib.c xfgnulib.c
-
-fpgnulib.c: $(srcdir)/config/m68k/fpgnulib.c
- cp $(srcdir)/config/m68k/fpgnulib.c fpgnulib.c
-xfgnulib.c: $(srcdir)/config/m68k/fpgnulib.c
- echo '#define EXTFLOAT' > xfgnulib.c
- cat $(srcdir)/config/m68k/fpgnulib.c >> xfgnulib.c
+LIB1ASMFUNCS = _mulsi3 _udivsi3 _divsi3 _umodsi3 _modsi3
+
+
+TARGET_LIBGCC2_CFLAGS = -msoft-float
+# We want fine grained libraries, so use the new code to build the
+# floating point emulation libraries.
+FPBIT = fp-bit.c
+DPBIT = dp-bit.c
+
+
+dp-bit.c: $(srcdir)/config/fp-bit.c
+ cat $(srcdir)/config/fp-bit.c >> dp-bit.c
+
+fp-bit.c: $(srcdir)/config/fp-bit.c
+ echo '#define FLOAT' > fp-bit.c
+ cat $(srcdir)/config/fp-bit.c >> fp-bit.c
+
MULTILIB_OPTIONS = m68000/m5200/m5206e/m528x/m5307/m5407/mcpu32 msep-data/mid-shared-library
MULTILIB_DIRNAMES =
-MULTILIB_MATCHES = m68000=mc68000 m68000=m68302 mcpu32=m68332 m5206e=m5272
+MULTILIB_MATCHES = m68000=mc68000 m68000=m68302 mcpu32=m68332 m5206e=m5272
MULTILIB_EXCEPTIONS = m68000/msep-data* m68000/mid-shared-library* msep-data* mid-shared-library*
LIBGCC = stmp-multilib
diff -urNa gcc-3.4.0/libiberty/strsignal.c gcc-3.4.0.mod/libiberty/strsignal.c
--- gcc-3.4.0/libiberty/strsignal.c 2003-04-15 22:36:33.000000000 +0200
+++ gcc-3.4.0.mod/libiberty/strsignal.c 2005-12-17 13:05:10.000000000 +0100
@@ -557,9 +557,10 @@
#ifndef HAVE_PSIGNAL
void
-psignal (signo, message)
+psignal (int signo, const char * message)
+/*psignal (signo, message)
unsigned signo;
- char *message;
+ char *message;*/
{
if (signal_names == NULL)
{
/* A C version of Kahan's Floating Point Test "Paranoia"
Thos Sumner, UCSF, Feb. 1985
David Gay, BTL, Jan. 1986
This is a rewrite from the Pascal version by
B. A. Wichmann, 18 Jan. 1985
(and does NOT exhibit good C programming style).
Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
compile with -DKR_headers or insert
#define KR_headers
at the beginning if you have an old-style C compiler.
(C) Apr 19 1983 in BASIC version by:
Professor W. M. Kahan,
567 Evans Hall
Electrical Engineering & Computer Science Dept.
University of California
Berkeley, California 94720
USA
converted to Pascal by:
B. A. Wichmann
National Physical Laboratory
Teddington Middx
TW11 OLW
UK
converted to C by:
David M. Gay and Thos Sumner
AT&T Bell Labs Computer Center, Rm. U-76
600 Mountain Avenue University of California
Murray Hill, NJ 07974 San Francisco, CA 94143
USA USA
with simultaneous corrections to the Pascal source (reflected
in the Pascal source available over netlib).
[A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
Reports of results on various systems from all the versions
of Paranoia are being collected by Richard Karpinski at the
same address as Thos Sumner. This includes sample outputs,
bug reports, and criticisms.
You may copy this program freely if you acknowledge its source.
Comments on the Pascal version to NPL, please.
The C version catches signals from floating-point exceptions.
If signal(SIGFPE,...) is unavailable in your environment, you may
#define NOSIGNAL to comment out the invocations of signal.
This source file is too big for some C compilers, but may be split
into pieces. Comments containing "SPLIT" suggest convenient places
for this splitting. At the end of these comments is an "ed script"
(for the UNIX(tm) editor ed) that will do this splitting.
By #defining Single when you compile this source, you may obtain
a single-precision C version of Paranoia.
The following is from the introductory commentary from Wichmann's work:
The BASIC program of Kahan is written in Microsoft BASIC using many
facilities which have no exact analogy in Pascal. The Pascal
version below cannot therefore be exactly the same. Rather than be
a minimal transcription of the BASIC program, the Pascal coding
follows the conventional style of block-structured languages. Hence
the Pascal version could be useful in producing versions in other
structured languages.
Rather than use identifiers of minimal length (which therefore have
little mnemonic significance), the Pascal version uses meaningful
identifiers as follows [Note: A few changes have been made for C]:
BASIC C BASIC C BASIC C
A J S StickyBit
A1 AInverse J0 NoErrors T
B Radix [Failure] T0 Underflow
B1 BInverse J1 NoErrors T2 ThirtyTwo
B2 RadixD2 [SeriousDefect] T5 OneAndHalf
B9 BMinusU2 J2 NoErrors T7 TwentySeven
C [Defect] T8 TwoForty
C1 CInverse J3 NoErrors U OneUlp
D [Flaw] U0 UnderflowThreshold
D4 FourD K PageNo U1
E0 L Milestone U2
E1 M V
E2 Exp2 N V0
E3 N1 V8
E5 MinSqEr O Zero V9
E6 SqEr O1 One W
E7 MaxSqEr O2 Two X
E8 O3 Three X1
E9 O4 Four X8
F1 MinusOne O5 Five X9 Random1
F2 Half O8 Eight Y
F3 Third O9 Nine Y1
F6 P Precision Y2
F9 Q Y9 Random2
G1 GMult Q8 Z
G2 GDiv Q9 Z0 PseudoZero
G3 GAddSub R Z1
H R1 RMult Z2
H1 HInverse R2 RDiv Z9
I R3 RAddSub
IO NoTrials R4 RSqrt
I3 IEEE R9 Random9
SqRWrng
All the variables in BASIC are true variables and in consequence,
the program is more difficult to follow since the "constants" must
be determined (the glossary is very helpful). The Pascal version
uses Real constants, but checks are added to ensure that the values
are correctly converted by the compiler.
The major textual change to the Pascal version apart from the
identifiersis that named procedures are used, inserting parameters
wherehelpful. New procedures are also introduced. The
correspondence is as follows:
BASIC Pascal
lines
90- 140 Pause
170- 250 Instructions
380- 460 Heading
480- 670 Characteristics
690- 870 History
2940-2950 Random
3710-3740 NewD
4040-4080 DoesYequalX
4090-4110 PrintIfNPositive
4640-4850 TestPartialUnderflow
=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
Below is an "ed script" that splits para.c into 10 files
of the form part[1-8].c, subs.c, and msgs.c, plus a header
file, paranoia.h, that these files require.
r paranoia.c
$
?SPLIT
.d
+d
-,$w msgs.c
-,$d
?SPLIT
.d
+d
-,$w subs.c
-,$d
?part8
+d
?include
.,$w part8.c
.,$d
-d
?part7
+d
?include
.,$w part7.c
.,$d
-d
?part6
+d
?include
.,$w part6.c
.,$d
-d
?part5
+d
?include
.,$w part5.c
.,$d
-d
?part4
+d
?include
.,$w part4.c
.,$d
-d
?part3
+d
?include
.,$w part3.c
.,$d
-d
?part2
+d
?include
.,$w part2.c
.,$d
?SPLIT
.d
1,/^#include/-1d
1,$w part1.c
/Computed constants/,$d
1,$s/^int/extern &/
1,$s/^FLOAT/extern &/
1,$s/^char/extern &/
1,$s! = .*!;!
/^Guard/,/^Round/s/^/extern /
/^jmp_buf/s/^/extern /
/^Sig_type/s/^/extern /
s/$/\
extern void sigfpe(INT);/
w paranoia.h
q
*/
#include <stdio.h>
#ifndef NOSIGNAL
#include <signal.h>
#endif
#include <setjmp.h>
#ifdef Single
#define FLOAT float
#define FABS(x) (float)fabs((double)(x))
#define FLOOR(x) (float)floor((double)(x))
#define LOG(x) (float)log((double)(x))
#define POW(x,y) (float)pow((double)(x),(double)(y))
#define SQRT(x) (float)sqrt((double)(x))
#else
#define FLOAT double
#define FABS(x) fabs(x)
#define FLOOR(x) floor(x)
#define LOG(x) log(x)
#define POW(x,y) pow(x,y)
#define SQRT(x) sqrt(x)
#endif
jmp_buf ovfl_buf;
#ifdef KR_headers
#define VOID /* void */
#define INT /* int */
#define FP /* FLOAT */
#define CHARP /* char * */
#define CHARPP /* char ** */
extern double fabs(), floor(), log(), pow(), sqrt();
extern void exit();
typedef void (*Sig_type)();
FLOAT Sign(), Random();
extern void BadCond();
extern void SqXMinX();
extern void TstCond();
extern void notify();
extern int read();
#else
#define VOID void
#define INT int
#define FP FLOAT
#define CHARP char *
#define CHARPP char **
#ifdef __STDC__
#include <stdlib.h>
#include <math.h>
#else
#ifdef __cplusplus
extern "C" {
#endif
extern double fabs(double), floor(double), log(double);
extern double pow(double,double), sqrt(double);
extern void exit(INT);
#ifdef __cplusplus
}
#endif
#endif
typedef void (*Sig_type)(int);
FLOAT Sign(FLOAT), Random(void);
extern void BadCond(int, char*);
extern void SqXMinX(int);
extern void TstCond(int, int, char*);
extern void notify(char*);
extern int read(int, char*, int);
#endif
#undef V9
extern void Characteristics(VOID);
extern void Heading(VOID);
extern void History(VOID);
extern void Instructions(VOID);
extern void IsYeqX(VOID);
extern void NewD(VOID);
extern void Pause(VOID);
extern void PrintIfNPositive(VOID);
extern void SR3750(VOID);
extern void SR3980(VOID);
extern void TstPtUf(VOID);
Sig_type sigsave;
#define KEYBOARD 0
FLOAT Radix, BInvrse, RadixD2, BMinusU2;
/*Small floating point constants.*/
FLOAT Zero = 0.0;
FLOAT Half = 0.5;
FLOAT One = 1.0;
FLOAT Two = 2.0;
FLOAT Three = 3.0;
FLOAT Four = 4.0;
FLOAT Five = 5.0;
FLOAT Eight = 8.0;
FLOAT Nine = 9.0;
FLOAT TwentySeven = 27.0;
FLOAT ThirtyTwo = 32.0;
FLOAT TwoForty = 240.0;
FLOAT MinusOne = -1.0;
FLOAT OneAndHalf = 1.5;
/*Integer constants*/
int NoTrials = 20; /*Number of tests for commutativity. */
#define False 0
#define True 1
/* Definitions for declared types
Guard == (Yes, No);
Rounding == (Chopped, Rounded, Other);
Message == packed array [1..40] of char;
Class == (Flaw, Defect, Serious, Failure);
*/
#define Yes 1
#define No 0
#define Chopped 2
#define Rounded 1
#define Other 0
#define Flaw 3
#define Defect 2
#define Serious 1
#define Failure 0
typedef int Guard, Rounding, Class;
typedef char Message;
/* Declarations of Variables */
int Indx;
char ch[8];
FLOAT AInvrse, A1;
FLOAT C, CInvrse;
FLOAT D, FourD;
FLOAT E0, E1, Exp2, E3, MinSqEr;
FLOAT SqEr, MaxSqEr, E9;
FLOAT Third;
FLOAT F6, F9;
FLOAT H, HInvrse;
int I;
FLOAT StickyBit, J;
FLOAT MyZero;
FLOAT Precision;
FLOAT Q, Q9;
FLOAT R, Random9;
FLOAT T, Underflow, S;
FLOAT OneUlp, UfThold, U1, U2;
FLOAT V, V0, V9;
FLOAT W;
FLOAT X, X1, X2, X8, Random1;
FLOAT Y, Y1, Y2, Random2;
FLOAT Z, PseudoZero, Z1, Z2, Z9;
int ErrCnt[4];
int fpecount;
int Milestone;
int PageNo;
int M, N, N1;
Guard GMult, GDiv, GAddSub;
Rounding RMult, RDiv, RAddSub, RSqrt;
int Break, Done, NotMonot, Monot, Anomaly, IEEE,
SqRWrng, UfNGrad;
/* Computed constants. */
/*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
/*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
/* floating point exception receiver */
void
sigfpe(INT x)
{
fpecount++;
printf("\n* * * FLOATING-POINT ERROR %d * * *\n", x);
fflush(stdout);
if (sigsave) {
#ifndef NOSIGNAL
signal(SIGFPE, sigsave);
#endif
sigsave = 0;
longjmp(ovfl_buf, 1);
}
exit(1);
}
main(VOID)
{
/* First two assignments use integer right-hand sides. */
Zero = 0;
One = 1;
Two = One + One;
Three = Two + One;
Four = Three + One;
Five = Four + One;
Eight = Four + Four;
Nine = Three * Three;
TwentySeven = Nine * Three;
ThirtyTwo = Four * Eight;
TwoForty = Four * Five * Three * Four;
MinusOne = -One;
Half = One / Two;
OneAndHalf = One + Half;
ErrCnt[Failure] = 0;
ErrCnt[Serious] = 0;
ErrCnt[Defect] = 0;
ErrCnt[Flaw] = 0;
PageNo = 1;
/*=============================================*/
Milestone = 0;
/*=============================================*/
#ifndef NOSIGNAL
signal(SIGFPE, sigfpe);
#endif
Instructions();
Pause();
Heading();
Pause();
Characteristics();
Pause();
History();
Pause();
/*=============================================*/
Milestone = 7;
/*=============================================*/
printf("Program is now RUNNING tests on small integers:\n");
TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
&& (One > Zero) && (One + One == Two),
"0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
Z = - Zero;
if (Z != 0.0) {
ErrCnt[Failure] = ErrCnt[Failure] + 1;
printf("Comparison alleges that -0.0 is Non-zero!\n");
U2 = 0.001;
Radix = 1;
TstPtUf();
}
TstCond (Failure, (Three == Two + One) && (Four == Three + One)
&& (Four + Two * (- Two) == Zero)
&& (Four - Three - One == Zero),
"3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
TstCond (Failure, (MinusOne == (0 - One))
&& (MinusOne + One == Zero ) && (One + MinusOne == Zero)
&& (MinusOne + FABS(One) == Zero)
&& (MinusOne + MinusOne * MinusOne == Zero),
"-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
TstCond (Failure, Half + MinusOne + Half == Zero,
"1/2 + (-1) + 1/2 != 0");
/*=============================================*/
/*SPLIT
{
extern void part2(VOID), part3(VOID), part4(VOID),
part5(VOID), part6(VOID), part7(VOID);
int part8(VOID);
part2();
part3();
part4();
part5();
part6();
part7();
return part8();
}
}
#include "paranoia.h"
void part2(VOID){
*/
Milestone = 10;
/*=============================================*/
TstCond (Failure, (Nine == Three * Three)
&& (TwentySeven == Nine * Three) && (Eight == Four + Four)
&& (ThirtyTwo == Eight * Four)
&& (ThirtyTwo - TwentySeven - Four - One == Zero),
"9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
TstCond (Failure, (Five == Four + One) &&
(TwoForty == Four * Five * Three * Four)
&& (TwoForty / Three - Four * Four * Five == Zero)
&& ( TwoForty / Four - Five * Three * Four == Zero)
&& ( TwoForty / Five - Four * Three * Four == Zero),
"5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
if (ErrCnt[Failure] == 0) {
printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
printf("\n");
}
printf("Searching for Radix and Precision.\n");
W = One;
do {
W = W + W;
Y = W + One;
Z = Y - W;
Y = Z - One;
} while (MinusOne + FABS(Y) < Zero);
/*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
Precision = Zero;
Y = One;
do {
Radix = W + Y;
Y = Y + Y;
Radix = Radix - W;
} while ( Radix == Zero);
if (Radix < Two) Radix = One;
printf("Radix = %f .\n", Radix);
if (Radix != 1) {
W = One;
do {
Precision = Precision + One;
W = W * Radix;
Y = W + One;
} while ((Y - W) == One);
}
/*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
...*/
U1 = One / W;
U2 = Radix * U1;
printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
printf("Recalculating radix and precision\n ");
/*save old values*/
E0 = Radix;
E1 = U1;
E9 = U2;
E3 = Precision;
X = Four / Three;
Third = X - One;
F6 = Half - Third;
X = F6 + F6;
X = FABS(X - Third);
if (X < U2) X = U2;
/*... now X = (unknown no.) ulps of 1+...*/
do {
U2 = X;
Y = Half * U2 + ThirtyTwo * U2 * U2;
Y = One + Y;
X = Y - One;
} while ( ! ((U2 <= X) || (X <= Zero)));
/*... now U2 == 1 ulp of 1 + ... */
X = Two / Three;
F6 = X - Half;
Third = F6 + F6;
X = Third - Half;
X = FABS(X + F6);
if (X < U1) X = U1;
/*... now X == (unknown no.) ulps of 1 -... */
do {
U1 = X;
Y = Half * U1 + ThirtyTwo * U1 * U1;
Y = Half - Y;
X = Half + Y;
Y = Half - X;
X = Half + Y;
} while ( ! ((U1 <= X) || (X <= Zero)));
/*... now U1 == 1 ulp of 1 - ... */
if (U1 == E1) printf("confirms closest relative separation U1 .\n");
else printf("gets better closest relative separation U1 = %.7e .\n", U1);
W = One / U1;
F9 = (Half - U1) + Half;
Radix = FLOOR(0.01 + U2 / U1);
if (Radix == E0) printf("Radix confirmed.\n");
else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
TstCond (Defect, Radix <= Eight + Eight,
"Radix is too big: roundoff problems");
TstCond (Flaw, (Radix == Two) || (Radix == 10)
|| (Radix == One), "Radix is not as good as 2 or 10");
/*=============================================*/
Milestone = 20;
/*=============================================*/
TstCond (Failure, F9 - Half < Half,
"(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
X = F9;
I = 1;
Y = X - Half;
Z = Y - Half;
TstCond (Failure, (X != One)
|| (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
X = One + U2;
I = 0;
/*=============================================*/
Milestone = 25;
/*=============================================*/
/*... BMinusU2 = nextafter(Radix, 0) */
BMinusU2 = Radix - One;
BMinusU2 = (BMinusU2 - U2) + One;
/* Purify Integers */
if (Radix != One) {
X = - TwoForty * LOG(U1) / LOG(Radix);
Y = FLOOR(Half + X);
if (FABS(X - Y) * Four < One) X = Y;
Precision = X / TwoForty;
Y = FLOOR(Half + Precision);
if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
}
if ((Precision != FLOOR(Precision)) || (Radix == One)) {
printf("Precision cannot be characterized by an Integer number\n");
printf("of significant digits but, by itself, this is a minor flaw.\n");
}
if (Radix == One)
printf("logarithmic encoding has precision characterized solely by U1.\n");
else printf("The number of significant digits of the Radix is %f .\n",
Precision);
TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
"Precision worse than 5 decimal figures ");
/*=============================================*/
Milestone = 30;
/*=============================================*/
/* Test for extra-precise subepressions */
X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
do {
Z2 = X;
X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
} while ( ! ((Z2 <= X) || (X <= Zero)));
X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
do {
Z1 = Z;
Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
+ One / Two)) + One / Two;
} while ( ! ((Z1 <= Z) || (Z <= Zero)));
do {
do {
Y1 = Y;
Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
)) + Half;
} while ( ! ((Y1 <= Y) || (Y <= Zero)));
X1 = X;
X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
} while ( ! ((X1 <= X) || (X <= Zero)));
if ((X1 != Y1) || (X1 != Z1)) {
BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
printf("are symptoms of inconsistencies introduced\n");
printf("by extra-precise evaluation of arithmetic subexpressions.\n");
notify("Possibly some part of this");
if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
"That feature is not tested further by this program.\n") ;
}
else {
if ((Z1 != U1) || (Z2 != U2)) {
if ((Z1 >= U1) || (Z2 >= U2)) {
BadCond(Failure, "");
notify("Precision");
printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
}
else {
if ((Z1 <= Zero) || (Z2 <= Zero)) {
printf("Because of unusual Radix = %f", Radix);
printf(", or exact rational arithmetic a result\n");
printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
notify("of an\nextra-precision");
}
if (Z1 != Z2 || Z1 > Zero) {
X = Z1 / U1;
Y = Z2 / U2;
if (Y > X) X = Y;
Q = - LOG(X);
printf("Some subexpressions appear to be calculated extra\n");
printf("precisely with about %g extra B-digits, i.e.\n",
(Q / LOG(Radix)));
printf("roughly %g extra significant decimals.\n",
Q / LOG(10.));
}
printf("That feature is not tested further by this program.\n");
}
}
}
Pause();
/*=============================================*/
/*SPLIT
}
#include "paranoia.h"
void part3(VOID){
*/
Milestone = 35;
/*=============================================*/
if (Radix >= Two) {
X = W / (Radix * Radix);
Y = X + One;
Z = Y - X;
T = Z + U2;
X = T - Z;
TstCond (Failure, X == U2,
"Subtraction is not normalized X=Y,X+Z != Y+Z!");
if (X == U2) printf(
"Subtraction appears to be normalized, as it should be.");
}
printf("\nChecking for guard digit in *, /, and -.\n");
Y = F9 * One;
Z = One * F9;
X = F9 - Half;
Y = (Y - Half) - X;
Z = (Z - Half) - X;
X = One + U2;
T = X * Radix;
R = Radix * X;
X = T - Radix;
X = X - Radix * U2;
T = R - Radix;
T = T - Radix * U2;
X = X * (Radix - One);
T = T * (Radix - One);
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
else {
GMult = No;
TstCond (Serious, False,
"* lacks a Guard Digit, so 1*X != X");
}
Z = Radix * U2;
X = One + Z;
Y = FABS((X + Z) - X * X) - U2;
X = One - U2;
Z = FABS((X - U2) - X * X) - U1;
TstCond (Failure, (Y <= Zero)
&& (Z <= Zero), "* gets too many final digits wrong.\n");
Y = One - U2;
X = One + U2;
Z = One / Y;
Y = Z - X;
X = One / Three;
Z = Three / Nine;
X = X - Z;
T = Nine / TwentySeven;
Z = Z - T;
TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
"Division lacks a Guard Digit, so error can exceed 1 ulp\n\
or 1/3 and 3/9 and 9/27 may disagree");
Y = F9 / One;
X = F9 - Half;
Y = (Y - Half) - X;
X = One + U2;
T = X / One;
X = T - X;
if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
else {
GDiv = No;
TstCond (Serious, False,
"Division lacks a Guard Digit, so X/1 != X");
}
X = One / (One + U2);
Y = X - Half - Half;
TstCond (Serious, Y < Zero,
"Computed value of 1/1.000..1 >= 1");
X = One - U2;
Y = One + Radix * U2;
Z = X * Radix;
T = Y * Radix;
R = Z / Radix;
StickyBit = T / Radix;
X = R - X;
Y = StickyBit - Y;
TstCond (Failure, X == Zero && Y == Zero,
"* and/or / gets too many last digits wrong");
Y = One - U1;
X = One - F9;
Y = One - Y;
T = Radix - U2;
Z = Radix - BMinusU2;
T = Radix - T;
if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
else {
GAddSub = No;
TstCond (Serious, False,
"- lacks Guard Digit, so cancellation is obscured");
}
if (F9 != One && F9 - One >= Zero) {
BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
printf(" such precautions against division by zero as\n");
printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
}
if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
" *, /, and - appear to have guard digits, as they should.\n");
/*=============================================*/
Milestone = 40;
/*=============================================*/
Pause();
printf("Checking rounding on multiply, divide and add/subtract.\n");
RMult = Other;
RDiv = Other;
RAddSub = Other;
RadixD2 = Radix / Two;
A1 = Two;
Done = False;
do {
AInvrse = Radix;
do {
X = AInvrse;
AInvrse = AInvrse / A1;
} while ( ! (FLOOR(AInvrse) != AInvrse));
Done = (X == One) || (A1 > Three);
if (! Done) A1 = Nine + One;
} while ( ! (Done));
if (X == One) A1 = Radix;
AInvrse = One / A1;
X = A1;
Y = AInvrse;
Done = False;
do {
Z = X * Y - Half;
TstCond (Failure, Z == Half,
"X * (1/X) differs from 1");
Done = X == Radix;
X = Radix;
Y = One / X;
} while ( ! (Done));
Y2 = One + U2;
Y1 = One - U2;
X = OneAndHalf - U2;
Y = OneAndHalf + U2;
Z = (X - U2) * Y2;
T = Y * Y1;
Z = Z - X;
T = T - X;
X = X * Y2;
Y = (Y + U2) * Y1;
X = X - OneAndHalf;
Y = Y - OneAndHalf;
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
X = (OneAndHalf + U2) * Y2;
Y = OneAndHalf - U2 - U2;
Z = OneAndHalf + U2 + U2;
T = (OneAndHalf - U2) * Y1;
X = X - (Z + U2);
StickyBit = Y * Y1;
S = Z * Y2;
T = T - Y;
Y = (U2 - Y) + StickyBit;
Z = S - (Z + U2 + U2);
StickyBit = (Y2 + U2) * Y1;
Y1 = Y2 * Y1;
StickyBit = StickyBit - Y2;
Y1 = Y1 - Half;
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
&& ( StickyBit == Zero) && (Y1 == Half)) {
RMult = Rounded;
printf("Multiplication appears to round correctly.\n");
}
else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
&& (T < Zero) && (StickyBit + U2 == Zero)
&& (Y1 < Half)) {
RMult = Chopped;
printf("Multiplication appears to chop.\n");
}
else printf("* is neither chopped nor correctly rounded.\n");
if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
}
else printf("* is neither chopped nor correctly rounded.\n");
/*=============================================*/
Milestone = 45;
/*=============================================*/
Y2 = One + U2;
Y1 = One - U2;
Z = OneAndHalf + U2 + U2;
X = Z / Y2;
T = OneAndHalf - U2 - U2;
Y = (T - U2) / Y1;
Z = (Z + U2) / Y2;
X = X - OneAndHalf;
Y = Y - T;
T = T / Y1;
Z = Z - (OneAndHalf + U2);
T = (U2 - OneAndHalf) + T;
if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
X = OneAndHalf / Y2;
Y = OneAndHalf - U2;
Z = OneAndHalf + U2;
X = X - Y;
T = OneAndHalf / Y1;
Y = Y / Y1;
T = T - (Z + U2);
Y = Y - Z;
Z = Z / Y2;
Y1 = (Y2 + U2) / Y2;
Z = Z - OneAndHalf;
Y2 = Y1 - Y2;
Y1 = (F9 - U1) / F9;
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
&& (Y2 == Zero) && (Y2 == Zero)
&& (Y1 - Half == F9 - Half )) {
RDiv = Rounded;
printf("Division appears to round correctly.\n");
if (GDiv == No) notify("Division");
}
else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
&& (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
RDiv = Chopped;
printf("Division appears to chop.\n");
}
}
if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
BInvrse = One / Radix;
TstCond (Failure, (BInvrse * Radix - Half == Half),
"Radix * ( 1 / Radix ) differs from 1");
/*=============================================*/
/*SPLIT
}
#include "paranoia.h"
void part4(VOID){
*/
Milestone = 50;
/*=============================================*/
TstCond (Failure, ((F9 + U1) - Half == Half)
&& ((BMinusU2 + U2 ) - One == Radix - One),
"Incomplete carry-propagation in Addition");
X = One - U1 * U1;
Y = One + U2 * (One - U2);
Z = F9 - Half;
X = (X - Half) - Z;
Y = Y - One;
if ((X == Zero) && (Y == Zero)) {
RAddSub = Chopped;
printf("Add/Subtract appears to be chopped.\n");
}
if (GAddSub == Yes) {
X = (Half + U2) * U2;
Y = (Half - U2) * U2;
X = One + X;
Y = One + Y;
X = (One + U2) - X;
Y = One - Y;
if ((X == Zero) && (Y == Zero)) {
X = (Half + U2) * U1;
Y = (Half - U2) * U1;
X = One - X;
Y = One - Y;
X = F9 - X;
Y = One - Y;
if ((X == Zero) && (Y == Zero)) {
RAddSub = Rounded;
printf("Addition/Subtraction appears to round correctly.\n");
if (GAddSub == No) notify("Add/Subtract");
}
else printf("Addition/Subtraction neither rounds nor chops.\n");
}
else printf("Addition/Subtraction neither rounds nor chops.\n");
}
else printf("Addition/Subtraction neither rounds nor chops.\n");
S = One;
X = One + Half * (One + Half);
Y = (One + U2) * Half;
Z = X - Y;
T = Y - X;
StickyBit = Z + T;
if (StickyBit != Zero) {
S = Zero;
BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
}
StickyBit = Zero;
if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
&& (RMult == Rounded) && (RDiv == Rounded)
&& (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
printf("Checking for sticky bit.\n");
X = (Half + U1) * U2;
Y = Half * U2;
Z = One + Y;
T = One + X;
if ((Z - One <= Zero) && (T - One >= U2)) {
Z = T + Y;
Y = Z - X;
if ((Z - T >= U2) && (Y - T == Zero)) {
X = (Half + U1) * U1;
Y = Half * U1;
Z = One - Y;
T = One - X;
if ((Z - One == Zero) && (T - F9 == Zero)) {
Z = (Half - U1) * U1;
T = F9 - Z;
Q = F9 - Y;
if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
Z = (One + U2) * OneAndHalf;
T = (OneAndHalf + U2) - Z + U2;
X = One + Half / Radix;
Y = One + Radix * U2;
Z = X * Y;
if (T == Zero && X + Radix * U2 - Z == Zero) {
if (Radix != Two) {
X = Two + U2;
Y = X / Two;
if ((Y - One == Zero)) StickyBit = S;
}
else StickyBit = S;
}
}
}
}
}
}
if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
else printf("Sticky bit used incorrectly or not at all.\n");
TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
RMult == Other || RDiv == Other || RAddSub == Other),
"lack(s) of guard digits or failure(s) to correctly round or chop\n\
(noted above) count as one flaw in the final tally below");
/*=============================================*/
Milestone = 60;
/*=============================================*/
printf("\n");
printf("Does Multiplication commute? ");
printf("Testing on %d random pairs.\n", NoTrials);
Random9 = SQRT(3.0);
Random1 = Third;
I = 1;
do {
X = Random();
Y = Random();
Z9 = Y * X;
Z = X * Y;
Z9 = Z - Z9;
I = I + 1;
} while ( ! ((I > NoTrials) || (Z9 != Zero)));
if (I == NoTrials) {
Random1 = One + Half / Three;
Random2 = (U2 + U1) + One;
Z = Random1 * Random2;
Y = Random2 * Random1;
Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
Three) * ((U2 + U1) + One);
}
if (! ((I == NoTrials) || (Z9 == Zero)))
BadCond(Defect, "X * Y == Y * X trial fails.\n");
else printf(" No failures found in %d integer pairs.\n", NoTrials);
/*=============================================*/
Milestone = 70;
/*=============================================*/
printf("\nRunning test of square root(x).\n");
TstCond (Failure, (Zero == SQRT(Zero))
&& (- Zero == SQRT(- Zero))
&& (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
MinSqEr = Zero;
MaxSqEr = Zero;
J = Zero;
X = Radix;
OneUlp = U2;
SqXMinX (Serious);
X = BInvrse;
OneUlp = BInvrse * U1;
SqXMinX (Serious);
X = U1;
OneUlp = U1 * U1;
SqXMinX (Serious);
if (J != Zero) Pause();
printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
J = Zero;
X = Two;
Y = Radix;
if ((Radix != One)) do {
X = Y;
Y = Radix * Y;
} while ( ! ((Y - X >= NoTrials)));
OneUlp = X * U2;
I = 1;
while (I <= NoTrials) {
X = X + One;
SqXMinX (Defect);
if (J > Zero) break;
I = I + 1;
}
printf("Test for sqrt monotonicity.\n");
I = - 1;
X = BMinusU2;
Y = Radix;
Z = Radix + Radix * U2;
NotMonot = False;
Monot = False;
while ( ! (NotMonot || Monot)) {
I = I + 1;
X = SQRT(X);
Q = SQRT(Y);
Z = SQRT(Z);
if ((X > Q) || (Q > Z)) NotMonot = True;
else {
Q = FLOOR(Q + Half);
if (!(I > 0 || Radix == Q * Q)) Monot = True;
else if (I > 0) {
if (I > 1) Monot = True;
else {
Y = Y * BInvrse;
X = Y - U1;
Z = Y + U1;
}
}
else {
Y = Q;
X = Y - U2;
Z = Y + U2;
}
}
}
if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
else {
BadCond(Defect, "");
printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
}
/*=============================================*/
/*SPLIT
}
#include "paranoia.h"
void part5(VOID){
*/
Milestone = 80;
/*=============================================*/
MinSqEr = MinSqEr + Half;
MaxSqEr = MaxSqEr - Half;
Y = (SQRT(One + U2) - One) / U2;
SqEr = (Y - One) + U2 / Eight;
if (SqEr > MaxSqEr) MaxSqEr = SqEr;
SqEr = Y + U2 / Eight;
if (SqEr < MinSqEr) MinSqEr = SqEr;
Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
SqEr = Y + U1 / Eight;
if (SqEr > MaxSqEr) MaxSqEr = SqEr;
SqEr = (Y + One) + U1 / Eight;
if (SqEr < MinSqEr) MinSqEr = SqEr;
OneUlp = U2;
X = OneUlp;
for( Indx = 1; Indx <= 3; ++Indx) {
Y = SQRT((X + U1 + X) + F9);
Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
SqEr = (Y + Half) + Z;
if (SqEr < MinSqEr) MinSqEr = SqEr;
SqEr = (Y - Half) + Z;
if (SqEr > MaxSqEr) MaxSqEr = SqEr;
if (((Indx == 1) || (Indx == 3)))
X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
else {
OneUlp = U1;
X = - OneUlp;
}
}
/*=============================================*/
Milestone = 85;
/*=============================================*/
SqRWrng = False;
Anomaly = False;
RSqrt = Other; /* ~dgh */
if (Radix != One) {
printf("Testing whether sqrt is rounded or chopped.\n");
D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
/* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
X = D / Radix;
Y = D / A1;
if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
Anomaly = True;
}
else {
X = Zero;
Z2 = X;
Y = One;
Y2 = Y;
Z1 = Radix - One;
FourD = Four * D;
do {
if (Y2 > Z2) {
Q = Radix;
Y1 = Y;
do {
X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
Q = Y1;
Y1 = X1;
} while ( ! (X1 <= Zero));
if (Q <= One) {
Z2 = Y2;
Z = Y;
}
}
Y = Y + Two;
X = X + Eight;
Y2 = Y2 + X;
if (Y2 >= FourD) Y2 = Y2 - FourD;
} while ( ! (Y >= D));
X8 = FourD - Z2;
Q = (X8 + Z * Z) / FourD;
X8 = X8 / Eight;
if (Q != FLOOR(Q)) Anomaly = True;
else {
Break = False;
do {
X = Z1 * Z;
X = X - FLOOR(X / Radix) * Radix;
if (X == One)
Break = True;
else
Z1 = Z1 - One;
} while ( ! (Break || (Z1 <= Zero)));
if ((Z1 <= Zero) && (! Break)) Anomaly = True;
else {
if (Z1 > RadixD2) Z1 = Z1 - Radix;
do {
NewD();
} while ( ! (U2 * D >= F9));
if (D * Radix - D != W - D) Anomaly = True;
else {
Z2 = D;
I = 0;
Y = D + (One + Z) * Half;
X = D + Z + Q;
SR3750();
Y = D + (One - Z) * Half + D;
X = D - Z + D;
X = X + Q + X;
SR3750();
NewD();
if (D - Z2 != W - Z2) Anomaly = True;
else {
Y = (D - Z2) + (Z2 + (One - Z) * Half);
X = (D - Z2) + (Z2 - Z + Q);
SR3750();
Y = (One + Z) * Half;
X = Q;
SR3750();
if (I == 0) Anomaly = True;
}
}
}
}
}
if ((I == 0) || Anomaly) {
BadCond(Failure, "Anomalous arithmetic with Integer < ");
printf("Radix^Precision = %.7e\n", W);
printf(" fails test whether sqrt rounds or chops.\n");
SqRWrng = True;
}
}
if (! Anomaly) {
if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
RSqrt = Rounded;
printf("Square root appears to be correctly rounded.\n");
}
else {
if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
|| (MinSqEr + Radix < Half)) SqRWrng = True;
else {
RSqrt = Chopped;
printf("Square root appears to be chopped.\n");
}
}
}
if (SqRWrng) {
printf("Square root is neither chopped nor correctly rounded.\n");
printf("Observed errors run from %.7e ", MinSqEr - Half);
printf("to %.7e ulps.\n", Half + MaxSqEr);
TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
"sqrt gets too many last digits wrong");
}
/*=============================================*/
Milestone = 90;
/*=============================================*/
Pause();
printf("Testing powers Z^i for small Integers Z and i.\n");
N = 0;
/* ... test powers of zero. */
I = 0;
Z = -Zero;
M = 3;
Break = False;
do {
X = One;
SR3980();
if (I <= 10) {
I = 1023;
SR3980();
}
if (Z == MinusOne) Break = True;
else {
Z = MinusOne;
/* .. if(-1)^N is invalid, replace MinusOne by One. */
I = - 4;
}
} while ( ! Break);
PrintIfNPositive();
N1 = N;
N = 0;
Z = A1;
M = (int)FLOOR(Two * LOG(W) / LOG(A1));
Break = False;
do {
X = Z;
I = 1;
SR3980();
if (Z == AInvrse) Break = True;
else Z = AInvrse;
} while ( ! (Break));
/*=============================================*/
Milestone = 100;
/*=============================================*/
/* Powers of Radix have been tested, */
/* next try a few primes */
M = NoTrials;
Z = Three;
do {
X = Z;
I = 1;
SR3980();
do {
Z = Z + Two;
} while ( Three * FLOOR(Z / Three) == Z );
} while ( Z < Eight * Three );
if (N > 0) {
printf("Errors like this may invalidate financial calculations\n");
printf("\tinvolving interest rates.\n");
}
PrintIfNPositive();
N += N1;
if (N == 0) printf("... no discrepancies found.\n");
if (N > 0) Pause();
else printf("\n");
/*=============================================*/
/*SPLIT
}
#include "paranoia.h"
void part6(VOID){
*/
Milestone = 110;
/*=============================================*/
printf("Seeking Underflow thresholds UfThold and E0.\n");
D = U1;
if (Precision != FLOOR(Precision)) {
D = BInvrse;
X = Precision;
do {
D = D * BInvrse;
X = X - One;
} while ( X > Zero);
}
Y = One;
Z = D;
/* ... D is power of 1/Radix < 1. */
do {
C = Y;
Y = Z;
Z = Y * Y;
} while ((Y > Z) && (Z + Z > Z));
Y = C;
Z = Y * D;
do {
C = Y;
Y = Z;
Z = Y * D;
} while ((Y > Z) && (Z + Z > Z));
if (Radix < Two) HInvrse = Two;
else HInvrse = Radix;
H = One / HInvrse;
/* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
CInvrse = One / C;
E0 = C;
Z = E0 * H;
/* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
do {
Y = E0;
E0 = Z;
Z = E0 * H;
} while ((E0 > Z) && (Z + Z > Z));
UfThold = E0;
E1 = Zero;
Q = Zero;
E9 = U2;
S = One + E9;
D = C * S;
if (D <= C) {
E9 = Radix * U2;
S = One + E9;
D = C * S;
if (D <= C) {
BadCond(Failure, "multiplication gets too many last digits wrong.\n");
Underflow = E0;
Y1 = Zero;
PseudoZero = Z;
Pause();
}
}
else {
Underflow = D;
PseudoZero = Underflow * H;
UfThold = Zero;
do {
Y1 = Underflow;
Underflow = PseudoZero;
if (E1 + E1 <= E1) {
Y2 = Underflow * HInvrse;
E1 = FABS(Y1 - Y2);
Q = Y1;
if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
}
PseudoZero = PseudoZero * H;
} while ((Underflow > PseudoZero)
&& (PseudoZero + PseudoZero > PseudoZero));
}
/* Comment line 4530 .. 4560 */
if (PseudoZero != Zero) {
printf("\n");
Z = PseudoZero;
/* ... Test PseudoZero for "phoney- zero" violates */
/* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
... */
if (PseudoZero <= Zero) {
BadCond(Failure, "Positive expressions can underflow to an\n");
printf("allegedly negative value\n");
printf("PseudoZero that prints out as: %g .\n", PseudoZero);
X = - PseudoZero;
if (X <= Zero) {
printf("But -PseudoZero, which should be\n");
printf("positive, isn't; it prints out as %g .\n", X);
}
}
else {
BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
printf("value PseudoZero that prints out as %g .\n", PseudoZero);
}
TstPtUf();
}
/*=============================================*/
Milestone = 120;
/*=============================================*/
if (CInvrse * Y > CInvrse * Y1) {
S = H * S;
E0 = Underflow;
}
if (! ((E1 == Zero) || (E1 == E0))) {
BadCond(Defect, "");
if (E1 < E0) {
printf("Products underflow at a higher");
printf(" threshold than differences.\n");
if (PseudoZero == Zero)
E0 = E1;
}
else {
printf("Difference underflows at a higher");
printf(" threshold than products.\n");
}
}
printf("Smallest strictly positive number found is E0 = %g .\n", E0);
Z = E0;
TstPtUf();
Underflow = E0;
if (N == 1) Underflow = Y;
I = 4;
if (E1 == Zero) I = 3;
if (UfThold == Zero) I = I - 2;
UfNGrad = True;
switch (I) {
case 1:
UfThold = Underflow;
if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
UfThold = Y;
BadCond(Failure, "Either accuracy deteriorates as numbers\n");
printf("approach a threshold = %.17e\n", UfThold);;
printf(" coming down from %.17e\n", C);
printf(" or else multiplication gets too many last digits wrong.\n");
}
Pause();
break;
case 2:
BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
printf("Q == Y while denying that |Q - Y| == 0; these values\n");
printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
UfThold = Q;
break;
case 3:
X = X;
break;
case 4:
if ((Q == UfThold) && (E1 == E0)
&& (FABS( UfThold - E1 / E9) <= E1)) {
UfNGrad = False;
printf("Underflow is gradual; it incurs Absolute Error =\n");
printf("(roundoff in UfThold) < E0.\n");
Y = E0 * CInvrse;
Y = Y * (OneAndHalf + U2);
X = CInvrse * (One + U2);
Y = Y / X;
IEEE = (Y == E0);
}
}
if (UfNGrad) {
printf("\n");
sigsave = sigfpe;
if (setjmp(ovfl_buf)) {
printf("Underflow / UfThold failed!\n");
R = H + H;
}
else R = SQRT(Underflow / UfThold);
sigsave = 0;
if (R <= H) {
Z = R * UfThold;
X = Z * (One + R * H * (One + H));
}
else {
Z = UfThold;
X = Z * (One + H * H * (One + H));
}
if (! ((X == Z) || (X - Z != Zero))) {
BadCond(Flaw, "");
printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
Z9 = X - Z;
printf("yet X - Z yields %.17e .\n", Z9);
printf(" Should this NOT signal Underflow, ");
printf("this is a SERIOUS DEFECT\nthat causes ");
printf("confusion when innocent statements like\n");;
printf(" if (X == Z) ... else");
printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
printf("encounter Division by Zero although actually\n");
sigsave = sigfpe;
if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
sigsave = 0;
}
}
printf("The Underflow threshold is %.17e, %s\n", UfThold,
" below which");
printf("calculation may suffer larger Relative error than ");
printf("merely roundoff.\n");
Y2 = U1 * U1;
Y = Y2 * Y2;
Y2 = Y * U1;
if (Y2 <= UfThold) {
if (Y > E0) {
BadCond(Defect, "");
I = 5;
}
else {
BadCond(Serious, "");
I = 4;
}
printf("Range is too narrow; U1^%d Underflows.\n", I);
}
/*=============================================*/
/*SPLIT
}
#include "paranoia.h"
void part7(VOID){
*/
Milestone = 130;
/*=============================================*/
Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
Y2 = Y + Y;
printf("Since underflow occurs below the threshold\n");
printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n",
HInvrse, Y2);
printf("actually calculating yields:");
if (setjmp(ovfl_buf)) {
sigsave = 0;
BadCond(Serious, "trap on underflow.\n");
}
else {
sigsave = sigfpe;
V9 = POW(HInvrse, Y2);
sigsave = 0;
printf(" %.17e .\n", V9);
if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
BadCond(Serious, "this is not between 0 and underflow\n");
printf(" threshold = %.17e .\n", UfThold);
}
else if (! (V9 > UfThold * (One + E9)))
printf("This computed value is O.K.\n");
else {
BadCond(Defect, "this is not between 0 and underflow\n");
printf(" threshold = %.17e .\n", UfThold);
}
}
/*=============================================*/
Milestone = 140;
/*=============================================*/
printf("\n");
/* ...calculate Exp2 == exp(2) == 7.389056099... */
X = Zero;
I = 2;
Y = Two * Three;
Q = Zero;
N = 0;
do {
Z = X;
I = I + 1;
Y = Y / (I + I);
R = Y + Q;
X = Z + R;
Q = (Z - X) + R;
} while(X > Z);
Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
X = Z * Z;
Exp2 = X * X;
X = F9;
Y = X - U1;
printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
Exp2);
for(I = 1;;) {
Z = X - BInvrse;
Z = (X + One) / (Z - (One - BInvrse));
Q = POW(X, Z) - Exp2;
if (FABS(Q) > TwoForty * U2) {
N = 1;
V9 = (X - BInvrse) - (One - BInvrse);
BadCond(Defect, "Calculated");
printf(" %.17e for\n", POW(X,Z));
printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
printf("\tdiffers from correct value by %.17e .\n", Q);
printf("\tThis much error may spoil financial\n");
printf("\tcalculations involving tiny interest rates.\n");
break;
}
else {
Z = (Y - X) * Two + Y;
X = Y;
Y = Z;
Z = One + (X - F9)*(X - F9);
if (Z > One && I < NoTrials) I++;
else {
if (X > One) {
if (N == 0)
printf("Accuracy seems adequate.\n");
break;
}
else {
X = One + U2;
Y = U2 + U2;
Y += X;
I = 1;
}
}
}
}
/*=============================================*/
Milestone = 150;
/*=============================================*/
printf("Testing powers Z^Q at four nearly extreme values.\n");
N = 0;
Z = A1;
Q = FLOOR(Half - LOG(C) / LOG(A1));
Break = False;
do {
X = CInvrse;
Y = POW(Z, Q);
IsYeqX();
Q = - Q;
X = C;
Y = POW(Z, Q);
IsYeqX();
if (Z < One) Break = True;
else Z = AInvrse;
} while ( ! (Break));
PrintIfNPositive();
if (N == 0) printf(" ... no discrepancies found.\n");
printf("\n");
/*=============================================*/
Milestone = 160;
/*=============================================*/
Pause();
printf("Searching for Overflow threshold:\n");
printf("This may generate an error.\n");
Y = - CInvrse;
V9 = HInvrse * Y;
sigsave = sigfpe;
if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
do {
V = Y;
Y = V9;
V9 = HInvrse * Y;
} while(V9 < Y);
I = 1;
overflow:
sigsave = 0;
Z = V9;
printf("Can `Z = -Y' overflow?\n");
printf("Trying it on Y = %.17e .\n", Y);
V9 = - Y;
V0 = V9;
if (V - Y == V + V0) printf("Seems O.K.\n");
else {
printf("finds a ");
BadCond(Flaw, "-(-Y) differs from Y.\n");
}
if (Z != Y) {
BadCond(Serious, "");
printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
}
if (I) {
Y = V * (HInvrse * U2 - HInvrse);
Z = Y + ((One - HInvrse) * U2) * V;
if (Z < V0) Y = Z;
if (Y < V0) V = Y;
if (V0 - V < V0) V = V0;
}
else {
V = Y * (HInvrse * U2 - HInvrse);
V = V + ((One - HInvrse) * U2) * Y;
}
printf("Overflow threshold is V = %.17e .\n", V);
if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
else printf("There is no saturation value because \
the system traps on overflow.\n");
V9 = V * One;
printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
V9 = V / One;
printf(" nor for V / 1 = %.17e .\n", V9);
printf("Any overflow signal separating this * from the one\n");
printf("above is a DEFECT.\n");
/*=============================================*/
Milestone = 170;
/*=============================================*/
if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
BadCond(Failure, "Comparisons involving ");
printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
V, V0, UfThold);
}
/*=============================================*/
Milestone = 175;
/*=============================================*/
printf("\n");
for(Indx = 1; Indx <= 3; ++Indx) {
switch (Indx) {
case 1: Z = UfThold; break;
case 2: Z = E0; break;
case 3: Z = PseudoZero; break;
}
if (Z != Zero) {
V9 = SQRT(Z);
Y = V9 * V9;
if (Y / (One - Radix * E9) < Z
|| Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
if (V9 > U1) BadCond(Serious, "");
else BadCond(Defect, "");
printf("Comparison alleges that what prints as Z = %.17e\n", Z);
printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
}
}
}
/*=============================================*/
Milestone = 180;
/*=============================================*/
for(Indx = 1; Indx <= 2; ++Indx) {
if (Indx == 1) Z = V;
else Z = V0;
V9 = SQRT(Z);
X = (One - Radix * E9) * V9;
V9 = V9 * X;
if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
Y = V9;
if (X < W) BadCond(Serious, "");
else BadCond(Defect, "");
printf("Comparison alleges that Z = %17e\n", Z);
printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
}
}
/*=============================================*/
/*SPLIT
}
#include "paranoia.h"
int part8(VOID){
*/
Milestone = 190;
/*=============================================*/
Pause();
X = UfThold * V;
Y = Radix * Radix;
if (X*Y < One || X > Y) {
if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
else BadCond(Flaw, "");
printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
X, "is too far from 1.\n");
}
/*=============================================*/
Milestone = 200;
/*=============================================*/
for (Indx = 1; Indx <= 5; ++Indx) {
X = F9;
switch (Indx) {
case 2: X = One + U2; break;
case 3: X = V; break;
case 4: X = UfThold; break;
case 5: X = Radix;
}
Y = X;
sigsave = sigfpe;
if (setjmp(ovfl_buf))
printf(" X / X traps when X = %g\n", X);
else {
V9 = (Y / X - Half) - Half;
if (V9 == Zero) continue;
if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
else BadCond(Serious, "");
printf(" X / X differs from 1 when X = %.17e\n", X);
printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
}
sigsave = 0;
}
/*=============================================*/
Milestone = 210;
/*=============================================*/
MyZero = Zero;
printf("\n");
printf("What message and/or values does Division by Zero produce?\n") ;
#ifndef NOPAUSE
printf("This can interupt your program. You can ");
printf("skip this part if you wish.\n");
printf("Do you wish to compute 1 / 0? ");
fflush(stdout);
read (KEYBOARD, ch, 8);
if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif
sigsave = sigfpe;
printf(" Trying to compute 1 / 0 produces ...");
if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero);
sigsave = 0;
#ifndef NOPAUSE
}
else printf("O.K.\n");
printf("\nDo you wish to compute 0 / 0? ");
fflush(stdout);
read (KEYBOARD, ch, 80);
if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif
sigsave = sigfpe;
printf("\n Trying to compute 0 / 0 produces ...");
if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero);
sigsave = 0;
#ifndef NOPAUSE
}
else printf("O.K.\n");
#endif
/*=============================================*/
Milestone = 220;
/*=============================================*/
Pause();
printf("\n");
{
static char *msg[] = {
"FAILUREs encountered =",
"SERIOUS DEFECTs discovered =",
"DEFECTs discovered =",
"FLAWs discovered =" };
int i;
for(i = 0; i < 4; i++) if (ErrCnt[i])
printf("The number of %-29s %d.\n",
msg[i], ErrCnt[i]);
}
printf("\n");
if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
+ ErrCnt[Flaw]) > 0) {
if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
Defect] == 0) && (ErrCnt[Flaw] > 0)) {
printf("The arithmetic diagnosed seems ");
printf("Satisfactory though flawed.\n");
}
if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
&& ( ErrCnt[Defect] > 0)) {
printf("The arithmetic diagnosed may be Acceptable\n");
printf("despite inconvenient Defects.\n");
}
if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
printf("The arithmetic diagnosed has ");
printf("unacceptable Serious Defects.\n");
}
if (ErrCnt[Failure] > 0) {
printf("Potentially fatal FAILURE may have spoiled this");
printf(" program's subsequent diagnoses.\n");
}
}
else {
printf("No failures, defects nor flaws have been discovered.\n");
if (! ((RMult == Rounded) && (RDiv == Rounded)
&& (RAddSub == Rounded) && (RSqrt == Rounded)))
printf("The arithmetic diagnosed seems Satisfactory.\n");
else {
if (StickyBit >= One &&
(Radix - Two) * (Radix - Nine - One) == Zero) {
printf("Rounding appears to conform to ");
printf("the proposed IEEE standard P");
if ((Radix == Two) &&
((Precision - Four * Three * Two) *
( Precision - TwentySeven -
TwentySeven + One) == Zero))
printf("754");
else printf("854");
if (IEEE) printf(".\n");
else {
printf(",\nexcept for possibly Double Rounding");
printf(" during Gradual Underflow.\n");
}
}
printf("The arithmetic diagnosed appears to be Excellent!\n");
}
}
if (fpecount)
printf("\nA total of %d floating point exceptions were registered.\n",
fpecount);
printf("END OF TEST.\n");
return 0;
}
/*SPLIT subs.c
#include "paranoia.h"
*/
FLOAT
Sign (FP X)
#ifdef KR_headers
FLOAT X;
#endif
{ return X >= 0. ? 1.0 : -1.0; }
void
Pause(VOID)
{
#ifndef NOPAUSE
char ch[8];
printf("\nTo continue, press RETURN");
fflush(stdout);
read(KEYBOARD, ch, 8);
#endif
printf("\nDiagnosis resumes after milestone Number %d", Milestone);
printf(" Page: %d\n\n", PageNo);
++Milestone;
++PageNo;
}
void
TstCond (INT K, INT Valid, CHARP T)
#ifdef KR_headers
int K, Valid;
char *T;
#endif
{ if (! Valid) { BadCond(K,T); printf(".\n"); } }
void
BadCond(INT K, CHARP T)
#ifdef KR_headers
int K;
char *T;
#endif
{
static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
ErrCnt [K] = ErrCnt [K] + 1;
printf("%s: %s", msg[K], T);
}
FLOAT
Random(VOID)
/* Random computes
X = (Random1 + Random9)^5
Random1 = X - FLOOR(X) + 0.000005 * X;
and returns the new value of Random1
*/
{
FLOAT X, Y;
X = Random1 + Random9;
Y = X * X;
Y = Y * Y;
X = X * Y;
Y = X - FLOOR(X);
Random1 = Y + X * 0.000005;
return(Random1);
}
void
SqXMinX (INT ErrKind)
#ifdef KR_headers
int ErrKind;
#endif
{
FLOAT XA, XB;
XB = X * BInvrse;
XA = X - XB;
SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
if (SqEr != Zero) {
if (SqEr < MinSqEr) MinSqEr = SqEr;
if (SqEr > MaxSqEr) MaxSqEr = SqEr;
J = J + 1.0;
BadCond(ErrKind, "\n");
printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr);
printf("\tinstead of correct value 0 .\n");
}
}
void
NewD(VOID)
{
X = Z1 * Q;
X = FLOOR(Half - X / Radix) * Radix + X;
Q = (Q - X * Z) / Radix + X * X * (D / Radix);
Z = Z - Two * X * D;
if (Z <= Zero) {
Z = - Z;
Z1 = - Z1;
}
D = Radix * D;
}
void
SR3750(VOID)
{
if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
I = I + 1;
X2 = SQRT(X * D);
Y2 = (X2 - Z2) - (Y - Z2);
X2 = X8 / (Y - Half);
X2 = X2 - Half * X2 * X2;
SqEr = (Y2 + Half) + (Half - X2);
if (SqEr < MinSqEr) MinSqEr = SqEr;
SqEr = Y2 - X2;
if (SqEr > MaxSqEr) MaxSqEr = SqEr;
}
}
void
IsYeqX(VOID)
{
if (Y != X) {
if (N <= 0) {
if (Z == Zero && Q <= Zero)
printf("WARNING: computing\n");
else BadCond(Defect, "computing\n");
printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
printf("\tyielded %.17e;\n", Y);
printf("\twhich compared unequal to correct %.17e ;\n",
X);
printf("\t\tthey differ by %.17e .\n", Y - X);
}
N = N + 1; /* ... count discrepancies. */
}
}
void
SR3980(VOID)
{
do {
Q = (FLOAT) I;
Y = POW(Z, Q);
IsYeqX();
if (++I > M) break;
X = Z * X;
} while ( X < W );
}
void
PrintIfNPositive(VOID)
{
if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
}
void
TstPtUf(VOID)
{
N = 0;
if (Z != Zero) {
printf("Since comparison denies Z = 0, evaluating ");
printf("(Z + Z) / Z should be safe.\n");
sigsave = sigfpe;
if (setjmp(ovfl_buf)) goto very_serious;
Q9 = (Z + Z) / Z;
printf("What the machine gets for (Z + Z) / Z is %.17e .\n",
Q9);
if (FABS(Q9 - Two) < Radix * U2) {
printf("This is O.K., provided Over/Underflow");
printf(" has NOT just been signaled.\n");
}
else {
if ((Q9 < One) || (Q9 > Two)) {
very_serious:
N = 1;
ErrCnt [Serious] = ErrCnt [Serious] + 1;
printf("This is a VERY SERIOUS DEFECT!\n");
}
else {
N = 1;
ErrCnt [Defect] = ErrCnt [Defect] + 1;
printf("This is a DEFECT!\n");
}
}
sigsave = 0;
V9 = Z * One;
Random1 = V9;
V9 = One * Z;
Random2 = V9;
V9 = Z / One;
if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
if (N > 0) Pause();
}
else {
N = 1;
BadCond(Defect, "What prints as Z = ");
printf("%.17e\n\tcompares different from ", Z);
if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
if (! ((Z == Random2)
|| (Random2 == Random1)))
printf("1 * Z == %g\n", Random2);
if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
if (Random2 != Random1) {
ErrCnt [Defect] = ErrCnt [Defect] + 1;
BadCond(Defect, "Multiplication does not commute!\n");
printf("\tComparison alleges that 1 * Z = %.17e\n",
Random2);
printf("\tdiffers from Z * 1 = %.17e\n", Random1);
}
Pause();
}
}
}
void
notify(CHARP s)
#ifdef KR_headers
char *s;
#endif
{
printf("%s test appears to be inconsistent...\n", s);
printf(" PLEASE NOTIFY KARPINKSI!\n");
}
/*SPLIT msgs.c
#include "paranoia.h"
*/
void
msglist(CHARPP s)
#ifdef KR_headers
char **s;
#endif
{ while(*s) printf("%s\n", *s++); }
void
Instructions(VOID)
{
static char *instr[] = {
"Lest this program stop prematurely, i.e. before displaying\n",
" `END OF TEST',\n",
"try to persuade the computer NOT to terminate execution when an",
"error like Over/Underflow or Division by Zero occurs, but rather",
"to persevere with a surrogate value after, perhaps, displaying some",
"warning. If persuasion avails naught, don't despair but run this",
"program anyway to see how many milestones it passes, and then",
"amend it to make further progress.\n",
"Answer questions with Y, y, N or n (unless otherwise indicated).\n",
0};
msglist(instr);
}
void
Heading(VOID)
{
static char *head[] = {
"Users are invited to help debug and augment this program so it will",
"cope with unanticipated and newly uncovered arithmetic pathologies.\n",
"Please send suggestions and interesting results to",
"\tRichard Karpinski",
"\tComputer Center U-76",
"\tUniversity of California",
"\tSan Francisco, CA 94143-0704, USA\n",
"In doing so, please include the following information:",
#ifdef Single
"\tPrecision:\tsingle;",
#else
"\tPrecision:\tdouble;",
#endif
"\tVersion:\t10 February 1989;",
"\tComputer:\n",
"\tCompiler:\n",
"\tOptimization level:\n",
"\tOther relevant compiler options:",
0};
msglist(head);
}
void
Characteristics(VOID)
{
static char *chars[] = {
"Running this program should reveal these characteristics:",
" Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
" Precision = number of significant digits carried.",
" U2 = Radix/Radix^Precision = One Ulp",
"\t(OneUlpnit in the Last Place) of 1.000xxx .",
" U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
" Adequacy of guard digits for Mult., Div. and Subt.",
" Whether arithmetic is chopped, correctly rounded, or something else",
"\tfor Mult., Div., Add/Subt. and Sqrt.",
" Whether a Sticky Bit used correctly for rounding.",
" UnderflowThreshold = an underflow threshold.",
" E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
" V = an overflow threshold, roughly.",
" V0 tells, roughly, whether Infinity is represented.",
" Comparisions are checked for consistency with subtraction",
"\tand for contamination with pseudo-zeros.",
" Sqrt is tested. Y^X is not tested.",
" Extra-precise subexpressions are revealed but NOT YET tested.",
" Decimal-Binary conversion is NOT YET tested for accuracy.",
0};
msglist(chars);
}
void
History(VOID)
{ /* History */
/* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
with further massaging by David M. Gay. */
static char *hist[] = {
"The program attempts to discriminate among",
" FLAWs, like lack of a sticky bit,",
" Serious DEFECTs, like lack of a guard digit, and",
" FAILUREs, like 2+2 == 5 .",
"Failures may confound subsequent diagnoses.\n",
"The diagnostic capabilities of this program go beyond an earlier",
"program called `MACHAR', which can be found at the end of the",
"book `Software Manual for the Elementary Functions' (1980) by",
"W. J. Cody and W. Waite. Although both programs try to discover",
"the Radix, Precision and range (over/underflow thresholds)",
"of the arithmetic, this program tries to cope with a wider variety",
"of pathologies, and to say how well the arithmetic is implemented.",
"\nThe program is based upon a conventional radix representation for",
"floating-point numbers, but also allows logarithmic encoding",
"as used by certain early WANG machines.\n",
"BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
"see source comments for more history.",
0};
msglist(hist);
}
#--------------------------------------------------------------------------*/
# Copyright (c) Qiagen Instruments AG */
# Title : uClinux on Coldfire test program */
# Project Name : */
# Project : FIW-00000039-000-A */
# Subproject : WEB Maintenance Interface */
# File : Makefile */
#--------------------------------------------------------------------------*/
# Created by Felix Daners Engineering, 8200 Schaffhausen, Switzerland */
#--------------------------------------------------------------------------*/
#
# $Author: fd $
# $Revision: 1.1 $
# $Date: 2005/12/18 17:12:42 $
# $Name: $
# $RCSfile: Makefile,v $
# $Locker: $
#
#--------------------------------------------------------------------------*/
# Comments: */
override CROSS_COMPILE=m68k-uclinux-
export CROSS_COMPILE
PLATFORMFLAGS = -D__uClinux__ -D_CONSOLE\
-D_REENTRANT -D_CAN_ -D_FB_COP_ -fsigned-char\
$(TESTFLAGS) $(CONFIG_FLAGS)
LDFLAGS = -Wl,-elf2flt
FLTFLAGS = '-s 16000'
CPUFLAGS=-m528x
EXE1=test-float
EXE2=test-double
EXE3=test-rint
# Set UPLOADDIR to your (witeable) directory for uploads:
CFLAGS = -g $(PLATFORMFLAGS) $(CPUFLAGS) -D_GNU_SOURCE
OBJS =
LDLIBS += -lc -lm
override CC=$(CROSS_COMPILE)gcc
override AR=$(CROSS_COMPILE)ar
# targets
all: $(EXE1) $(EXE2) $(EXE3)
$(EXE1): paranoia.c $(OBJS)
$(CC) $(CFLAGS) $(LDFLAGS) -DSingle -o $@ $< $(OBJS) -Wl,-znodefaultlib $(LDLIBS)
$(EXE2): paranoia.c $(OBJS)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ $< $(OBJS) $(LDLIBS)
$(EXE3): rintf.c $(OBJS)
$(CC) $(CFLAGS) $(LDFLAGS) -o $@ $< $(OBJS) $(LDLIBS)
clean:
-rm -f $(EXE1)
-rm -f $(EXE2)
-rm -f $(EXE3)