math/single/rredf.c

/*
 * rredf.c - trigonometric range reduction function
 *
 * Copyright (c) 2009-2018, Arm Limited.
 * SPDX-License-Identifier: MIT
 */

/*
 * This code is intended to be used as the second half of a range
 * reducer whose first half is an inline function defined in
 * rredf.h. Each trig function performs range reduction by invoking
 * that, which handles the quickest and most common cases inline
 * before handing off to this function for everything else. Thus a
 * reasonable compromise is struck between speed and space. (I
 * hope.) In particular, this approach avoids a function call
 * overhead in the common case.
 */

#include "math_private.h"

#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */

/*
 * Input values to this function:
 *  - x is the original user input value, unchanged by the
 *    first-tier reducer in the case where it hands over to us.
 *  - q is still the place where the caller expects us to leave the
 *    quadrant code.
 *  - k is the IEEE bit pattern of x (which it would seem a shame to
 *    recompute given that the first-tier reducer already went to
 *    the effort of extracting it from the VFP). FIXME: in softfp,
 *    on the other hand, it's unconscionably wasteful to replicate
 *    this value into a second register and we should change the
 *    prototype!
 */
float __mathlib_rredf2(float x, int *q, unsigned k)
{
    /*
     * First, weed out infinities and NaNs, and deal with them by
     * returning a negative q.
     */
    if ((k << 1) >= 0xFF000000) {
        *q = -1;
        return x;
    }
    /*
     * We do general range reduction by multiplying by 2/pi, and
     * retaining the bottom two bits of the integer part and an
     * initial chunk of the fraction below that. The integer bits
     * are directly output as *q; the fraction is then multiplied
     * back up by pi/2 before returning it.
     *
     * To get this right, we don't have to multiply by the _whole_
     * of 2/pi right from the most significant bit downwards:
     * instead we can discard any bit of 2/pi with a place value
     * high enough that multiplying it by the LSB of x will yield a
     * place value higher than 2. Thus we can bound the required
     * work by a reasonably small constant regardless of the size of
     * x (unlike, for instance, the IEEE remainder operation).
     *
     * At the other end, however, we must take more care: it isn't
     * adequate just to acquire two integer bits and 24 fraction
     * bits of (2/pi)x, because if a lot of those fraction bits are
     * zero then we will suffer significance loss. So we must keep
     * computing fraction bits as far down as 23 bits below the
     * _highest set fraction bit_.
     *
     * The immediate question, therefore, is what the bound on this
     * end of the job will be. In other words: what is the smallest
     * difference between an integer multiple of pi/2 and a
     * representable IEEE single precision number larger than the
     * maximum size handled by rredf.h?
     *
     * The most difficult cases for each exponent can readily be
     * found by Tim Peters's modular minimisation algorithm, and are
     * tabulated in mathlib/tests/directed/rredf.tst. The single
     * worst case is the IEEE single-precision number 0x6F79BE45,
     * whose numerical value is in the region of 7.7*10^28; when
     * reduced mod pi/2, it attains the value 0x30DDEEA9, or about
     * 0.00000000161. The highest set bit of this value is the one
     * with place value 2^-30; so its lowest is 2^-53. Hence, to be
     * sure of having enough fraction bits to output at full single
     * precision, we must be prepared to collect up to 53 bits of
     * fraction in addition to our two bits of integer part.
     *
     * To begin with, this means we must store the value of 2/pi to
     * a precision of 128+53 = 181 bits. That's six 32-bit words.
     * (Hardly a chore, unlike the equivalent problem in double
     * precision!)
     */
    {
        static const unsigned twooverpi[] = {
            /* We start with a zero word, because that takes up less
             * space than the array bounds checking and special-case
             * handling that would have to occur in its absence. */
            0,
            /* 2/pi in hex is 0.a2f9836e... */
            0xa2f9836e, 0x4e441529, 0xfc2757d1,
            0xf534ddc0, 0xdb629599, 0x3c439041,
            /* Again, to avoid array bounds overrun, we store a spare
             * word at the end. And it would be a shame to fill it
             * with zeroes when we could use more bits of 2/pi... */
            0xfe5163ab
        };

        /*
         * Multiprecision multiplication of this nature is more
         * readily done in integers than in VFP, since we can use
         * UMULL (on CPUs that support it) to multiply 32 by 32 bits
         * at a time whereas the VFP would only be able to do 12x12
         * without losing accuracy.
         *
         * So extract the mantissa of the input number as a 32-bit
         * integer.
         */
        unsigned mantissa = 0x80000000 | (k << 8);

        /*
         * Now work out which part of our stored value of 2/pi we're
         * supposed to be multiplying by.
         *
         * Let the IEEE exponent field of x be e. With its bias
         * removed, (e-127) is the index of the set bit at the top
         * of 'mantissa' (i.e. that set bit has real place value
         * 2^(e-127)). So the lowest set bit in 'mantissa', 23 bits
         * further down, must have place value 2^(e-150).
         *
         * We begin taking an interest in the value of 2/pi at the
         * bit which multiplies by _that_ to give something with
         * place value at most 2. In other words, the highest bit of
         * 2/pi we're interested in is the one with place value
         * 2/(2^(e-150)) = 2^(151-e).
         *
         * The bit at the top of the first (zero) word of the above
         * array has place value 2^31. Hence, the bit we want to put
         * at the top of the first word we extract from that array
         * is the one at bit index n, where 31-n = 151-e and hence
         * n=e-120.
         */
        int topbitindex = ((k >> 23) & 0xFF) - 120;
        int wordindex = topbitindex >> 5;
        int shiftup = topbitindex & 31;
        int shiftdown = 32 - shiftup;
        unsigned word1, word2, word3;
        if (shiftup) {
            word1 = (twooverpi[wordindex] << shiftup) | (twooverpi[wordindex+1] >> shiftdown);
            word2 = (twooverpi[wordindex+1] << shiftup) | (twooverpi[wordindex+2] >> shiftdown);
            word3 = (twooverpi[wordindex+2] << shiftup) | (twooverpi[wordindex+3] >> shiftdown);
        } else {
            word1 = twooverpi[wordindex];
            word2 = twooverpi[wordindex+1];
            word3 = twooverpi[wordindex+2];
        }

        /*
         * Do the multiplications, and add them together.
         */
        unsigned long long mult1 = (unsigned long long)word1 * mantissa;
        unsigned long long mult2 = (unsigned long long)word2 * mantissa;
        unsigned long long mult3 = (unsigned long long)word3 * mantissa;

        unsigned /* bottom3 = (unsigned)mult3, */ top3 = (unsigned)(mult3 >> 32);
        unsigned bottom2 = (unsigned)mult2, top2 = (unsigned)(mult2 >> 32);
        unsigned bottom1 = (unsigned)mult1, top1 = (unsigned)(mult1 >> 32);

        unsigned out3, out2, out1, carry;

        out3 = top3 + bottom2; carry = (out3 < top3);
        out2 = top2 + bottom1 + carry; carry = carry ? (out2 <= top2) : (out2 < top2);
        out1 = top1 + carry;

        /*
         * The two words we multiplied to get mult1 had their top
         * bits at (respectively) place values 2^(151-e) and
         * 2^(e-127). The value of those two bits multiplied
         * together will have ended up in bit 62 (the
         * topmost-but-one bit) of mult1, i.e. bit 30 of out1.
         * Hence, that bit has place value 2^(151-e+e-127) = 2^24.
         * So the integer value that we want to output as q,
         * consisting of the bits with place values 2^1 and 2^0,
         * must be 23 and 24 bits below that, i.e. in bits 7 and 6
         * of out1.
         *
         * Or, at least, it will be once we add 1/2, to round to the
         * _nearest_ multiple of pi/2 rather than the next one down.
         */
        *q = (out1 + (1<<5)) >> 6;

        /*
         * Now we construct the output fraction, which is most
         * simply done in the VFP. We just extract three consecutive
         * bit strings from our chunk of binary data, convert them
         * to integers, equip each with an appropriate FP exponent,
         * add them together, and (don't forget) multiply back up by
         * pi/2. That way we don't have to work out ourselves where
         * the highest fraction bit ended up.
         *
         * Since our displacement from the nearest multiple of pi/2
         * can be positive or negative, the topmost of these three
         * values must be arranged with its 2^-1 bit at the very top
         * of the word, and then treated as a _signed_ integer.
         */
        {
            int i1 = (out1 << 26) | ((out2 >> 19) << 13);
            unsigned i2 = out2 << 13;
            unsigned i3 = out3;
            float f1 = i1, f2 = i2 * (1.0f/524288.0f), f3 = i3 * (1.0f/524288.0f/524288.0f);

            /*
             * Now f1+f2+f3 is a representation, potentially to
             * twice double precision, of 2^32 times ((2/pi)*x minus
             * some integer). So our remaining job is to multiply
             * back down by (pi/2)*2^-32, and convert back to one
             * single-precision output number.
             */

            /* Normalise to a prec-and-a-half representation... */
            float ftop = CLEARBOTTOMHALF(f1+f2+f3), fbot = f3-((ftop-f1)-f2);

            /* ... and multiply by a prec-and-a-half value of (pi/2)*2^-32. */
            float ret = (ftop * 0x1.92p-32F) + (ftop * 0x1.fb5444p-44F + fbot * 0x1.921fb6p-32F);

            /* Just before we return, take the input sign into account. */
            if (k & 0x80000000) {
                *q = 0x10000000 - *q;
                ret = -ret;
            }
            return ret;
        }
    }
}

#ifdef __cplusplus
} /* end of extern "C" */
#endif /* __cplusplus */

/* end of rredf.c */