home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Geek Gadgets 1
/
ADE-1.bin
/
ade-dist
/
ixemul-45.0-src.tgz
/
tar.out
/
contrib
/
ixemul
/
general
/
arith.c
< prev
next >
Wrap
C/C++ Source or Header
|
1996-09-28
|
13KB
|
777 lines
#define ALL
#include "common.h"
#include "defs.h"
#include <proto/mathieeedoubbas.h>
#include <proto/mathieeedoubtrans.h>
#include <proto/mathieeesingbas.h>
/* Special patch to work around bug in IEEEDPCmp:
* if the first 32 bits of both doubles are equal, and
* both doubles are negative, then the result can no longer
* be trusted.
*
* This is the output of a small test program:
*
* test -2.000001 -2.0000009
* a = -2.0000009999999996956888, b = -2.0000008999999998593466
* (0xc0000000 0x8637bd05, 0xc0000000 0x78cbc3b8)
* cmp(a,b) = 0, cmp(b,a) = 1
*
* test -2.0000001 -2.0000002
* a = -2.00000009999999983634211, b = -2.0000001999999996726842
* (0xc0000000 0xd6bf94d, 0xc0000000 0x1ad7f29a)
* cmp(a,b) = 1, cmp(b,a) = 1
*
* As you can see, the results are wrong.
*
* So, we just make both variables positive and exchange them before
* passing them to IEEEDPCmp.
*
* This bug was discovered by Bart Van Assche, thanks!
*/
static int ieeedpcmp(double a, double b)
{
if (*((char *)&a) & *((char *)&b) & 0x80) /* both doubles negative? */
{
*((char *)&a) &= 0x7f; /* yes, make positive */
*((char *)&b) &= 0x7f;
return IEEEDPCmp(b, a); /* pass them to IEEEDPCmp the other way round */
}
return IEEEDPCmp(a, b);
}
#if defined(L__eqdf2) || defined(ALL)
int __eqdf2 (double a, double b)
{
return ieeedpcmp (a, b);
}
#endif
#if defined(L__eqsf2) || defined(ALL)
int __eqsf2 (FLOAT a, FLOAT b)
{
return IEEESPCmp (a, b);
}
#endif
#if defined(L__fixsfsi) || defined(ALL)
SItype __fixsfsi (FLOAT a)
{
return IEEESPFix(a);
}
#endif
#if defined(L__floatsisf) || defined(ALL)
SFVALUE __floatsisf (SItype a)
{
return IEEESPFlt(a);
}
#endif
#if defined(L__gedf2) || defined(ALL)
int __gedf2 (double a, double b)
{
return ieeedpcmp (a, b);
}
#endif
#if defined(L__gesf2) || defined(ALL)
int __gesf2 (FLOAT a, FLOAT b)
{
return IEEESPCmp (a, b);
}
#endif
#if defined(L__gtdf2) || defined(ALL)
int __gtdf2 (double a, double b)
{
return ieeedpcmp (a, b);
}
#endif
#if defined(L__gtsf2) || defined(ALL)
int __gtsf2 (FLOAT a, FLOAT b)
{
return IEEESPCmp (a, b);
}
#endif
#if defined(L__ledf2) || defined(ALL)
int __ledf2 (double a, double b)
{
return ieeedpcmp (a, b);
}
#endif
#if defined(L__lesf2) || defined(ALL)
int __lesf2 (FLOAT a, FLOAT b)
{
return IEEESPCmp (a, b);
}
#endif
#if defined(L__ltdf2) || defined(ALL)
int __ltdf2 (double a, double b)
{
return ieeedpcmp (a, b);
}
#endif
#if defined(L__ltsf2) || defined(ALL)
int __ltsf2 (FLOAT a, FLOAT b)
{
return IEEESPCmp (a, b);
}
#endif
#if defined(L__nedf2) || defined(ALL)
int __nedf2 (double a, double b)
{
return !!ieeedpcmp (a, b);
}
#endif
#if defined(L__nesf2) || defined(ALL)
int __nesf2 (FLOAT a, FLOAT b)
{
return !!IEEESPCmp (a, b);
}
#endif
#if defined (mc68020) || defined (mc68030) || defined (mc68040) || defined (mc68060)
/* int / int */
#if defined(L__divsi3) || defined(ALL)
ENTRY(__divsi3)
asm("
movel sp@(4),d0
divsl sp@(8),d0
rts
");
#endif
/* int % int */
#if defined(L__modsi3) || defined(ALL)
ENTRY(__modsi3)
asm("
movel sp@(4),d1
divsll sp@(8),d0:d1
rts
");
#endif
/* int * int */
#if defined(L__mulsi3) || defined(ALL)
ENTRY(__mulsi3)
asm("
movel sp@(4),d0
mulsl sp@(8),d0
rts
");
#endif
/* unsigned / unsigned */
#if defined(L__udivsi3) || defined(ALL)
ENTRY(__udivsi3)
asm("
movel sp@(4),d0
divul sp@(8),d0
rts
");
#endif
/* unsigned % unsigned */
#if defined(L__umodsi3) || defined(ALL)
ENTRY(__umodsi3)
asm("
movel sp@(4),d1
divull sp@(8),d0:d1
rts
");
#endif
/* unsigned * unsigned */
#if defined(L__umulsi3) || defined(ALL)
ENTRY(__umulsi3)
asm("
movel sp@(4),d0
mulul sp@(8),d0
rts
");
#endif
#else
#if defined(L__divsi3) || defined(ALL)
SItype __divsi3 (SItype a, SItype b)
{
unsigned SItype q, r;
int neg = (a < 0) != (b < 0);
if (a < 0) a = -a;
if (b < 0) b = -b;
divmodu (q, r, a, b);
return neg ? -q : q;
}
#endif
#if defined(L__modsi3) || defined(ALL)
SItype __modsi3 (SItype a, SItype b)
{
unsigned SItype q, r;
int neg = (a < 0);
if (a < 0) a = -a;
if (b < 0) b = -b;
divmodu (q, r, a, b);
return neg ? -r : r;
}
#endif
#if defined(L__mulsi3) || defined(ALL)
SItype __mulsi3 (SItype a, SItype b)
{
int neg = (a < 0) != (b < 0);
SItype res;
if (a < 0) a = -a;
if (b < 0) b = -b;
res = mulu (a,b);
return neg ? -res : res;
}
#endif
#if defined(L__udivsi3) || defined(ALL)
unsigned SItype __udivsi3 (unsigned SItype a, unsigned SItype b)
{
unsigned SItype q, r;
divmodu (q, r, a, b);
return q;
}
#endif
#if defined(L__umodsi3) || defined(ALL)
unsigned SItype __umodsi3 (unsigned SItype a, unsigned SItype b)
{
unsigned SItype q, r;
divmodu (q, r, a, b);
return r;
}
#endif
#endif
/* -double */
#if defined(L__negdf2) || defined(ALL)
ENTRY(__negdf2)
asm("
movel sp@(4),d0
movel sp@(8),d1
bchg #31,d0
rts
");
#endif
/* -single */
#if defined(L__negsf2) || defined(ALL)
ENTRY(__negsf2)
asm("
movel sp@(4),d0
bchg #31,d0
rts
");
#endif
#ifdef __HAVE_68881__
/* double + double */
#if defined(L__adddf3) || defined(ALL)
ENTRY(__adddf3)
asm("
fmoved sp@(4),fp0
faddd sp@(12),fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/* single + single */
#if defined(L__addsf3) || defined(ALL)
ENTRY(__addsf3)
asm("
fmoves sp@(4),fp0
fadds sp@(8),fp0
fmoves fp0,d0
rts
");
#endif
/* double > double: 1 */
/* double < double: -1 */
/* double == double: 0 */
#if defined(L__cmpdf2) || defined(ALL)
ENTRY(__cmpdf2)
asm("
fmoved sp@(4),fp0
fcmpd sp@(12),fp0
fbgt Lagtb1
fslt d0
extbl d0
rts
Lagtb1:
moveq #1,d0
rts
");
#endif
/* single > single: 1 */
/* single < single: -1 */
/* single == single: 0 */
#if defined(L__cmpsf2) || defined(ALL)
ENTRY(__cmpsf2)
asm("
fmoves sp@(4),fp0
fcmps sp@(8),fp0
fbgt Lagtb2
fslt d0
extbl d0
rts
Lagtb2:
moveq #1,d0
rts
");
#endif
/* double / double */
#if defined(L__divdf3) || defined(ALL)
ENTRY(__divdf3)
asm("
fmoved sp@(4),fp0
fdivd sp@(12),fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/* single / single */
#if defined(L__divsf3) || defined(ALL)
ENTRY(__divsf3)
asm("
fmoves sp@(4),fp0
fdivs sp@(8),fp0
fmoves fp0,d0
rts
");
#endif
/* (double) float */
#if defined(L__extendsfdf2) || defined(ALL)
ENTRY(__extendsfdf2)
asm("
fmoves sp@(4),fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
#if defined(Lfabs) || defined(ALL)
ENTRY(fabs)
asm("
fmoved sp@(4),fp0
fjnlt L1
fnegx fp0
L1:
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/* (int) double */
#if defined(L__fixdfsi) || defined(ALL)
ENTRY(__fixdfsi)
asm("
fintrzd sp@(4),fp0
fmovel fp0,d0
rts
");
#endif
/* (double) int */
#if defined(L__floatsidf) || defined(ALL)
ENTRY(__floatsidf)
asm("
fmovel sp@(4),fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/*
* double frexp(val, eptr)
* returns: x s.t. val = x * (2 ** n), with n stored in *eptr
*/
#if defined(Lfrexp) || defined(ALL)
ENTRY(frexp)
asm("
fmoved sp@(4),fp1
fgetmanx fp1,fp0
fgetexpx fp1
fmovel fp1,d0
addql #1,d0
movel sp@(12),a0
movel d0,a0@
fdivl #2,fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/*
* double ldexp(val, exp)
* returns: val * (2**exp), for integer exp
*/
#if defined(Lldexp) || defined(ALL)
ENTRY(ldexp)
asm("
fmoved sp@(4),fp0
fbeq Ldone
ftwotoxl sp@(12),fp1
fmulx fp1,fp0
Ldone:
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/*
* double modf(val, iptr)
* returns: xxx and n (in *iptr) where val == n.xxx
*/
#if defined(Lmodf) || defined(ALL)
ENTRY(modf)
asm("
fmoved sp@(4),fp0
movel sp@(12),a0
fintrzx fp0,fp1
fmoved fp1,a0@
fsubx fp1,fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/* double * double */
#if defined(L__muldf3) || defined(ALL)
ENTRY(__muldf3)
asm("
fmoved sp@(4),fp0
fmuld sp@(12),fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/* single * single */
#if defined(L__mulsf3) || defined(ALL)
ENTRY(__mulsf3)
asm("
fmoves sp@(4),fp0
fmuls sp@(8),fp0
fmoves fp0,d0
rts
");
#endif
/* double - double */
#if defined(L__subdf3) || defined(ALL)
ENTRY(__subdf3)
asm("
fmoved sp@(4),fp0
fsubd sp@(12),fp0
fmoved fp0,sp@-
movel sp@+,d0
movel sp@+,d1
rts
");
#endif
/* single - single */
#if defined(L__subsf3) || defined(ALL)
ENTRY(__subsf3)
asm("
fmoves sp@(4),fp0
fsubs sp@(8),fp0
fmoves fp0,d0
rts
");
#endif
/* (float) double */
#if defined(L__truncdfsf2) || defined(ALL)
ENTRY(__truncdfsf2)
asm("
fmoved sp@(4),fp0
fmoves fp0,d0
rts
");
#endif
#else /* __HAVE_68881__ */
#if defined(L__adddf3) || defined(ALL)
double __adddf3 (double a, double b)
{
return IEEEDPAdd (a, b);
}
#endif
#if defined(L__addsf3) || defined(ALL)
SFVALUE __addsf3 (FLOAT a, FLOAT b)
{
return IEEESPAdd (a, b);
}
#endif
#if defined(L__cmpdf2) || defined(ALL)
int __cmpdf2 (double a, double b)
{
return ieeedpcmp (a, b);
}
#endif
#if defined(L__cmpsf2) || defined(ALL)
int __cmpsf2 (FLOAT a, FLOAT b)
{
return IEEESPCmp (a, b);
}
#endif
#if defined(L__divdf3) || defined(ALL)
double __divdf3 (double a, double b)
{
return IEEEDPDiv (a, b);
}
#endif
#if defined(L__divsf3) || defined(ALL)
SFVALUE __divsf3 (FLOAT a, FLOAT b)
{
return IEEESPDiv (a, b);
}
#endif
#if defined(L__extendsfdf2) || defined(ALL)
double __extendsfdf2 (FLOAT a)
{
return IEEEDPFieee(a);
}
#endif
#if defined(Lfabs) || defined(ALL)
double fabs (double d)
{
return d < 0.0 ? -d : d;
}
#endif
#if defined(L__fixdfsi) || defined(ALL)
SItype __fixdfsi (double a)
{
return IEEEDPFix (a);
}
#endif
#if defined(L__floatsidf) || defined(ALL)
double __floatsidf (SItype a)
{
return IEEEDPFlt (a);
}
#endif
/*
* the call
* x = frexp(arg,&exp);
* must return a double fp quantity x which is <1.0
* and the corresponding binary exponent "exp".
* such that
* arg = x*2^exp
* if the argument is 0.0, return 0.0 mantissa and 0 exponent.
*/
#if defined(Lfrexp) || defined(ALL)
double frexp (double x, int *i)
{
int neg = 0;
int j = 0;
if (x < 0)
{
x = -x;
neg = 1;
}
if (x >= 1.0)
while (x >= 1.0)
{
j = j + 1;
x = x / 2;
}
else if (x < 0.5 && x != 0.0)
while (x < 0.5)
{
j = j - 1;
x = 2 * x;
}
*i = j;
if (neg)
x = -x;
return x;
}
#endif
#if defined(Lldexp) || defined(ALL)
/*
* ldexp returns the quanity "value" * 2 ^ "exp"
*
* For the mc68000 using IEEE format the double precision word format is:
*
* WORD N => SEEEEEEEEEEEMMMM
* WORD N+1 => MMMMMMMMMMMMMMMM
* WORD N+2 => MMMMMMMMMMMMMMMM
* WORD N+3 => MMMMMMMMMMMMMMMM
*
* Where: S => Sign bit
* E => Exponent
* X => Ignored (set to 0)
* M => Mantissa bit
*
* NOTE: Beware of 0.0; on some machines which use excess 128 notation for the
* exponent, if the mantissa is zero the exponent is also.
*
*/
#define MANT_MASK 0x800FFFFF /* Mantissa extraction mask */
#define ZPOS_MASK 0x3FF00000 /* Positive # mask for exp = 0 */
#define ZNEG_MASK 0x3FF00000 /* Negative # mask for exp = 0 */
#define EXP_MASK 0x7FF00000 /* Mask for exponent */
#define EXP_SHIFTS 20 /* Shifts to get into LSB's */
#define EXP_BIAS 1023 /* Exponent bias */
union dtol
{
double dval;
int ival[2];
};
double ldexp (double value, int exp)
{
union dtol number;
int *iptr, cexp;
if (value == 0.0)
return (0.0);
number.dval = value;
iptr = &number.ival[0];
cexp = (((*iptr) & EXP_MASK) >> EXP_SHIFTS) - EXP_BIAS;
*iptr &= ~EXP_MASK;
exp += EXP_BIAS;
*iptr |= ((exp + cexp) << EXP_SHIFTS) & EXP_MASK;
return (number.dval);
}
#endif
/*
* modf(value, iptr): return fractional part of value, and stores the
* integral part into iptr (a pointer to double).
*/
#if defined(Lmodf) || defined(ALL)
double modf (double value, double *iptr)
{
/* if value negative */
if (IEEEDPTst (value) < 0)
{
/* in that case, the integer part is calculated by ceil() */
*iptr = IEEEDPCeil (value);
return IEEEDPSub (*iptr, value);
}
else
{
/* if positive, we go for the floor() */
*iptr = IEEEDPFloor (value);
return IEEEDPSub (value, *iptr);
}
}
#endif
#if defined(L__muldf3) || defined(ALL)
double __muldf3 (double a, double b)
{
return IEEEDPMul (a, b);
}
#endif
#if defined(L__mulsf3) || defined(ALL)
SFVALUE __mulsf3 (FLOAT a, FLOAT b)
{
return IEEESPMul (a, b);
}
#endif
#if defined(L__subdf3) || defined(ALL)
double __subdf3 (double a, double b)
{
return IEEEDPSub (a, b);
}
#endif
#if defined(L__subsf3) || defined(ALL)
SFVALUE __subsf3 (FLOAT a, FLOAT b)
{
return IEEESPSub (a, b);
}
#endif
#if defined(L__truncdfsf2) || defined(ALL)
SFVALUE __truncdfsf2 (double a)
{
return IEEEDPTieee(a);
}
#endif
#endif /* __HAVE_68881__ */