/*
 *  Calcule signe(a*d-b*c), de manière exacte, avec des doubles.
 *  © Sylvain Pion, le 19/9/96.
 *
 *  Méthode:
 *    - filtre
 *    - au besoin, on utilise ensuite l'arithmétique modulaire.
 *
 *  Conditions d'utilisation:
 *    - architecture 32 bits supportant le calcul flottant IEEE 754.
 *    - les entrées a,b,c,d doivent être des entiers, de valeur absolue
 *      inférieure à 2^53, codés par des doubles.
 */

/*
 * Quelques constantes utiles:
 * CST_2_24 == 2^24
 * CST_2_32 == 2^32
 * CST_3_75 == 3 * 2^75  (75 = 51 + 24)
 * CST_3_83 == 3 * 2^83  (83 = 51 + 32)
 */

#define CST_2_24 (65536.0*256.0)
#define CST_2_32 (CST_2_24*256.0)
#define CST_3_51 (3.0*8.0*CST_2_24*CST_2_24)
#define CST_3_75 (3.0*8.0*CST_2_24*CST_2_24*CST_2_24)
#define CST_3_83 (256.0*CST_3_75)

/*
 * Le choix de la décomposition:
 *  * 2^32 et (2^32-1), c'est bien, sauf que le calcul en flottant
 *    pour 2^32-1 ne rentre pas...
 *  * 2^32 et (2^24-1), ça rentre, mais la décomposition est moins simple.
 *  * 2^27 et (2^27-1), la décomposition est facile, ça rentre,
 *    mais le calul est moins simple (celui modulo 2^27 vs 2^32)
 *    => faire les 2 en flottants => est-ce que ça apporte quelque
 *       chose par rapport à une décomposition normale ???
 *
 * Je choisis: 2^32 et 2^24-1.
 */

#define sgn(a) ( ((a)>0) ? 1 : ( ((a)<0) ? -1 : 0) )

#define dec_32(e) ((signed int) ((e) - (((e) + CST_3_83) - CST_3_83)))

#define dec_24_1(e) ((e) - (CST_2_24-1) * \
                           (((e)/(CST_2_24-1)+CST_3_51)-CST_3_51))


int modulo_signe(double a, double b, double c, double d)
{
	double g = a*d - b*c;

	if (g != 0)  {return( (g>0) ? 1 : -1 );}
	else {

	/* Le déterminant modulo 2^24-1. */
	/* Le calcul est effectué avec des doubles. */

		g = dec_24_1(a)*dec_24_1(d) - dec_24_1(b)*dec_24_1(c);
		g = dec_24_1(g);

        /* Le déterminant modulo 2^32. */
	/* Le calcul est effectué avec des signed int. */
 
		a = (double) (dec_32(a)*dec_32(d) - dec_32(b)*dec_32(c));

	/* On recompose le déterminant (th. des restes chinois). */

		b = dec_24_1 (a);
		g = (g - b) * 65536.0;
		b = (dec_24_1 (g)) * CST_2_32 + a;

		return sgn(b);
	};
}