//===-- lib/divsf3.c - Single-precision division ------------------*- C -*-===//

//

//                     The LLVM Compiler Infrastructure

//

// This file is dual licensed under the MIT and the University of Illinois Open

// Source Licenses. See LICENSE.TXT for details.

//

//===----------------------------------------------------------------------===//

//

// This file implements single-precision soft-float division

// with the IEEE-754 default rounding (to nearest, ties to even).

//

// For simplicity, this implementation currently flushes denormals to zero.

// It should be a fairly straightforward exercise to implement gradual

// underflow with correct rounding.

//

//===----------------------------------------------------------------------===//

#define SINGLE_PRECISION

#include "fp_lib.h"

ARM_EABI_FNALIAS(fdiv, divsf3)

COMPILER_RT_ABI fp_t

__divsf3(fp_t a, fp_t b) {

    int scale = 0;

    // Detect if a or b is zero, denormal, infinity, or NaN.

        // NaN / anything = qNaN

        // anything / NaN = qNaN

            // infinity / infinity = NaN

            // infinity / anything else = +/- infinity

        }

        // anything else / infinity = +/- 0

            // zero / zero = NaN

            // zero / anything else = +/- zero

        }

        // anything else / zero = +/- infinity

        // one or both of a or b is denormal, the other (if applicable) is a

        // normal number.  Renormalize one or both of a and b, and set scale to

        // include the necessary exponent adjustment.

    }

    // Or in the implicit significand bit.  (If we fell through from the

    // denormal path it was already set by normalize( ), but setting it twice

    // won't hurt anything.)

    // Align the significand of b as a Q31 fixed-point number in the range

    // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax

    // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2.  This

    // is accurate to about 3.5 binary digits.

    // Now refine the reciprocal estimate using a Newton-Raphson iteration:

    //

    //     x1 = x0 * (2 - x0 * b)

    //

    // This doubles the number of correct binary digits in the approximation

    // with each iteration, so after three iterations, we have about 28 binary

    // digits of accuracy.

    uint32_t correction;

    // Exhaustive testing shows that the error in reciprocal after three steps

    // is in the interval [-0x1.f58108p-31, 0x1.d0e48cp-29], in line with our

    // expectations.  We bump the reciprocal by a tiny value to force the error

    // to be strictly positive (in the range [0x1.4fdfp-37,0x1.287246p-29], to

    // be specific).  This also causes 1/1 to give a sensible approximation

    // instead of zero (due to overflow).

    // The numerical reciprocal is accurate to within 2^-28, lies in the

    // interval [0x1.000000eep-1, 0x1.fffffffcp-1], and is strictly smaller

    // than the true reciprocal of b.  Multiplying a by this reciprocal thus

    // gives a numerical q = a/b in Q24 with the following properties:

    //

    //    1. q < a/b

    //    2. q is in the interval [0x1.000000eep-1, 0x1.fffffffcp0)

    //    3. the error in q is at most 2^-24 + 2^-27 -- the 2^24 term comes

    //       from the fact that we truncate the product, and the 2^27 term

    //       is the error in the reciprocal of b scaled by the maximum

    //       possible value of a.  As a consequence of this error bound,

    //       either q or nextafter(q) is the correctly rounded

    // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).

    // In either case, we are going to compute a residual of the form

    //

    //     r = a - q*b

    //

    // We know from the construction of q that r satisfies:

    //

    //     0 <= r < ulp(q)*b

    //

    // if r is greater than 1/2 ulp(q)*b, then q rounds up.  Otherwise, we

    // already have the correct result.  The exact halfway case cannot occur.

    // We also take this time to right shift quotient if it falls in the [1,2)

    // range and adjust the exponent accordingly.

    rep_t residual;

    }

        // If we have overflowed the exponent, return infinity.

    }

        // Flush denormals to zero.  In the future, it would be nice to add

        // code to round them correctly.

    }

    else {

        // Clear the implicit bit

        // Insert the exponent

        // Round

        // Insert the sign and return

    }

}