/*
 * SoftFP Library
 * 
 * Copyright (c) 2016 Fabrice Bellard
 *
 * 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.
 */
#if F_SIZE == 32
#define F_UINT uint32_t
#define F_ULONG uint64_t
#define MANT_SIZE 23
#define EXP_SIZE 8
#elif F_SIZE == 64
#define F_UHALF uint32_t
#define F_UINT uint64_t
#ifdef HAVE_INT128
#define F_ULONG uint128_t
#endif
#define MANT_SIZE 52
#define EXP_SIZE 11
#elif F_SIZE == 128
#define F_UHALF uint64_t
#define F_UINT uint128_t
#define MANT_SIZE 112
#define EXP_SIZE 15
#else
#error unsupported F_SIZE
#endif

#define EXP_MASK ((1 << EXP_SIZE) - 1)
#define MANT_MASK (((F_UINT)1 << MANT_SIZE) - 1)
#define SIGN_MASK ((F_UINT)1 << (F_SIZE - 1))
#define IMANT_SIZE (F_SIZE - 2) /* internal mantissa size */
#define RND_SIZE (IMANT_SIZE - MANT_SIZE)
#define QNAN_MASK ((F_UINT)1 << (MANT_SIZE - 1))

/* quiet NaN */
#define F_QNAN glue(F_QNAN, F_SIZE)
#define clz glue(clz, F_SIZE)
#define pack_sf glue(pack_sf, F_SIZE)
#define unpack_sf glue(unpack_sf, F_SIZE)
#define rshift_rnd glue(rshift_rnd, F_SIZE)
#define round_pack_sf glue(roundpack_sf, F_SIZE)
#define normalize_sf glue(normalize_sf, F_SIZE)
#define normalize2_sf glue(normalize2_sf, F_SIZE)
#define issignan_sf glue(issignan_sf, F_SIZE)
#define isnan_sf glue(isnan_sf, F_SIZE)
#define add_sf glue(add_sf, F_SIZE)
#define mul_sf glue(mul_sf, F_SIZE)
#define fma_sf glue(fma_sf, F_SIZE)
#define div_sf glue(div_sf, F_SIZE)
#define sqrt_sf glue(sqrt_sf, F_SIZE)
#define normalize_subnormal_sf glue(normalize_subnormal_sf, F_SIZE)
#define divrem_u glue(divrem_u, F_SIZE)
#define sqrtrem_u glue(sqrtrem_u, F_SIZE)
#define mul_u glue(mul_u, F_SIZE)
#define cvt_sf32_sf glue(cvt_sf32_sf, F_SIZE)
#define cvt_sf64_sf glue(cvt_sf64_sf, F_SIZE)

static const F_UINT F_QNAN = (((F_UINT)EXP_MASK << MANT_SIZE) | ((F_UINT)1 << (MANT_SIZE - 1)));

static inline F_UINT pack_sf(uint32_t a_sign, uint32_t a_exp, F_UINT a_mant)
{
    return ((F_UINT)a_sign << (F_SIZE - 1)) |
        ((F_UINT)a_exp << MANT_SIZE) | 
        (a_mant & MANT_MASK);
}

static inline F_UINT unpack_sf(uint32_t *pa_sign, int32_t *pa_exp,
                               F_UINT a)
{
    *pa_sign = a >> (F_SIZE - 1);
    *pa_exp = (a >> MANT_SIZE) & EXP_MASK;
    return a & MANT_MASK;
} 

static F_UINT rshift_rnd(F_UINT a, int d)
{
    F_UINT mask;
    if (d != 0) {
        if (d >= F_SIZE) {
            a = (a != 0);
        } else {
            mask = ((F_UINT)1 << d) - 1;
            a = (a >> d) | ((a & mask) != 0);
        }
    }
    return a;
}

