#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__ */