Tcl Source Code

Check-in [50ecdc3e87]
Login
EuroTcl/OpenACS 11 - 12 JULY 2024, VIENNA

Many hyperlinks are disabled.
Use anonymous login to enable hyperlinks.

Overview
Comment:Restore FP control word on conversion of zero values.

Sneak path existed that failed to restore the floating point control word when scanning a zero value. The result was that floating point rounding was incorrectly left set to 53-bit significance, round-to-even, which is incompatible with the math library from musl.

Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | bug-713653b951
Files: files | file ages | folders
SHA3-256: 50ecdc3e87fb5d772490b5490090d3ce706d9ceeb6b9725f1e61b049dee2fd5a
User & Date: kbk 2022-07-15 17:11:41
Context
2022-07-17
12:57
Fix [713653b951]: Floating point precision problems on x86 musl check-in: 5e9b6b560e user: jan.nijtmans tags: core-8-5-branch
2022-07-15
17:11
Restore FP control word on conversion of zero values.

Sneak path existed that failed to restore the... Closed-Leaf check-in: 50ecdc3e87 user: kbk tags: bug-713653b951

2022-07-12
08:20
Code cleanup (use {} in if/else statemenets) check-in: ba232277ee user: jan.nijtmans tags: core-8-6-branch
Changes
Hide Diffs Unified Diffs Ignore Whitespace Patch

Changes to generic/tclStrToD.c.

1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638

1639
1640
1641
1642
1643
1644

1645
1646
1647
1648
1649







1650
1651
1652
1653
1654
1655
1656
static double
MakeLowPrecisionDouble(
    int signum,			/* 1 if the number is negative, 0 otherwise */
    Tcl_WideUInt significand,	/* Significand of the number */
    int numSigDigs,		/* Number of digits in the significand */
    long exponent)		/* Power of ten */
{
    double retval;		/* Value of the number. */
    mp_int significandBig;	/* Significand expressed as a bignum. */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits.

     */

    TCL_IEEE_DOUBLE_ROUNDING;

    /*
     * Test for the easy cases.

     */

    if (significand == 0) {
	return copysign(0.0, -signum);
    }







    if (numSigDigs <= QUICK_MAX) {
	if (exponent >= 0) {
	    if (exponent <= mmaxpow) {
		/*
		 * The significand is an exact integer, and so is
		 * 10**exponent. The product will be correct to within 1/2 ulp
		 * without special handling.







<







|
>

|
<


|
>

<



>
>
>
>
>
>
>







1623
1624
1625
1626
1627
1628
1629

1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640

1641
1642
1643
1644
1645

1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
static double
MakeLowPrecisionDouble(
    int signum,			/* 1 if the number is negative, 0 otherwise */
    Tcl_WideUInt significand,	/* Significand of the number */
    int numSigDigs,		/* Number of digits in the significand */
    long exponent)		/* Power of ten */
{

    mp_int significandBig;	/* Significand expressed as a bignum. */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits. We also make
     * 'retval' volatile, so that it doesn't get promoted to a register.
     */
    volatile double retval;		/* Value of the number. */


    /*
     * Test for zero significand, which requires explicit construction
     * of -0.0. (Unary minus returns a positive zero.)
     */

    if (significand == 0) {
	return copysign(0.0, -signum);
    }

    /*
     * Set the FP control word for 53 bits, WARNING: It must be reset
     * before returning.
     */
    TCL_IEEE_DOUBLE_ROUNDING;

    if (numSigDigs <= QUICK_MAX) {
	if (exponent >= 0) {
	    if (exponent <= mmaxpow) {
		/*
		 * The significand is an exact integer, and so is
		 * 10**exponent. The product will be correct to within 1/2 ulp
		 * without special handling.
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755


1756

1757












1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
static double
MakeHighPrecisionDouble(
    int signum,			/* 1=negative, 0=nonnegative */
    mp_int *significand,	/* Exact significand of the number */
    int numSigDigs,		/* Number of significant digits */
    long exponent)		/* Power of 10 by which to multiply */
{
    double retval;
    int machexp;		/* Machine exponent of a power of 10. */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits.


     */














    TCL_IEEE_DOUBLE_ROUNDING;

    /*
     * Quick checks for zero, and over/underflow. Be careful to avoid
     * integer overflow when calculating with 'exponent'.
     */

    if (mp_iszero(significand)) {
	return copysign(0.0, -signum);
    }
    if (exponent >= 0 && exponent-1 > maxDigits-numSigDigs) {
	retval = HUGE_VAL;
	goto returnValue;
    } else if (exponent < 0 && numSigDigs+exponent < minDigits+1) {
	retval = 0.0;
	goto returnValue;
    }







<







|
>
>

>

>
>
>
>
>
>
>
>
>
>
>
>



|


<
<
<
<







1746
1747
1748
1749
1750
1751
1752

1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783




1784
1785
1786
1787
1788
1789
1790
static double
MakeHighPrecisionDouble(
    int signum,			/* 1=negative, 0=nonnegative */
    mp_int *significand,	/* Exact significand of the number */
    int numSigDigs,		/* Number of significant digits */
    long exponent)		/* Power of 10 by which to multiply */
{

    int machexp;		/* Machine exponent of a power of 10. */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits. We also make
     * 'retval' volatile to make sure that it doesn't get promoted to a
     * register.
     */
    volatile double retval;

    /* 
     * A zero significand requires explicit construction of -0.0.
     * (Unary minus returns positive zero.)
     */
    if (mp_iszero(significand)) {
	return copysign(0.0, -signum);
    }

    /*
     * Set the 53-bit rounding mode. WARNING: It must be reset before
     * returning.
     */
    TCL_IEEE_DOUBLE_ROUNDING;

    /*
     * Make quick checks for over/underflow. Be careful to avoid
     * integer overflow when calculating with 'exponent'.
     */




    if (exponent >= 0 && exponent-1 > maxDigits-numSigDigs) {
	retval = HUGE_VAL;
	goto returnValue;
    } else if (exponent < 0 && numSigDigs+exponent < minDigits+1) {
	retval = 0.0;
	goto returnValue;
    }