/* a_mant is considered to have its MSB at F_SIZE - 2 bits */
static F_UINT round_pack_sf(uint32_t a_sign, int a_exp, F_UINT a_mant,
                            RoundingModeEnum rm, uint32_t *pfflags)
{
    int diff;
    uint32_t addend, rnd_bits;

    switch(rm) {
    case RM_RNE:
    case RM_RMM:
        addend = (1 << (RND_SIZE - 1));
        break;
    case RM_RTZ:
        addend = 0;
        break;
    default:
    case RM_RDN:
    case RM_RUP:
        //        printf("s=%d rm=%d m=%x\n", a_sign, rm, a_mant);
        if (a_sign ^ (rm & 1))
            addend = (1 << RND_SIZE) - 1;
        else
            addend = 0;
        break;
    }

    /* potentially subnormal */
    if (a_exp <= 0) {
        BOOL is_subnormal;
        /* Note: we set the underflow flag if the rounded result
           is subnormal and inexact */
        is_subnormal = (a_exp < 0 || 
                        (a_mant + addend) < ((F_UINT)1 << (F_SIZE - 1)));
        diff = 1 - a_exp;
        a_mant = rshift_rnd(a_mant, diff);
        rnd_bits = a_mant & ((1 << RND_SIZE ) - 1);
        if (is_subnormal && rnd_bits != 0) {
            *pfflags |= FFLAG_UNDERFLOW;
        }
        a_exp = 1;
    } else {
        rnd_bits = a_mant & ((1 << RND_SIZE ) - 1);
    }
    if (rnd_bits != 0)
        *pfflags |= FFLAG_INEXACT;
    a_mant = (a_mant + addend) >> RND_SIZE;
    /* half way: select even result */
    if (rm == RM_RNE && rnd_bits == (1 << (RND_SIZE - 1)))
        a_mant &= ~1;
    /* Note the rounding adds at least 1, so this is the maximum
       value */
    a_exp += a_mant >> (MANT_SIZE + 1);
    if (a_mant <= MANT_MASK) {
        /* denormalized or zero */
        a_exp = 0;
    } else if (a_exp >= EXP_MASK) {
        /* overflow */
        if (addend == 0) {
            a_exp = EXP_MASK - 1;
            a_mant = MANT_MASK;
        } else {
            /* infinity */
            a_exp = EXP_MASK;
            a_mant = 0;
        }
        *pfflags |= FFLAG_OVERFLOW | FFLAG_INEXACT;
    }
    return pack_sf(a_sign, a_exp, a_mant);
}

/* a_mant is considered to have at most F_SIZE - 1 bits */
static F_UINT normalize_sf(uint32_t a_sign, int a_exp, F_UINT a_mant,
                           RoundingModeEnum rm, uint32_t *pfflags)
{
    int shift;
    shift = clz(a_mant) - (F_SIZE - 1 - IMANT_SIZE);
    assert(shift >= 0);
    a_exp -= shift;
    a_mant <<= shift;
    return round_pack_sf(a_sign, a_exp, a_mant, rm, pfflags);
}

/* same as normalize_sf() but with a double word mantissa. a_mant1 is
   considered to have at most F_SIZE - 1 bits */
static F_UINT normalize2_sf(uint32_t a_sign, int a_exp, F_UINT a_mant1, F_UINT a_mant0,
                            RoundingModeEnum rm, uint32_t *pfflags)
{
    int l, shift;
    if (a_mant1 == 0) {
        l = F_SIZE + clz(a_mant0);
    } else {
        l = clz(a_mant1);
    }
    shift = l - (F_SIZE - 1 - IMANT_SIZE);
    assert(shift >= 0);
    a_exp -= shift;
    if (shift == 0) {
        a_mant1 |= (a_mant0 != 0);
    } else if (shift < F_SIZE) {
        a_mant1 = (a_mant1 << shift) | (a_mant0 >> (F_SIZE - shift));
        a_mant0 <<= shift;
        a_mant1 |= (a_mant0 != 0);
    } else {
        a_mant1 = a_mant0 << (shift - F_SIZE);
    }
    return round_pack_sf(a_sign, a_exp, a_mant1, rm, pfflags);
}

BOOL issignan_sf(F_UINT a)
{
    uint32_t a_exp1;
    F_UINT a_mant;
    a_exp1 = (a >> (MANT_SIZE - 1)) & ((1 << (EXP_SIZE + 1)) - 1);
    a_mant = a & MANT_MASK;
    return (a_exp1 == (2 * EXP_MASK) && a_mant != 0);
}

BOOL isnan_sf(F_UINT a)
{
    uint32_t a_exp;
    F_UINT a_mant;
    a_exp = (a >> MANT_SIZE) & EXP_MASK;
    a_mant = a & MANT_MASK;
    return (a_exp == EXP_MASK && a_mant != 0);
}


