2878 lines · 108.8 KB · Java
RawDownload
                    B2 += bigIntExp;                } else {                    D2 -= bigIntExp;                }                Ulp2 = B2;                // shift bigB and bigD left by a number s. t.                // halfUlp is still an integer.                int hulpbias;                if ( bigIntExp+bigIntNBits <= -expBias+1 ){                    // This is going to be a denormalized number                    // (if not actually zero).                    // half an ULP is at 2^-(expBias+expShift+1)                    hulpbias = bigIntExp+ expBias + expShift;                } else {                    hulpbias = expShift + 2 - bigIntNBits;                }                B2 += hulpbias;                D2 += hulpbias;                // if there are common factors of 2, we might just as well                // factor them out, as they add nothing useful.                int common2 = Math.min( B2, Math.min( D2, Ulp2 ) );                B2 -= common2;                D2 -= common2;                Ulp2 -= common2;                // do multiplications by powers of 5 and 2                bigB = multPow52( bigB, B5, B2 );                FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 );                //                // to recap:                // bigB is the scaled-big-int version of our floating-point                // candidate.                // bigD is the scaled-big-int version of the exact value                // as we understand it.                // halfUlp is 1/2 an ulp of bigB, except for special cases                // of exact powers of 2                //                // the plan is to compare bigB with bigD, and if the difference                // is less than halfUlp, then we're satisfied. Otherwise,                // use the ratio of difference to halfUlp to calculate a fudge                // factor to add to the floating value, then go 'round again.                //                FDBigInt diff;                int cmpResult;                boolean overvalue;                if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){                    overvalue = true; // our candidate is too big.                    diff = bigB.sub( bigD );                    if ( (bigIntNBits == 1) && (bigIntExp > -expBias+1) ){                        // candidate is a normalized exact power of 2 and                        // is too big. We will be subtracting.                        // For our purposes, ulp is the ulp of the                        // next smaller range.                        Ulp2 -= 1;                        if ( Ulp2 < 0 ){                            // rats. Cannot de-scale ulp this far.                            // must scale diff in other direction.                            Ulp2 = 0;                            diff.lshiftMe( 1 );                        }                    }                } else if ( cmpResult < 0 ){                    overvalue = false; // our candidate is too small.                    diff = bigD.sub( bigB );                } else {                    // the candidate is exactly right!                    // this happens with surprising frequency                    break correctionLoop;                }                FDBigInt halfUlp = constructPow52( B5, Ulp2 );                if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){                    // difference is small.                    // this is close enough                    if (mustSetRoundDir) {                        roundDir = overvalue ? -1 : 1;                    }                    break correctionLoop;                } else if ( cmpResult == 0 ){                    // difference is exactly half an ULP                    // round to some other value maybe, then finish                    dValue += 0.5*ulp( dValue, overvalue );                    // should check for bigIntNBits == 1 here??                    if (mustSetRoundDir) {                        roundDir = overvalue ? -1 : 1;                    }                    break correctionLoop;                } else {                    // difference is non-trivial.                    // could scale addend by ratio of difference to                    // halfUlp here, if we bothered to compute that difference.                    // Most of the time ( I hope ) it is about 1 anyway.                    dValue += ulp( dValue, overvalue );                    if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY )                        break correctionLoop; // oops. Fell off end of range.                    continue; // try again.                }             }            return (isNegative)? -dValue : dValue;        }    }     /*     * Take a FloatingDecimal, which we presumably just scanned in,     * and find out what its value is, as a float.     * This is distinct from doubleValue() to avoid the extremely     * unlikely case of a double rounding error, wherein the conversion     * to double has one rounding error, and the conversion of that double     * to a float has another rounding error, IN THE WRONG DIRECTION,     * ( because of the preference to a zero low-order bit ).     */     public strictfp float floatValue(){        int     kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 );        int     iValue;        float   fValue;         // First, check for NaN and Infinity values        if(digits == infinity || digits == notANumber) {            if(digits == notANumber)                return Float.NaN;            else                return (isNegative?Float.NEGATIVE_INFINITY:Float.POSITIVE_INFINITY);        }        else {            /*             * convert the lead kDigits to an integer.             */            iValue = (int)digits[0]-(int)'0';            for ( int i=1; i < kDigits; i++ ){                iValue = iValue*10 + (int)digits[i]-(int)'0';            }            fValue = (float)iValue;            int exp = decExponent-kDigits;            /*             * iValue now contains an integer with the value of             * the first kDigits digits of the number.             * fValue contains the (float) of the same.             */             if ( nDigits <= singleMaxDecimalDigits ){                /*                 * possibly an easy case.                 * We know that the digits can be represented                 * exactly. And if the exponent isn't too outrageous,                 * the whole thing can be done with one operation,                 * thus one rounding error.                 * Note that all our constructors trim all leading and                 * trailing zeros, so simple values (including zero)                 * will always end up here.                 */                if (exp == 0 || fValue == 0.0f)                    return (isNegative)? -fValue : fValue; // small floating integer                else if ( exp >= 0 ){                    if ( exp <= singleMaxSmallTen ){                        /*                         * Can get the answer with one operation,                         * thus one roundoff.                         */                        fValue *= singleSmall10pow[exp];                        return (isNegative)? -fValue : fValue;                    }                    int slop = singleMaxDecimalDigits - kDigits;                    if ( exp <= singleMaxSmallTen+slop ){                        /*                         * We can multiply dValue by 10^(slop)                         * and it is still "small" and exact.                         * Then we can multiply by 10^(exp-slop)                         * with one rounding.                         */                        fValue *= singleSmall10pow[slop];                        fValue *= singleSmall10pow[exp-slop];                        return (isNegative)? -fValue : fValue;                    }                    /*                     * Else we have a hard case with a positive exp.                     */                } else {                    if ( exp >= -singleMaxSmallTen ){                        /*                         * Can get the answer in one division.                         */                        fValue /= singleSmall10pow[-exp];                        return (isNegative)? -fValue : fValue;                    }                    /*                     * Else we have a hard case with a negative exp.                     */                }            } else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){                /*                 * In double-precision, this is an exact floating integer.                 * So we can compute to double, then shorten to float                 * with one round, and get the right answer.                 *                 * First, finish accumulating digits.                 * Then convert that integer to a double, multiply                 * by the appropriate power of ten, and convert to float.                 */                long lValue = (long)iValue;                for ( int i=kDigits; i < nDigits; i++ ){                    lValue = lValue*10L + (long)((int)digits[i]-(int)'0');                }                double dValue = (double)lValue;                exp = decExponent-nDigits;                dValue *= small10pow[exp];                fValue = (float)dValue;                return (isNegative)? -fValue : fValue;             }            /*             * Harder cases:             * The sum of digits plus exponent is greater than             * what we think we can do with one error.             *             * Start by weeding out obviously out-of-range             * results, then convert to double and go to             * common hard-case code.             */            if ( decExponent > singleMaxDecimalExponent+1 ){                /*                 * Lets face it. This is going to be                 * Infinity. Cut to the chase.                 */                return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;            } else if ( decExponent < singleMinDecimalExponent-1 ){                /*                 * Lets face it. This is going to be                 * zero. Cut to the chase.                 */                return (isNegative)? -0.0f : 0.0f;            }             /*             * Here, we do 'way too much work, but throwing away             * our partial results, and going and doing the whole             * thing as double, then throwing away half the bits that computes             * when we convert back to float.             *             * The alternative is to reproduce the whole multiple-precision             * algorithm for float precision, or to try to parameterize it             * for common usage. The former will take about 400 lines of code,             * and the latter I tried without success. Thus the semi-hack             * answer here.             */            mustSetRoundDir = !fromHex;            double dValue = doubleValue();            return stickyRound( dValue );        }    }      /*     * All the positive powers of 10 that can be     * represented exactly in double/float.     */    private static final double small10pow[] = {        1.0e0,        1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5,        1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10,        1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15,        1.0e16, 1.0e17, 1.0e18, 1.0e19, 1.0e20,        1.0e21, 1.0e22    };     private static final float singleSmall10pow[] = {        1.0e0f,        1.0e1f, 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f,        1.0e6f, 1.0e7f, 1.0e8f, 1.0e9f, 1.0e10f    };     private static final double big10pow[] = {        1e16, 1e32, 1e64, 1e128, 1e256 };    private static final double tiny10pow[] = {        1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };     private static final int maxSmallTen = small10pow.length-1;    private static final int singleMaxSmallTen = singleSmall10pow.length-1;     private static final int small5pow[] = {        1,        5,        5*5,        5*5*5,        5*5*5*5,        5*5*5*5*5,        5*5*5*5*5*5,        5*5*5*5*5*5*5,        5*5*5*5*5*5*5*5,        5*5*5*5*5*5*5*5*5,        5*5*5*5*5*5*5*5*5*5,        5*5*5*5*5*5*5*5*5*5*5,        5*5*5*5*5*5*5*5*5*5*5*5,        5*5*5*5*5*5*5*5*5*5*5*5*5    };      private static final long long5pow[] = {        1L,        5L,        5L*5,        5L*5*5,        5L*5*5*5,        5L*5*5*5*5,        5L*5*5*5*5*5,        5L*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,    };     // approximately ceil( log2( long5pow[i] ) )    private static final int n5bits[] = {        0,        3,        5,        7,        10,        12,        14,        17,        19,        21,        24,        26,        28,        31,        33,        35,        38,        40,        42,        45,        47,        49,        52,        54,        56,        59,        61,    };     private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' };    private static final char notANumber[] = { 'N', 'a', 'N' };    private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' };      /*     * Grammar is compatible with hexadecimal floating-point constants     * described in section 6.4.4.2 of the C99 specification.     */    private static Pattern hexFloatPattern = null;    private static synchronized Pattern getHexFloatPattern() {        if (hexFloatPattern == null) {           hexFloatPattern = Pattern.compile(                   //1           234                   56                7                   8      9                    "([-+])?0[xX](((\\p{XDigit}+)\\.?)|((\\p{XDigit}*)\\.(\\p{XDigit}+)))[pP]([-+])?(\\p{Digit}+)[fFdD]?"                    );        }        return hexFloatPattern;    }     /*     * Convert string s to a suitable floating decimal; uses the     * double constructor and set the roundDir variable appropriately     * in case the value is later converted to a float.     */   static FloatingDecimal parseHexString(String s) {        // Verify string is a member of the hexadecimal floating-point        // string language.        Matcher m = getHexFloatPattern().matcher(s);        boolean validInput = m.matches();         if (!validInput) {            // Input does not match pattern            throw new NumberFormatException("For input string: \"" + s + "\"");        } else { // validInput            /*             * We must isolate the sign, significand, and exponent             * fields.  The sign value is straightforward.  Since             * floating-point numbers are stored with a normalized             * representation, the significand and exponent are             * interrelated.             *             * After extracting the sign, we normalized the             * significand as a hexadecimal value, calculating an             * exponent adjust for any shifts made during             * normalization.  If the significand is zero, the             * exponent doesn't need to be examined since the output             * will be zero.             *             * Next the exponent in the input string is extracted.             * Afterwards, the significand is normalized as a *binary*             * value and the input value's normalized exponent can be             * computed.  The significand bits are copied into a             * double significand; if the string has more logical bits             * than can fit in a double, the extra bits affect the             * round and sticky bits which are used to round the final             * value.             */             //  Extract significand sign            String group1 = m.group(1);            double sign = (( group1 == null ) || group1.equals("+"))? 1.0 : -1.0;              //  Extract Significand magnitude            /*             * Based on the form of the significand, calculate how the             * binary exponent needs to be adjusted to create a             * normalized *hexadecimal* floating-point number; that             * is, a number where there is one nonzero hex digit to             * the left of the (hexa)decimal point.  Since we are             * adjusting a binary, not hexadecimal exponent, the             * exponent is adjusted by a multiple of 4.             *             * There are a number of significand scenarios to consider;             * letters are used in indicate nonzero digits:             *             * 1. 000xxxx       =>      x.xxx   normalized             *    increase exponent by (number of x's - 1)*4             *             * 2. 000xxx.yyyy =>        x.xxyyyy        normalized             *    increase exponent by (number of x's - 1)*4             *             * 3. .000yyy  =>   y.yy    normalized             *    decrease exponent by (number of zeros + 1)*4             *             * 4. 000.00000yyy => y.yy normalized             *    decrease exponent by (number of zeros to right of point + 1)*4             *             * If the significand is exactly zero, return a properly             * signed zero.             */             String significandString =null;            int signifLength = 0;            int exponentAdjust = 0;            {                int leftDigits  = 0; // number of meaningful digits to                                     // left of "decimal" point                                     // (leading zeros stripped)                int rightDigits = 0; // number of digits to right of                                     // "decimal" point; leading zeros                                     // must always be accounted for                /*                 * The significand is made up of either                 *                 * 1. group 4 entirely (integer portion only)                 *                 * OR                 *                 * 2. the fractional portion from group 7 plus any                 * (optional) integer portions from group 6.                 */                String group4;                if( (group4 = m.group(4)) != null) {  // Integer-only significand                    // Leading zeros never matter on the integer portion                    significandString = stripLeadingZeros(group4);                    leftDigits = significandString.length();                }                else {                    // Group 6 is the optional integer; leading zeros                    // never matter on the integer portion                    String group6 = stripLeadingZeros(m.group(6));                    leftDigits = group6.length();                     // fraction                    String group7 = m.group(7);                    rightDigits = group7.length();                     // Turn "integer.fraction" into "integer"+"fraction"                    significandString =                        ((group6 == null)?"":group6) + // is the null                        // check necessary?                        group7;                }                 significandString = stripLeadingZeros(significandString);                signifLength  = significandString.length();                 /*                 * Adjust exponent as described above                 */                if (leftDigits >= 1) {  // Cases 1 and 2

Lines 1501–2000 of 2878. Previous lines