🚀 go-pugleaf

/*
 * Unsigned divide 64-bit / 32-bit
 *
 * Input:
 *     uint D    32-bit divisor
 *     uint H    upper 32-bits of dividend
 *     uint L    lower 32-bits of dividend
 *     uint *RP   pointer to 32 bit remainder
 *     uint *QP   pointer to 32 bit quotient
 *
 * Output:
 *     int       overflow flag
 *     *RP       32-bit remainder
 *     *QP       32-bit quotient
 *
 * Notes:
 *     1) If quotient does not fit in 32-bits then 1 is returned
 *
 *     2) The D2INT macro truncates a double to an integer.
 *        It is coded for BIG ENDIAN.
 *        The rounding mode must be set to round towards zero.
 *
 * Authors: Mauricio Breternitz Jr.     (Motorola)
 *          Arturo Martin-de-Nicolas    (IBM)
 *          Joseph M. Orost             (AT&T Bell Labs)
 *
 * Algorithm:
 *      1) Verify that the quotient will fit in 32-bits.
 *      2) Compute the reciprocal of the divisor using
 *         the floating point unit.
 *      3) Compute an approximation of the quotient by
 *         multiplying each half of the dividend by
 *         the reciprocal.
 *      4) Truncate the approximate quotient to an
 *         integer.  Due to rounding, this value may be
 *         smaller than the actual quotient by one.
 *      5) Compute the remainder based on the approximate
 *         quotient.
 *      6) If necessary, increment quotient and adjust
 *         remainder.  There are two cases of an incorrect
 *         quotient:
 *              a) The computed remainder is greater than
 *                 or equal to the divisor.
 *              b) The computed remainder is smaller than
 *                 an approximate remainder computed by
 *                 multiplying the fractional part of the
 *                 quotient with the divisor.
 */
/************************************************************************/
/*  WE PROVIDE THIS SOURCE CODE EXAMPLE "AS IS," WITHOUT WARRANTY OF    */
/*  ANY KIND EITHER EXPRESSED OR IMPLIED INCLUDING BUT NOT LIMITED TO   */
/*  THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A         */
/*  PARTICULAR PURPOSE.  The entire risk as to the quality and          */
/*  performance of this example is with you.                            */
/************************************************************************/

/*
 *
 *      Values for the IEEE Rounding Modes (ANSI Encoding)
 *
 *      RZ = Round toward zero
 *      RN = Round toward nearest (default)
 *      RP = Round toward plus infinity
 *      RM = Round toward minus infinity
 *
 */
#define FP_RND_RZ       0
#define FP_RND_RN       1
#define FP_RND_RP       2
#define FP_RND_RM       3


#define TWO32 ((double)(1<<16)*(double)(1<<16))
#define TWO52 ((double)(1<<26)*(double)(1<<26))

#define uint unsigned int

#define D2INT(DOUBLE,INT) { double convert = (DOUBLE) + TWO52; \
                            INT = ((uint *)(&convert))[1]; }


int divu_64_32( uint D, uint H, uint L, uint *RP, uint *QP )
{
	uint Q;

	/* Determine if quotient fits in 32-bits */
	if (H >= D)
		return 1;         /* Overflow */
	else {
		unsigned short save = fp_swap_rnd( FP_RND_RZ );
		double rd = 1 / (double)D;
		double Q_f = rd * (double)H * TWO32 + rd * (double)L;
		D2INT( Q_f, Q );
		{
			uint R_aprox;
			uint R_qd = L - Q * D;
			if (R_qd >= D)
				Q++, *RP = R_qd-D;
			else {
				D2INT((Q_f - (double)Q) * (double)D, R_aprox);
				if (R_aprox > R_qd)
					Q++, *RP = R_qd-D;
				else
					*RP = R_qd;
			}
		}
		(void) fp_swap_rnd( save );
	}

	*QP = Q;
	return 0;
}

/* The rest of this is written by Joe Orost (joe@babel.ho.att.com) */

#define sign(a) ((a) < 0)
#define abs(a)  (sign(a) ? -(a) : (a))

int
divs_64_32(d, n_0, n_1, RP, QP)
long d, n_0;
unsigned long n_1;
long *RP, *QP;
{
    register long s_n, s_d, ov;
    unsigned long a_n_0, a_d;
    register long long nn, a_nn;			/* Compile with gcc */

    nn = (((long long)(n_0)) << 32) + n_1;
    s_n = sign(nn);
    a_nn = abs(nn);
    a_n_0 = (long)(a_nn >> 32);
    n_1 = (long)(a_nn & 0xffffffff);
    s_d = sign(d);
    a_d = abs(d);
    ov = divu_64_32(a_d, a_n_0, n_1, RP, QP);
    if(ov)
        return 1;
    if(s_n != s_d) {
        *QP = -(*QP);
        if(*QP > 0)
            return 1;
    } else if(*QP < 0) {
        return 1;
    }
    if(s_n)
        *RP = -(*RP);
    return 0;
}