F_UINT add_sf(F_UINT a, F_UINT b, RoundingModeEnum rm,
              uint32_t *pfflags)
{
    uint32_t a_sign, b_sign, a_exp, b_exp;
    F_UINT tmp, a_mant, b_mant;

    /* swap so that  abs(a) >= abs(b) */
    if ((a & ~SIGN_MASK) < (b & ~SIGN_MASK)) {
        tmp = a;
        a = b;
        b = tmp;
    }
    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);
    a_exp = (a >> MANT_SIZE) & EXP_MASK;
    b_exp = (b >> MANT_SIZE) & EXP_MASK;
    a_mant = (a & MANT_MASK) << 3;
    b_mant = (b & MANT_MASK) << 3;
    if (unlikely(a_exp == EXP_MASK)) {
        if (a_mant != 0) {
            /* NaN result */
            if (!(a_mant & (QNAN_MASK << 3)) || issignan_sf(b))
                *pfflags |= FFLAG_INVALID_OP;
            return F_QNAN;
        } else if (b_exp == EXP_MASK && a_sign != b_sign) {
            *pfflags |= FFLAG_INVALID_OP;
            return F_QNAN;
        } else {
            /* infinity */
            return a;
        }
    }
    if (a_exp == 0) {
        a_exp = 1;
    } else {
        a_mant |= (F_UINT)1 << (MANT_SIZE + 3);
    }
    if (b_exp == 0) {
        b_exp = 1;
    } else {
        b_mant |= (F_UINT)1 << (MANT_SIZE + 3);
    }
    b_mant = rshift_rnd(b_mant, a_exp - b_exp);
    if (a_sign == b_sign) {
        /* same signs : add the absolute values  */
        a_mant += b_mant;
    } else {
        /* different signs : subtract the absolute values  */
        a_mant -= b_mant;
        if (a_mant == 0) {
            /* zero result : the sign needs a specific handling */
            a_sign = (rm == RM_RDN);
        }
    }
    a_exp += (RND_SIZE - 3);
    return normalize_sf(a_sign, a_exp, a_mant, rm, pfflags);
}

F_UINT glue(sub_sf, F_SIZE)(F_UINT a, F_UINT b, RoundingModeEnum rm,
                            uint32_t *pfflags)
{
    return add_sf(a, b ^ SIGN_MASK, rm, pfflags);
}

static inline F_UINT normalize_subnormal_sf(int32_t *pa_exp, F_UINT a_mant)
{
    int shift;
    shift = MANT_SIZE - ((F_SIZE - 1 - clz(a_mant)));
    *pa_exp = 1 - shift;
    return a_mant << shift;
}

#ifdef F_ULONG

static F_UINT mul_u(F_UINT *plow, F_UINT a, F_UINT b)
{
    F_ULONG r;
    r = (F_ULONG)a * (F_ULONG)b;
    *plow = r;
    return r >> F_SIZE;
}

#else

#define FH_SIZE (F_SIZE / 2)

static F_UINT mul_u(F_UINT *plow, F_UINT a, F_UINT b)
{
    F_UHALF a0, a1, b0, b1, r0, r1, r2, r3;
    F_UINT r00, r01, r10, r11, c;
    a0 = a;
    a1 = a >> FH_SIZE;
    b0 = b;
    b1 = b >> FH_SIZE;

    r00 = (F_UINT)a0 * (F_UINT)b0;
    r01 = (F_UINT)a0 * (F_UINT)b1;
    r10 = (F_UINT)a1 * (F_UINT)b0;
    r11 = (F_UINT)a1 * (F_UINT)b1;
    
    r0 = r00;
    c = (r00 >> FH_SIZE) + (F_UHALF)r01 + (F_UHALF)r10;
    r1 = c;
    c = (c >> FH_SIZE) + (r01 >> FH_SIZE) + (r10 >> FH_SIZE) + (F_UHALF)r11;
    r2 = c;
    r3 = (c >> FH_SIZE) + (r11 >> FH_SIZE);

    *plow = ((F_UINT)r1 << FH_SIZE) | r0;
    return ((F_UINT)r3 << FH_SIZE) | r2;
}

#undef FH_SIZE

#endif

F_UINT mul_sf(F_UINT a, F_UINT b, RoundingModeEnum rm,
              uint32_t *pfflags)
{
    uint32_t a_sign, b_sign, r_sign;
    int32_t a_exp, b_exp, r_exp;
    F_UINT a_mant, b_mant, r_mant, r_mant_low;

    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);
    r_sign = a_sign ^ b_sign;
    a_exp = (a >> MANT_SIZE) & EXP_MASK;
    b_exp = (b >> MANT_SIZE) & EXP_MASK;
    a_mant = a & MANT_MASK;
    b_mant = b & MANT_MASK;
    if (a_exp == EXP_MASK || b_exp == EXP_MASK) {
        if (isnan_sf(a) || isnan_sf(b)) {
            if (issignan_sf(a) || issignan_sf(b)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN;
        } else {
            /* infinity */
            if ((a_exp == EXP_MASK && (b_exp == 0 && b_mant == 0)) ||
                (b_exp == EXP_MASK && (a_exp == 0 && a_mant == 0))) {
                *pfflags |= FFLAG_INVALID_OP;
                return F_QNAN;
            } else {
                return pack_sf(r_sign, EXP_MASK, 0);
            }
        }
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            return pack_sf(r_sign, 0, 0); /* zero */
        a_mant = normalize_subnormal_sf(&a_exp, a_mant);
    } else {
        a_mant |= (F_UINT)1 << MANT_SIZE;
    }
    if (b_exp == 0) {
        if (b_mant == 0)
            return pack_sf(r_sign, 0, 0); /* zero */
        b_mant = normalize_subnormal_sf(&b_exp, b_mant);
    } else {
        b_mant |= (F_UINT)1 << MANT_SIZE;
    }
    r_exp = a_exp + b_exp - (1 << (EXP_SIZE - 1)) + 2;
    
    r_mant = mul_u(&r_mant_low,a_mant << RND_SIZE, b_mant << (RND_SIZE + 1));
    r_mant |= (r_mant_low != 0);
    return normalize_sf(r_sign, r_exp, r_mant, rm, pfflags);
}

