amd-builtins/math32/remainderF.h

/*
 * Copyright (c) 2014 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#include "math32.h"

#if !defined(SUBNORMALS_SUPPORTED)
static inline float
scaleFullRangef32(float y, float t)
{
    float ay, ty, r = 0;
    int k, iiy, iy, exp_iy0, exp_iy, manty, signy, miy;
    int delta, shift, ir;

    ay = fabs(t);
    k = ay > 1024 ? 1024 : (int) ay;
    k = t < 0 ? -k : k;
    t = (float) k;

    iiy = as_int(y);
    iy = iiy & EXSIGNBIT_SP32;
    signy = iiy & SIGNBIT_SP32;
    ay = as_float(iy);

    exp_iy0 = iy & EXPBITS_SP32;
    manty = iy & MANTBITS_SP32;

    //sub-normal
    ty = exp_iy0 == 0 ? (float) manty : as_float(iy);
    k = exp_iy0 == 0 ? k - 149 : k;
    ay = ty;
    iy = as_int(ay);
    exp_iy0 = iy & EXPBITS_SP32;
    exp_iy = (exp_iy0 >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
    // add k to y's exponent
    r = as_float(iy + (k << EXPSHIFTBITS_SP32));
    r = (exp_iy + k) > 127 ? as_float(PINFBITPATT_SP32) : r;
    // add k to y's exponent
    delta = -126 - (exp_iy + k);

    // sub-normal
    miy = iy & MANTBITS_SP32;
    miy |= IMPBIT_SP32;
    shift = delta > 23 ? 24 : delta;
    shift = delta < 0 ? 0 : shift;
    miy >>= shift;
    r = delta > 0 ? as_float(miy) : r;
    r = t > (float) (2 * EMAX_SP32) ? as_float(PINFBITPATT_SP32) : r;
    ir = as_int(r);
    r = ir <= PINFBITPATT_SP32 ? as_float(as_int(r) | signy) : r;
    return r;
}

/* Scales the float x by 2.0**n.
Assumes 2*EMIN <= n <= 2*EMAX, though this condition is not checked. */
static inline float
scaleFloat_2(float x, int n)
{
    float t1, t2;
    int n1, n2;
    n1 = n / 2;
    n2 = n - n1;
    /* Construct the numbers t1 = 2.0**n1 and t2 = 2.0**n2 */
    t1 = as_float((n1 + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
    t2 = as_float((n2 + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
    return (x * t1) * t2;
}

/* Scales the float x by 2.0**n.
   Assumes EMIN <= n <= EMAX, though this condition is not checked. */
static inline float
scaleFloat_1(float x, int n)
{
    float t;
    /* Construct the number t = 2.0**n */
    t = as_float((n + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
    return x * t;
}

/* Computes the exact product of x and y, the result being the
nearly double length number (z,zz) */
static inline void
mul12f(float x, float y, float *z, float *zz)
{
    float hx, tx, hy, ty;
    // Split x into hx (head) and tx (tail). Do the same for y.
    uint u;
    u = as_uint(x);
    u &= 0xfffff000;
    hx = as_float(u);
    tx = x - hx;
    u = as_uint(y);
    u &= 0xfffff000;
    hy = as_float(u);
    ty = y - hy;
    *z = x * y;
    *zz = (((hx * hy - *z) + hx * ty) + tx * hy) + tx * ty;
}

#endif //SUBNORMALS_SUPPORTED

#if defined(COMPILING_FMOD)
__attribute__((overloadable)) float
fmod(float x, float y)
#elif defined(COMPILING_REMQUO)
__attribute__((overloadable)) float
remquo(float x, float y, int *quo)
#else
__attribute__((overloadable)) float
remainder(float x, float y)
#endif
{
#if !defined(SUBNORMALS_SUPPORTED)

    const int loop_scale = 12;
    const float fscale = 1.0f / (float) (1 << loop_scale);

    int ntimes;
    float ret = 0;
    int ui_x, ui_y, ui_ax, ui_ay, xexp, yexp, signx;
    float af_x, af_y, af_ybase, fx, fxp, fxm, fy, w, scale, t, c, cc, v;
    float yscale, scaled_w, saved_w, div, sdiv, ratio, sratio, fxexp, sub_fx;
    int iw_scaled, wexp, it, i, ifx, ex, ey;;
    float xr, xr0, xr_base, yr;
    uint q;

    ui_x = as_int(x);
    ui_y = as_int(y);
    ui_ax = ui_x & EXSIGNBIT_SP32;
    ui_ay = ui_y & EXSIGNBIT_SP32;

    /* special case handle */
#if defined(COMPILING_REMQUO)
    *quo = 0;
#endif
    if (ui_ax > PINFBITPATT_SP32)
	return x;
    if (ui_ax == PINFBITPATT_SP32)
	return as_float(QNANBITPATT_SP32);
    if (ui_ay > PINFBITPATT_SP32)
	return y;
    if (ui_ay == PINFBITPATT_SP32)
	return x;
    if (ui_ay == 0 && ui_ax == 0)
	return as_float(QNANBITPATT_SP32);
    if (ui_ax == 0)
	return x;
    if (ui_ay == 0)
	return as_float(QNANBITPATT_SP32);

    signx = ui_x & SIGNBIT_SP32;
#if defined(COMPILING_REMQUO)
    int signy = ui_y & SIGNBIT_SP32;
#endif
    af_x = as_float(ui_ax);
    af_ybase = af_y = as_float(ui_ay);
    yexp = (int) ((ui_y & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);

    yscale = (float) ((yexp < 48 && ui_ay != 0) ? (48 - yexp) : 0);
    if (yscale != 0) {
	af_y = scaleFullRangef32(af_ybase, yscale);
    }

    ui_y = as_int(af_y);
    yexp = (int) ((ui_y & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
    xexp = (int) ((ui_x & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
    fx = af_x;
    fy = af_y;

    /* Set ntimes to the number of times we need to do a
       partial remainder. If the exponent of x is an exact multiple
       of 24 larger than the exponent of y, and the mantissa of x is
       less than the mantissa of y, ntimes will be one too large
       but it doesn't matter - it just means that we'll go round
       the loop below one extra time. */
    ntimes = (xexp - yexp) / loop_scale;
    ntimes = xexp <= yexp ? 0 : ntimes;

    /* Set w = y * 2^(ntimes*loop_scale) */
    w = scaleFloat_2(fy, ntimes*loop_scale);
    w = ntimes == 0 ? fy : w;

    /* Set scale = 2^(-loop_scale) */
    scale = ntimes == 0 ? 1.0f : fscale;

    // make sure recip does not overflow
    wexp = (int) ((as_int(w) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
    saved_w = w;
    scaled_w = scaleFloat_1(w, -14);
    iw_scaled = wexp > 105 & wexp <= 127;
    w = iw_scaled & ntimes > 0 ? scaled_w : w;

    /* Each time round the loop we compute a partial remainder.
       This is done by subtracting a large multiple of w
       from x each time, where w is a scaled up version of y.
       The subtraction can be performed exactly when performed
       in double precision, and the result at each stage can
       fit exactly in a single precision number. */
    for (i = 0; i < ntimes; i++) {
	/* Set fx = fx - w * t, where t is equal to trunc(dx/w). */
	div = __amdil_improved_div_f32(fx, w);
	sdiv = scaleFloat_1(div, -14);
	div = iw_scaled ? sdiv : div;
	t = floor(div);
	w = saved_w;
	iw_scaled = 0;

	/* At this point, t may be one too large due to rounding of fx/w */

	/* Compute w * t in quad precision */
	mul12f(w, t, &c, &cc);

	/* Subtract w * t from fx */
	v = fx - c;
	fx = v + (((fx - v) - c) - cc);

	/* If t was one too large, fx will be negative. Add back one w */
	/* It might be possible to speed up this loop by finding
	   a way to compute correctly truncated t directly from fx and w.
	   We would then avoid the need for this check on negative fx. */
	fxp = fx + w;
	fxm = fx - w;
	fx = fx < 0.0f ? fxp : fx;
	fx = fx >= w ? fxm : fx;

	/* Scale w down by for the next iteration */
	w *= scale;
	saved_w = w;
    }

    /* One more time */
    // iw = as_int(w);
    ifx = as_int(fx);
    fxexp = (int) ((ifx & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
    // wexp = (int) ((iw & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
    sub_fx = fx;
    // make sure recip does not overflow
    wexp = (int) ((as_int(w) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
    saved_w = w;
    scaled_w = scaleFloat_1(w, -14);
    iw_scaled = wexp > 105 & wexp <= 127;
    w = iw_scaled ? scaled_w : w;
    ratio = __amdil_improved_div_f32(fx, w);
    sratio = scaleFloat_1(ratio, -14);
    ratio = iw_scaled ? sratio : ratio;
    t = floor(ratio);
    it = (int) t;

    w = saved_w;
    mul12f(w, t, &c, &cc);

    v = fx - c;
    fx = v + (((fx - v) - c) - cc);

    if (fx < 0.0f) {
	fx += w;
	it--;
    }

    if (fx >= w) {
	fx -= w;
	it++;
    }

    // sub-normal fax
    fx = fxexp == 0 ? sub_fx : fx;

#if !defined(COMPILING_FMOD)
    float scaleback = 0;
#endif

    // in case fx == 0 and we'got a divisor
    it = (yscale > 30) ? 0 : ((unsigned int) it << (int) yscale);

    if (as_int(fx) != 0 && yscale != 0) {
	xr = fx;
	xr_base = fx;
	yr = af_ybase;
	q = 0;
	ex = ilogb(fx);
	ey = ilogb(af_ybase);

	yr = (float) scaleFullRangef32(af_ybase, (float) -ey);
	xr = (float) scaleFullRangef32(fx, (float) -ex);

	for (i = ex - ey; i > 0; i--) {
	    q <<= 1;
	    xr0 = xr;
	    xr = (xr0 >= yr) ? xr0 - yr : xr0;
	    q = (xr0 >= yr) ? q + 1 : q;
	    xr += xr;
	}
	q <<= 1;
	xr0 = xr;
	xr = (xr0 >= yr) ? xr0 - yr : xr0;
	q = (xr0 >= yr) ? q + 1 : q;
	xr = scaleFullRangef32(xr, (float) ey);

	fx = (ex - ey >= 0) ? xr : xr_base;
#if !defined(COMPILING_FMOD)
	q = (ex - ey >= 0) ? q : 0;
	it += q;

	xexp = (int) ((as_int(fx) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);

	w = af_ybase;
	if (xexp < 24) {
	    fx = scaleFullRangef32(fx, 48);
	    w = scaleFullRangef32(af_ybase, 48);
	    scaleback = -48;
	}
#endif
    }
#if !defined(COMPILING_FMOD)
    /* At this point, dx lies in the range [0,dy) */
    /* For the remainder function, we need to adjust dx
       so that it lies in the range (-y/2, y/2] by carefully
       subtracting w (== fy == y) if necessary. */
    if (fx * 2.f > w || ((fx * 2.f == w) && (it & 1))) {
	fx -= w;
	it++;
    }
    if (scaleback != 0) {
	fx = scaleFullRangef32(fx, scaleback);
    }
#endif

    ret = (signx) ? as_float(as_int(fx) ^ SIGNBIT_SP32) : fx;
#if defined(COMPILING_REMQUO)
    it = (signx ^ signy) ? -it : it;
    *quo = it;
#endif

    return ret;


#else

    x = FTZ(x);
    y = FTZ(y);

    int ux = as_int(x);
    int ax = ux & EXSIGNBIT_SP32;
    float xa = as_float(ax);
    int sx = ux ^ ax;
    int ex = ax >> EXPSHIFTBITS_SP32;

    int uy = as_int(y);
    int ay = uy & EXSIGNBIT_SP32;
    float ya = as_float(ay);
#if defined COMPILING_REMQUO
    int sy = uy ^ ay;
#endif
    int ey = ay >> EXPSHIFTBITS_SP32;

    float xr = as_float(0x3f800000 | (ax & 0x007fffff));
    float yr = as_float(0x3f800000 | (ay & 0x007fffff));
    int c;
    int k = ex - ey;

#if defined COMPILING_FMOD
# define BIT c = xr >= yr; xr -= c ? yr : 0.0f; xr += xr
#else
    uint q = 0;
# define BIT c = xr >= yr; q = (q << 1) | c; xr -= c ? yr : 0.0f; xr += xr
#endif

    while (k > 3) {
	BIT;
	BIT;
	BIT;
	BIT;
	k -= 4;
    }

    while (k > 0) {
	BIT;
	--k;
    }

#if !defined COMPILING_FMOD
    c = xr > yr;
    q = (q << 1) | c;
#else
    c = xr >= yr;
#endif
    xr -= c ? yr : 0.0f;

    int lt = ex < ey;

#if !defined COMPILING_FMOD
    q = lt ? 0 : q;
#endif
    xr = lt ? xa : xr;
    yr = lt ? ya : yr;

#if !defined COMPILING_FMOD
    c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & (q & 0x1) == 0x1);
    xr -= c ? yr : 0.0f;
    q += c;
#endif

    float s = as_float(ey << EXPSHIFTBITS_SP32);
    xr *= lt ? 1.0f : s;

#if defined COMPILING_REMQUO
    int qsgn = sx == sy ? 1 : -1;
    int quot = (q & 0x7f) * qsgn;
#endif

    c = ax == ay;
#if defined COMPILING_REMQUO
    quot = c ? qsgn : quot;
#endif
    xr = c ? 0.0f : xr;

    xr = as_float(sx ^ as_int(xr));

    c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | ay == 0;
#if defined COMPILING_REMQUO
    quot = c ? 0 : quot;
#endif
    xr = c ? as_float(QNANBITPATT_SP32) : xr;

#if defined COMPILING_REMQUO
    *quo = quot;
#endif

    return xr;

#endif
}