//===-- lib/divdf3.c - Double-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 double-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 DOUBLE_PRECISION

#include "fp_lib.h"

ARM_EABI_FNALIAS(ddiv, divdf3)

COMPILER_RT_ABI fp_t

__divdf3(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 correction32;

    // recip32 might have overflowed to exactly zero in the preceding

    // computation if the high word of b is exactly 1.0.  This would sabotage

    // the full-width final stage of the computation that follows, so we adjust

    // recip32 downward by one bit.

    // We need to perform one more iteration to get us to 56 binary digits;

    // The last iteration needs to happen with extra precision.

    uint64_t correction, reciprocal;

    // We already adjusted the 32-bit estimate, now we need to adjust the final

    // 64-bit reciprocal estimate downward to ensure that it is strictly smaller

    // than the infinitely precise exact reciprocal.  Because the computation

    // of the Newton-Raphson step is truncating at every step, this adjustment

    // is small; most of the work is already done.

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

    // interval [0.5, 1.0), and is strictly smaller than the true reciprocal

    // of b.  Multiplying a by this reciprocal thus gives a numerical q = a/b

    // in Q53 with the following properties:

    //

    //    1. q < a/b

    //    2. q is in the interval [0.5, 2.0)

    //    3. the error in q is bounded away from 2^-53 (actually, we have a

    //       couple of bits to spare, but this is all we need).

    // We need a 64 x 64 multiply high to compute q, which isn't a basic

    // operation in C, so we need to be a little bit fussy.

    // 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

        return result;

    }

}