/* fused multiply and add */
F_UINT fma_sf(F_UINT a, F_UINT b, F_UINT c, RoundingModeEnum rm,
              uint32_t *pfflags)
{
    uint32_t a_sign, b_sign, c_sign, r_sign;
    int32_t a_exp, b_exp, c_exp, r_exp, shift;
    F_UINT a_mant, b_mant, c_mant, r_mant1, r_mant0, c_mant1, c_mant0, mask;

    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);
    c_sign = c >> (F_SIZE - 1);
    r_sign = a_sign ^ b_sign;
    a_exp = (a >> MANT_SIZE) & EXP_MASK;
    b_exp = (b >> MANT_SIZE) & EXP_MASK;
    c_exp = (c >> MANT_SIZE) & EXP_MASK;
    a_mant = a & MANT_MASK;
    b_mant = b & MANT_MASK;
    c_mant = c & MANT_MASK;
    if (a_exp == EXP_MASK || b_exp == EXP_MASK || c_exp == EXP_MASK) {
        if (isnan_sf(a) || isnan_sf(b) || isnan_sf(c)) {
            if (issignan_sf(a) || issignan_sf(b) || issignan_sf(c)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN;
        } else {
            /* infinities */
            if ((a_exp == EXP_MASK && (b_exp == 0 && b_mant == 0)) ||
                (b_exp == EXP_MASK && (a_exp == 0 && a_mant == 0)) ||
                ((a_exp == EXP_MASK || b_exp == EXP_MASK) &&
                 (c_exp == EXP_MASK && r_sign != c_sign))) {
                *pfflags |= FFLAG_INVALID_OP;
                return F_QNAN;
            } else if (c_exp == EXP_MASK) {
                return pack_sf(c_sign, EXP_MASK, 0);
            } else {
                return pack_sf(r_sign, EXP_MASK, 0);
            }
        }
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            goto mul_zero;
        a_mant = normalize_subnormal_sf(&a_exp, a_mant);
    } else {
        a_mant |= (F_UINT)1 << MANT_SIZE;
    }
    if (b_exp == 0) {
        if (b_mant == 0) {
        mul_zero:
            if (c_exp == 0 && c_mant == 0) {
                if (c_sign != r_sign)
                    r_sign = (rm == RM_RDN);
                return pack_sf(r_sign, 0, 0);
            } else {
                return c;
            }
        }
        b_mant = normalize_subnormal_sf(&b_exp, b_mant);
    } else {
        b_mant |= (F_UINT)1 << MANT_SIZE;
    }
    /* multiply */
    r_exp = a_exp + b_exp - (1 << (EXP_SIZE - 1)) + 3;
    
    r_mant1 = mul_u(&r_mant0, a_mant << RND_SIZE, b_mant << RND_SIZE);
    /* normalize to F_SIZE - 3 */
    if (r_mant1 < ((F_UINT)1 << (F_SIZE - 3))) {
        r_mant1 = (r_mant1 << 1) | (r_mant0 >> (F_SIZE - 1));
        r_mant0 <<= 1;
        r_exp--;
    }

    /* add */
    if (c_exp == 0) {
        if (c_mant == 0) {
            /* add zero */
            r_mant1 |= (r_mant0 != 0);
            return normalize_sf(r_sign, r_exp, r_mant1, rm, pfflags);
        }
        c_mant = normalize_subnormal_sf(&c_exp, c_mant);
    } else {
        c_mant |= (F_UINT)1 << MANT_SIZE;
    }
    c_exp++;
    c_mant1 = c_mant << (RND_SIZE - 1);
    c_mant0 = 0;

    //    printf("r_s=%d r_exp=%d r_mant=%08x %08x\n", r_sign, r_exp, (uint32_t)r_mant1, (uint32_t)r_mant0);
    //    printf("c_s=%d c_exp=%d c_mant=%08x %08x\n", c_sign, c_exp, (uint32_t)c_mant1, (uint32_t)c_mant0);

    /* ensure that abs(r) >= abs(c) */
    if (!(r_exp > c_exp || (r_exp == c_exp && r_mant1 >= c_mant1))) {
        F_UINT tmp;
        int32_t c_tmp;
        /* swap */
        tmp = r_mant1; r_mant1 = c_mant1; c_mant1 = tmp;
        tmp = r_mant0; r_mant0 = c_mant0; c_mant0 = tmp;
        c_tmp = r_exp; r_exp = c_exp; c_exp = c_tmp;
        c_tmp = r_sign; r_sign = c_sign; c_sign = c_tmp;
    }
    /* right shift c_mant */
    shift = r_exp - c_exp;
    if (shift >= 2 * F_SIZE) {
        c_mant0 = (c_mant0 | c_mant1) != 0;
        c_mant1 = 0;
    } else if (shift >= F_SIZE + 1) {
        c_mant0 = rshift_rnd(c_mant1, shift - F_SIZE);
        c_mant1 = 0;
    } else if (shift == F_SIZE) {
        c_mant0 = c_mant1 | (c_mant0 != 0);
        c_mant1 = 0;
    } else if (shift != 0) {
        mask = ((F_UINT)1 << shift) - 1;
        c_mant0 = (c_mant1 << (F_SIZE - shift)) | (c_mant0 >> shift) | ((c_mant0 & mask) != 0);
        c_mant1 = c_mant1 >> shift;
    }
    //    printf("  r_mant=%08x %08x\n", (uint32_t)r_mant1, (uint32_t)r_mant0);
    //    printf("  c_mant=%08x %08x\n", (uint32_t)c_mant1, (uint32_t)c_mant0);
    /* add or subtract */
    if (r_sign == c_sign) {
        r_mant0 += c_mant0;
        r_mant1 += c_mant1 + (r_mant0 < c_mant0);
    } else {
        F_UINT tmp;
        tmp = r_mant0;
        r_mant0 -= c_mant0;
        r_mant1 = r_mant1 - c_mant1 - (r_mant0 > tmp);
        if ((r_mant0 | r_mant1) == 0) {
            /* zero result : the sign needs a specific handling */
            r_sign = (rm == RM_RDN);
        }
    }
#if 0
    //    printf("  r1_mant=%08x %08x\n", (uint32_t)r_mant1, (uint32_t)r_mant0);
    /* normalize */
    if (r_mant1 == 0) {
        r_mant1 = r_mant0;
        r_exp -= F_SIZE;
    } else {
        shift = clz(r_mant1) - (F_SIZE - 1 - IMANT_SIZE);
        if (shift != 0) {
            r_mant1 = (r_mant1 << shift) | (r_mant0 >> (F_SIZE - shift));
            r_mant0 <<= shift;
            r_exp -= shift;
        }
        r_mant1 |= (r_mant0 != 0);
    }
    return normalize_sf(r_sign, r_exp, r_mant1, rm, pfflags);
#else
    return normalize2_sf(r_sign, r_exp, r_mant1, r_mant0, rm, pfflags);
#endif
}

#ifdef F_ULONG

static F_UINT divrem_u(F_UINT *pr, F_UINT ah, F_UINT al, F_UINT b)
{
    F_ULONG a;
    a = ((F_ULONG)ah << F_SIZE) | al;
    *pr = a % b;
    return a / b;
}

#else

/* XXX: optimize */
static F_UINT divrem_u(F_UINT *pr, F_UINT a1, F_UINT a0, F_UINT b)
{
    int i, qb, ab;

    assert(a1 < b);
    for(i = 0; i < F_SIZE; i++) {
        ab = a1 >> (F_SIZE - 1);
        a1 = (a1 << 1) | (a0 >> (F_SIZE - 1));
        if (ab || a1 >= b) {
            a1 -= b;
            qb = 1;
        } else {
            qb = 0;
        }
        a0 = (a0 << 1) | qb;
    }
    *pr = a1;
    return a0;
}

#endif

F_UINT div_sf(F_UINT a, F_UINT b, RoundingModeEnum rm,
              uint32_t *pfflags)
{
    uint32_t a_sign, b_sign, r_sign;
    int32_t a_exp, b_exp, r_exp;
    F_UINT a_mant, b_mant, r_mant, r;

    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);
    r_sign = a_sign ^ b_sign;
    a_exp = (a >> MANT_SIZE) & EXP_MASK;
    b_exp = (b >> MANT_SIZE) & EXP_MASK;
    a_mant = a & MANT_MASK;
    b_mant = b & MANT_MASK;
    if (a_exp == EXP_MASK) {
        if (a_mant != 0 || isnan_sf(b)) {
            if (issignan_sf(a) || issignan_sf(b)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN;
        } else if (b_exp == EXP_MASK) {
            *pfflags |= FFLAG_INVALID_OP;
            return F_QNAN;
        } else {
            return pack_sf(r_sign, EXP_MASK, 0);
        }
    } else if (b_exp == EXP_MASK) {
        if (b_mant != 0) {
            if (issignan_sf(a) || issignan_sf(b)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN;
        } else {
            return pack_sf(r_sign, 0, 0);
        }
    }

    if (b_exp == 0) {
        if (b_mant == 0) { 
            /* zero */
            if (a_exp == 0 && a_mant == 0) {
                *pfflags |= FFLAG_INVALID_OP;
                return F_QNAN;
            } else {
                *pfflags |= FFLAG_DIVIDE_ZERO;
                return pack_sf(r_sign, EXP_MASK, 0);
            }
        }
        b_mant = normalize_subnormal_sf(&b_exp, b_mant);
    } else {
        b_mant |= (F_UINT)1 << MANT_SIZE;
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            return pack_sf(r_sign, 0, 0); /* zero */
        a_mant = normalize_subnormal_sf(&a_exp, a_mant);
    } else {
        a_mant |= (F_UINT)1 << MANT_SIZE;
    }
    r_exp = a_exp - b_exp + (1 << (EXP_SIZE - 1)) - 1;
    r_mant = divrem_u(&r, a_mant, 0, b_mant << 2);
    if (r != 0)
        r_mant |= 1;
    return normalize_sf(r_sign, r_exp, r_mant, rm, pfflags);
}

#ifdef F_ULONG

/* compute sqrt(a) with a = ah*2^F_SIZE+al and a < 2^(F_SIZE - 2)
   return true if not exact square. */
static int sqrtrem_u(F_UINT *pr, F_UINT ah, F_UINT al)
{
    F_ULONG a, u, s;
    int l, inexact;

    /* 2^l >= a */
    if (ah != 0) {
        l = 2 * F_SIZE - clz(ah - 1);
    } else {
        if (al == 0) {
            *pr = 0;
            return 0;
        }
        l = F_SIZE - clz(al - 1);
    }
    a = ((F_ULONG)ah << F_SIZE) | al;
    u = (F_ULONG)1 << ((l + 1) / 2);
    for(;;) {
        s = u;
        u = ((a / s) + s) / 2;
        if (u >= s)
            break;
    }
    inexact = (a - s * s) != 0;
    *pr = s;
    return inexact;
}

#else

static int sqrtrem_u(F_UINT *pr, F_UINT a1, F_UINT a0)
{
    int l, inexact;
    F_UINT u, s, r, q, sq0, sq1;

    /* 2^l >= a */
    if (a1 != 0) {
        l = 2 * F_SIZE - clz(a1 - 1);
    } else {
        if (a0 == 0) {
            *pr = 0;
            return 0;
        }
        l = F_SIZE - clz(a0 - 1);
    }
    u = (F_UINT)1 << ((l + 1) / 2);
    for(;;) {
        s = u;
        q = divrem_u(&r, a1, a0, s);
        u = (q + s) / 2;
        if (u >= s)
            break;
    }
    sq1 = mul_u(&sq0, s, s);
    inexact = (sq0 != a0 || sq1 != a1);
    *pr = s;
    return inexact;
}

#endif

F_UINT sqrt_sf(F_UINT a, RoundingModeEnum rm,
               uint32_t *pfflags)
{
    uint32_t a_sign;
    int32_t a_exp;
    F_UINT a_mant;

    a_sign = a >> (F_SIZE - 1);
    a_exp = (a >> MANT_SIZE) & EXP_MASK;
    a_mant = a & MANT_MASK;
    if (a_exp == EXP_MASK) {
        if (a_mant != 0) {
            if (issignan_sf(a)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN;
        } else if (a_sign) {
            goto neg_error;
        } else {
            return a; /* +infinity */
        }
    }
    if (a_sign) {
        if (a_exp == 0 && a_mant == 0)
            return a; /* -zero */
    neg_error:
        *pfflags |= FFLAG_INVALID_OP;
        return F_QNAN;
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            return pack_sf(0, 0, 0); /* zero */
        a_mant = normalize_subnormal_sf(&a_exp, a_mant);
    } else {
        a_mant |= (F_UINT)1 << MANT_SIZE;
    }
    a_exp -= EXP_MASK / 2;
    /* simpler to handle an even exponent */
    if (a_exp & 1) {
        a_exp--;
        a_mant <<= 1;
    }
    a_exp = (a_exp >> 1) + EXP_MASK / 2;
    a_mant <<= (F_SIZE - 4 - MANT_SIZE);
    if (sqrtrem_u(&a_mant, a_mant, 0))
        a_mant |= 1;
    return normalize_sf(a_sign, a_exp, a_mant, rm, pfflags);
}

/* comparisons */

static F_UINT glue(min_max_nan_sf, F_SIZE)(F_UINT a, F_UINT b, uint32_t *pfflags, SoftFPMinMaxTypeEnum minmax_type)
{
    if (issignan_sf(a) || issignan_sf(b)) {
        *pfflags |= FFLAG_INVALID_OP;
        if (minmax_type == FMINMAX_IEEE754_2008)
            return F_QNAN;
    }
    if (minmax_type == FMINMAX_PROP) {
        return F_QNAN;
    } else {
        if (isnan_sf(a)) {
            if (isnan_sf(b)) 
                return F_QNAN;
            else
                return b;
        } else {
            return a;
        }
    }
}

F_UINT glue(min_sf, F_SIZE)(F_UINT a, F_UINT b, uint32_t *pfflags,
                            SoftFPMinMaxTypeEnum minmax_type)
{
    uint32_t a_sign, b_sign;

    if (isnan_sf(a) || isnan_sf(b)) {
        return glue(min_max_nan_sf, F_SIZE)(a, b, pfflags, minmax_type);
    }
    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);

    if (a_sign != b_sign) {
        if (a_sign)
            return a;
        else
            return b;
    } else {
        if ((a < b) ^ a_sign)
            return a;
        else
            return b;
    }
}

F_UINT glue(max_sf, F_SIZE)(F_UINT a, F_UINT b, uint32_t *pfflags,
                            SoftFPMinMaxTypeEnum minmax_type)
{
    uint32_t a_sign, b_sign;

    if (isnan_sf(a) || isnan_sf(b)) {
        return glue(min_max_nan_sf, F_SIZE)(a, b, pfflags, minmax_type);
    }
    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);

    if (a_sign != b_sign) {
        if (a_sign)
            return b;
        else
            return a;
    } else {
        if ((a < b) ^ a_sign)
            return b;
        else
            return a;
    }
}

int glue(eq_quiet_sf, F_SIZE)(F_UINT a, F_UINT b, uint32_t *pfflags)
{
    if (isnan_sf(a) || isnan_sf(b)) {
        if (issignan_sf(a) || issignan_sf(b)) {
            *pfflags |= FFLAG_INVALID_OP;
        }
        return 0;
    }

    if ((F_UINT)((a | b) << 1) == 0)
        return 1; /* zero case */
    return (a == b);
}

int glue(le_sf, F_SIZE)(F_UINT a, F_UINT b, uint32_t *pfflags)
{
    uint32_t a_sign, b_sign;

    if (isnan_sf(a) || isnan_sf(b)) {
        *pfflags |= FFLAG_INVALID_OP;
        return 0;
    }

    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);
    if (a_sign != b_sign) {
        return (a_sign || ((F_UINT)((a | b) << 1) == 0));
    } else {
        if (a_sign) {
            return (a >= b);
        } else {
            return (a <= b);
        }
    }
}

int glue(lt_sf, F_SIZE)(F_UINT a, F_UINT b, uint32_t *pfflags)
{
    uint32_t a_sign, b_sign;

    if (isnan_sf(a) || isnan_sf(b)) {
        *pfflags |= FFLAG_INVALID_OP;
        return 0;
    }

    a_sign = a >> (F_SIZE - 1);
    b_sign = b >> (F_SIZE - 1);
    if (a_sign != b_sign) {
        return (a_sign && ((F_UINT)((a | b) << 1) != 0));
    } else {
        if (a_sign) {
            return (a > b);
        } else {
            return (a < b);
        }
    }
}

uint32_t glue(fclass_sf, F_SIZE)(F_UINT a)
{
    uint32_t a_sign;
    int32_t a_exp;
    F_UINT a_mant;
    uint32_t ret;

    a_sign = a >> (F_SIZE - 1);
    a_exp = (a >> MANT_SIZE) & EXP_MASK;
    a_mant = a & MANT_MASK;
    if (a_exp == EXP_MASK) {
        if (a_mant != 0) {
            if (a_mant & QNAN_MASK)
                ret = FCLASS_QNAN;
            else
                ret = FCLASS_SNAN;
        } else {
            if (a_sign)
                ret = FCLASS_NINF;
            else
                ret = FCLASS_PINF;
        }
    } else if (a_exp == 0) {
        if (a_mant == 0) {
            if (a_sign)
                ret = FCLASS_NZERO;
            else
                ret = FCLASS_PZERO;
        } else {
            if (a_sign)
                ret = FCLASS_NSUBNORMAL;
            else
                ret = FCLASS_PSUBNORMAL;
        }
    } else {
        if (a_sign)
            ret = FCLASS_NNORMAL;
        else
            ret = FCLASS_PNORMAL;
    }
    return ret;
}

/* conversions between floats */

#if F_SIZE >= 64

F_UINT cvt_sf32_sf(uint32_t a, uint32_t *pfflags)
{
    uint32_t a_sign;
    int32_t a_exp;
    F_UINT a_mant;

    a_mant = unpack_sf32(&a_sign, &a_exp, a);
    if (a_exp == 0xff) {
        if (a_mant != 0) {
            /* NaN */
            if (issignan_sf32(a)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN;
        } else {
            /* infinity */
            return pack_sf(a_sign, EXP_MASK, 0);
        }
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            return pack_sf(a_sign, 0, 0); /* zero */
        a_mant = normalize_subnormal_sf32(&a_exp, a_mant);
    }
    /* convert the exponent value */
    a_exp = a_exp - 0x7f + (EXP_MASK / 2);
    /* shift the mantissa */
    a_mant <<= (MANT_SIZE - 23);
    /* We assume the target float is large enough to that no
       normalization is necessary */
    return pack_sf(a_sign, a_exp, a_mant);
}

uint32_t glue(glue(cvt_sf, F_SIZE), _sf32)(F_UINT a, RoundingModeEnum rm,
                                           uint32_t *pfflags)
{
    uint32_t a_sign;
    int32_t a_exp;
    F_UINT a_mant;

    a_mant = unpack_sf(&a_sign, &a_exp, a);
    if (a_exp == EXP_MASK) {
        if (a_mant != 0) {
            /* NaN */
            if (issignan_sf(a)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN32;
        } else {
            /* infinity */
            return pack_sf32(a_sign, 0xff, 0);
        }
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            return pack_sf32(a_sign, 0, 0); /* zero */
        normalize_subnormal_sf(&a_exp, a_mant);
    } else {
        a_mant |= (F_UINT)1 << MANT_SIZE;
    }
    /* convert the exponent value */
    a_exp = a_exp - (EXP_MASK / 2) + 0x7f;
    /* shift the mantissa */
    a_mant = rshift_rnd(a_mant, MANT_SIZE - (32 - 2));
    return normalize_sf32(a_sign, a_exp, a_mant, rm, pfflags);
}

#endif

#if F_SIZE >= 128

F_UINT cvt_sf64_sf(uint64_t a, uint32_t *pfflags)
{
    uint32_t a_sign;
    int32_t a_exp;
    F_UINT a_mant;

    a_mant = unpack_sf64(&a_sign, &a_exp, a);

    if (a_exp == 0x7ff) {
        if (a_mant != 0) {
            /* NaN */
            if (issignan_sf64(a)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN;
        } else {
            /* infinity */
            return pack_sf(a_sign, EXP_MASK, 0);
        }
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            return pack_sf(a_sign, 0, 0); /* zero */
        a_mant = normalize_subnormal_sf64(&a_exp, a_mant);
    }
    /* convert the exponent value */
    a_exp = a_exp - 0x3ff + (EXP_MASK / 2);
    /* shift the mantissa */
    a_mant <<= (MANT_SIZE - 52);
    return pack_sf(a_sign, a_exp, a_mant);
}

uint64_t glue(glue(cvt_sf, F_SIZE), _sf64)(F_UINT a, RoundingModeEnum rm,
                                               uint32_t *pfflags)
{
    uint32_t a_sign;
    int32_t a_exp;
    F_UINT a_mant;

    a_mant = unpack_sf(&a_sign, &a_exp, a);
    if (a_exp == EXP_MASK) {
        if (a_mant != 0) {
            /* NaN */
            if (issignan_sf(a)) {
                *pfflags |= FFLAG_INVALID_OP;
            }
            return F_QNAN64;
        } else {
            /* infinity */
            return pack_sf64(a_sign, 0x7ff, 0);
        }
    }
    if (a_exp == 0) {
        if (a_mant == 0)
            return pack_sf64(a_sign, 0, 0); /* zero */
        normalize_subnormal_sf(&a_exp, a_mant);
    } else {
        a_mant |= (F_UINT)1 << MANT_SIZE;
    }
    /* convert the exponent value */
    a_exp = a_exp - (EXP_MASK / 2) + 0x3ff;
    /* shift the mantissa */
    a_mant = rshift_rnd(a_mant, MANT_SIZE - (64 - 2));
    return normalize_sf64(a_sign, a_exp, a_mant, rm, pfflags);
}

#endif

#undef clz

#define ICVT_SIZE 32
#include "softfp_template_icvt.h"

#define ICVT_SIZE 64
#include "softfp_template_icvt.h"

#ifdef HAVE_INT128
#define ICVT_SIZE 128
#include "softfp_template_icvt.h"
#endif

#undef F_SIZE
#undef F_UINT
#undef F_ULONG
#undef F_UHALF
#undef MANT_SIZE
#undef EXP_SIZE
#undef EXP_MASK
#undef MANT_MASK
#undef SIGN_MASK
#undef IMANT_SIZE
#undef RND_SIZE
#undef QNAN_MASK
#undef F_QNAN

#undef pack_sf
#undef unpack_sf
#undef rshift_rnd
#undef round_pack_sf
#undef normalize_sf
#undef normalize2_sf
#undef issignan_sf
#undef isnan_sf
#undef add_sf
#undef mul_sf
#undef fma_sf
#undef div_sf
#undef sqrt_sf
#undef normalize_subnormal_sf
#undef divrem_u
#undef sqrtrem_u
#undef mul_u
#undef cvt_sf32_sf
#undef cvt_sf64_sf