Studying note of GCC-3.4.6 source (79)

5.6.1.1.2.2. Case of floating point number
Analysizing integer literal string is quite trouble, but floating point literal string is a much complexer case. First, it needs select the type node associating with the building REAL_CST node. By default, the floating point has type double (CPP_N_MEDIUM).

568 static tree
569 interpret_float (const cpp_token *token, unsigned int flags) in c-lex.c
570 {
571 tree type;
572 tree value;
573 REAL_VALUE_TYPE real;
574 char *copy;
575 size_t copylen;
576 const char *typename;
577
578 /* FIXME: make %T work in error/warning, then we don't need typename. */
579 if ((flags & CPP_N_WIDTH) == CPP_N_LARGE)
580 {
581 type = long_double_type_node;
582 typename = "long double";
583 }
584 else if ((flags & CPP_N_WIDTH) == CPP_N_SMALL
585 || flag_single_precision_constant)
586 {
587 type = float_type_node;
588 typename = "float";
589 }
590 else
591 {
592 type = double_type_node;
593 typename = "double";
594 }
595
596 /* Copy the constant to a nul-terminated buffer. If the constant
597 has any suffixes, cut them off; REAL_VALUE_ATOF/ REAL_VALUE_HTOF
598 can't handle them. */
599 copylen = token->val.str.len;
600 if ((flags & CPP_N_WIDTH) != CPP_N_MEDIUM)
601 /* Must be an F or L suffix. */
602 copylen--;
603 if (flags & CPP_N_IMAGINARY)
604 /* I or J suffix. */
605 copylen--;
606
607 copy = alloca (copylen + 1);
608 memcpy (copy, token->val.str.text, copylen);
609 copy[copylen] = '/0';
610
611 real_from_string (&real, copy);
612 real_convert (&real, TYPE_MODE (type), &real);
613
614 /* A diagnostic is required for "soft" overflow by some ISO C
615 testsuites. This is not pedwarn, because some people don't want
616 an error for this.
617 ??? That's a dubious reason... is this a mandatory diagnostic or
618 isn't it? -- zw, 2001-08-21. */
619 if (REAL_VALUE_ISINF (real) && pedantic)
620 warning ("floating constant exceeds range of /"%s/"", typename);
621
622 /* Create a node with determined type and value. */
623 value = build_real (type, real);
624 if (flags & CPP_N_IMAGINARY)
625 value = build_complex (NULL_TREE, convert (type, integer_zero_node), value);
626
627 return value;
628 }

Next, it decides the value of the lieral string. There are 2 representations in floating point literal string, one is Hex which has exponent part power of 2; the other is Decimal which has exponent part power of 10. No matter which kind, at line 1765, first clear r and set it as unsigned.

1759 void
1760 real_from_string (REAL_VALUE_TYPE *r, const char *str) in real.c
1761{
1762 int exp = 0;
1763 bool sign = false;
1764
1765 get_zero (r, 0);
1766
1767 if (*str == '-')
1768 {
1769 sign = true;
1770 str++;
1771 }
1772 else if (*str == '+')
1773 str++;
1774
1775 if (str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
1776 {
1777 /* Hexadecimal floating point. */
1778 int pos = SIGNIFICAND_BITS - 4, d;
1779
1780 str += 2;
1781
1782 while (*str == '0')
1783 str++;
1784 while (1)
1785 {
1786 d = hex_value (*str);
1787 if (d == _hex_bad)
1788 break;
1789 if (pos >= 0)
1790 {
1791 r->sig[pos / HOST_BITS_PER_LONG]
1792 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1793 pos -= 4;
1794 }
1795 exp += 4;
1796 str++;
1797 }
1798 if (*str == '.')
1799 {
1800 str++;
1801 if (pos == SIGNIFICAND_BITS - 4)
1802 {
1803 while (*str == '0')
1804 str++, exp -= 4;
1805 }
1806 while (1)
1807 {
1808 d = hex_value (*str);
1809 if (d == _hex_bad)
1810 break;
1811 if (pos >= 0)
1812 {
1813 r->sig[pos / HOST_BITS_PER_LONG]
1814 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1815 pos -= 4;
1816 }
1817 str++;
1818 }
1819 }
1820 if (*str == 'p' || *str == 'P')
1821 {
1822 bool exp_neg = false;
1823
1824 str++;
1825 if (*str == '-')
1826 {
1827 exp_neg = true;
1828 str++;
1829 }
1830 else if (*str == '+')
1831 str++;
1832
1833 d = 0;
1834 while (ISDIGIT (*str))
1835 {
1836 d *= 10;
1837 d += *str - '0';
1838 if (d > MAX_EXP)
1839 {
1840 /* Overflowed the exponent. */
1841 if (exp_neg)
1842 goto underflow;
1843 else
1844 goto overflow;
1845 }
1846 str++;
1847 }
1848 if (exp_neg)
1849 d = -d;
1850
1851 exp += d;
1852 }
1853
1854 r->class = rvc_normal;
1855 r->exp = exp;
1856
1857 normalize (r);
1858 }

For the first form, it begins with “0x/0X”. As we have seen in before (section Node of tree_real_cst), REAL_VALUE_TYPE is the type compiler uses for floating point intermally, in which SIGNIFICAND_BITS is the significand bits in used (it is 128 + HOST_BITS_PER_LONG = 160) and EXP_BITS is the bits of exponent (27). Above code reveals if the literal string is too long, after filling out the significand part, the compiler will ignore it, but increate exponent part (if no exponent part in the literal string, it will be leaved to normalize for checking). Pay attention to the adjustment of exp at line 1795 and 1804, which keeps the significand part in form of 0.xx…x. At line 1857, normalize normalizes the result to achieve as high precision as possible (it is almost nothing to do here).
And for the Deciaml form, can’t simply use shift to achieve carry, the method is a little bit trouble.

real_from_string (continue)

1859 else
1860 {
1861 /* Decimal floating point. */
1862 const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
1863 int d;
1864
1865 while (*str == '0')
1866 str++;
1867 while (ISDIGIT (*str))
1868 {
1869 d = *str++ - '0';
1870 do_multiply (r, r, ten);
1871 if (d)
1872 do_add (r, r, real_digit (d), 0);
1873 }
1874 if (*str == '.')
1875 {
1876 str++;
1877 if (r->class == rvc_zero)
1878 {
1879 while (*str == '0')
1880 str++, exp--;
1881 }
1882 while (ISDIGIT (*str))
1883 {
1884 d = *str++ - '0';
1885 do_multiply (r, r, ten);
1886 if (d)
1887 do_add (r, r, real_digit (d), 0);
1888 exp--;
1889 }
1890 }
1891
1892 if (*str == 'e' || *str == 'E')
1893 {
1894 bool exp_neg = false;
1895
1896 str++;
1897 if (*str == '-')
1898 {
1899 exp_neg = true;
1900 str++;
1901 }
1902 else if (*str == '+')
1903 str++;
1904
1905 d = 0;
1906 while (ISDIGIT (*str))
1907 {
1908 d *= 10;
1909 d += *str - '0';
1910 if (d > MAX_EXP)
1911 {
1912 /* Overflowed the exponent. */
1913 if (exp_neg)
1914 goto underflow;
1915 else
1916 goto overflow;
1917 }
1918 str++;
1919 }
1920 if (exp_neg)
1921 d = -d;
1922 exp += d;
1923 }
1924
1925 if (exp)
1926 times_pten (r, exp);
1927 }
1928
1929 r->sign = sign;
1930 return;
1931
1932 underflow:
1933 get_zero (r, sign);
1934 return;
1935
1936 overflow:
1937 get_inf (r, sign);
1938 return;
1939 }

Function ten_to_ptwo returns 10**2**N in form of REAL_VALUE_TYPE. At line 1862, we get 10. See that at line 2012, rvc_zero indicates the value of tens is not available, or it should be rvc_normal.

2004 static const REAL_VALUE_TYPE * in real.c
2005 ten_to_ptwo (int n)
2006 {
2007 static REAL_VALUE_TYPE tens[EXP_BITS];
2008
2009 if (n < 0 || n >= EXP_BITS)
2010 abort ();
2011
2012 if (tens
.class == rvc_zero)
2013 {
2014 if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2015 {
2016 HOST_WIDE_INT t = 10;
2017 int i;
2018
2019 for (i = 0; i < n; ++i)
2020 t *= t;
2021
2022 real_from_integer (&tens
, VOIDmode, t, 0, 1);
2023 }
2024 else
2025 {
2026 const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2027 do_multiply (&tens
, t, t);
2028 }
2029 }
2030
2031 return &tens
;
2032 }

Function do_multiply only accepts arguemnts type of REAL_VALUE_TYPE. That is why uses ten_to_ptwo to construct REAL_VALUE_TYPE of 10. do_multiply performs r = a * b. First of all, it checks the special cases. Below CLASS2 has definition: #define CLASS2(A, B) ((A) << 2 | (B)), which returns unique value for pairs.

662 static bool
663 do_multiply (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a, in real.c
664 const REAL_VALUE_TYPE *b)
665 {
666 REAL_VALUE_TYPE u, t, *rr;
667 unsigned int i, j, k;
668 int sign = a->sign ^ b->sign;
669 bool inexact = false;
670
671 switch (CLASS2 (a->class, b->class))
672 {
673 case CLASS2 (rvc_zero, rvc_zero):
674 case CLASS2 (rvc_zero, rvc_normal):
675 case CLASS2 (rvc_normal, rvc_zero):
676 /* +-0 * ANY = 0 with appropriate sign. */
677 get_zero (r, sign);
678 return false;
679
680 case CLASS2 (rvc_zero, rvc_nan):
681 case CLASS2 (rvc_normal, rvc_nan):
682 case CLASS2 (rvc_inf, rvc_nan):
683 case CLASS2 (rvc_nan, rvc_nan):
684 /* ANY * NaN = NaN. */
685 *r = *b;
686 r->sign = sign;
687 return false;
688
689 case CLASS2 (rvc_nan, rvc_zero):
690 case CLASS2 (rvc_nan, rvc_normal):
691 case CLASS2 (rvc_nan, rvc_inf):
692 /* NaN * ANY = NaN. */
693 *r = *a;
694 r->sign = sign;
695 return false;
696
697 case CLASS2 (rvc_zero, rvc_inf):
698 case CLASS2 (rvc_inf, rvc_zero):
699 /* 0 * Inf = NaN */
700 get_canonical_qnan (r, sign);
701 return false;
702
703 case CLASS2 (rvc_inf, rvc_inf):
704 case CLASS2 (rvc_normal, rvc_inf):
705 case CLASS2 (rvc_inf, rvc_normal):
706 /* Inf * Inf = Inf, R * Inf = Inf */
707 get_inf (r, sign);
708 return false;
709
710 case CLASS2 (rvc_normal, rvc_normal):
711 break;
712
713 default:
714 abort ();
715 }
716
717 if (r == a || r == b)
718 rr = &t;
719 else
720 rr = r;
721 get_zero (rr, 0);
722
723 /* Collect all the partial products. Since we don't have sure access
724 to a widening multiply, we split each long into two half-words.
725
726 Consider the long-hand form of a four half-word multiplication:
727
728 A B C D
729 * E F G H
730 --------------------
731 DE DF DG DH
732 CE CF CG CH
733 BE BF BG BH
734 AE AF AG AH
735
736 We construct partial products of the widened half-word products
737 that are known to not overlap, e.g. DF+DH. Each such partial
738 product is given its proper exponent, which allows us to sum them
739 and obtain the finished product. */
740
741 for (i = 0; i < SIGSZ * 2; ++i)
742 {
743 unsigned long ai = a->sig[i / 2];
744 if (i & 1)
745 ai >>= HOST_BITS_PER_LONG / 2;
746 else
747 ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
748
749 if (ai == 0)
750 continue;
751
752 for (j = 0; j < 2; ++j)
753 {
754 int exp = (a->exp - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
755 + (b->exp - (1-j)*(HOST_BITS_PER_LONG/2)));
756
757 if (exp > MAX_EXP)
758 {
759 get_inf (r, sign);
760 return true;
761 }
762 if (exp < -MAX_EXP)
763 {
764 /* Would underflow to zero, which we shouldn't bother adding. */
765 inexact = true;
766 continue;
767 }
768
769 memset (&u, 0, sizeof (u));
770 u.class = rvc_normal;
771 u.exp = exp;
772
773 for (k = j; k < SIGSZ * 2; k += 2)
774 {
775 unsigned long bi = b->sig[k / 2];
776 if (k & 1)
777 bi >>= HOST_BITS_PER_LONG / 2;
778 else
779 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
780
781 u.sig[k / 2] = ai * bi;
782 }
783
784 normalize (&u);
785 inexact |= do_add (rr, rr, &u, 0);
786 }
787 }
788
789 rr->sign = sign;
790 if (rr != r)
791 *r = t;
792
793 return inexact;
794 }

If it is normal floating point value, code after 717 line would be run. Note that though we are processing decimal floating point, but REAL_VALUE_TYPE uses Hex representation. Here the initial value 0 and multiple 10 are both in form of REAL_VALUE_TYPE. Plus lines 1872 and 1887 in real_from_string, real_digit also generates REAL_VALUE_TYPE for the readin digit.

2052 static const REAL_VALUE_TYPE *
2053 real_digit (int n) in real.c
2054 {
2055 static REAL_VALUE_TYPE num[10];
2056
2057 if (n < 0 || n > 9)
2058 abort ();
2059
2060 if (n > 0 && num
.class == rvc_zero)
2061 real_from_integer (&num
, VOIDmode, n, 0, 1);
2062
2063 return &num
;
2064 }

The significand part of REAL_VALUE_TYPE is array of long, and SIGSZ is the size of the array. So f lines rom 754 to 767, checks if the exponent part of partial products mentioned in comment from line 731 to 734 overflow, and normalize at line 784 checks the overflow for the whole partial product, and keeps precision as possible (it is very important).
Add and substract r = a + b or r = a – b in form of REAL_VALUE_TYPE is done by do_add.

522 static bool
523 do_add (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a, in real.c
524 const REAL_VALUE_TYPE *b, int subtract_p)
525 {
526 int dexp, sign, exp;
527 REAL_VALUE_TYPE t;
528 bool inexact = false;
529
530 /* Determine if we need to add or subtract. */
531 sign = a->sign;
532 subtract_p = (sign ^ b->sign) ^ subtract_p;
533
534 switch (CLASS2 (a->class, b->class))
535 {
536 case CLASS2 (rvc_zero, rvc_zero):
537 /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
538 get_zero (r, sign & !subtract_p);
539 return false;
540
541 case CLASS2 (rvc_zero, rvc_normal):
542 case CLASS2 (rvc_zero, rvc_inf):
543 case CLASS2 (rvc_zero, rvc_nan):
544 /* 0 + ANY = ANY. */
545 case CLASS2 (rvc_normal, rvc_nan):
546 case CLASS2 (rvc_inf, rvc_nan):
547 case CLASS2 (rvc_nan, rvc_nan):
548 /* ANY + NaN = NaN. */
549 case CLASS2 (rvc_normal, rvc_inf):
550 /* R + Inf = Inf. */
551 *r = *b;
552 r->sign = sign ^ subtract_p;
553 return false;
554
555 case CLASS2 (rvc_normal, rvc_zero):
556 case CLASS2 (rvc_inf, rvc_zero):
557 case CLASS2 (rvc_nan, rvc_zero):
558 /* ANY + 0 = ANY. */
559 case CLASS2 (rvc_nan, rvc_normal):
560 case CLASS2 (rvc_nan, rvc_inf):
561 /* NaN + ANY = NaN. */
562 case CLASS2 (rvc_inf, rvc_normal):
563 /* Inf + R = Inf. */
564 *r = *a;
565 return false;
566
567 case CLASS2 (rvc_inf, rvc_inf):
568 if (subtract_p)
569 /* Inf - Inf = NaN. */
570 get_canonical_qnan (r, 0);
571 else
572 /* Inf + Inf = Inf. */
573 *r = *a;
574 return false;
575
576 case CLASS2 (rvc_normal, rvc_normal):
577 break;
578
579 default:
580 abort ();
581 }
582
583 /* Swap the arguments such that A has the larger exponent. */
584 dexp = a->exp - b->exp;
585 if (dexp < 0)
586 {
587 const REAL_VALUE_TYPE *t;
588 t = a, a = b, b = t;
589 dexp = -dexp;
590 sign ^= subtract_p;
591 }
592 exp = a->exp;
593
594 /* If the exponents are not identical, we need to shift the
595 significand of B down. */
596 if (dexp > 0)
597 {
598 /* If the exponents are too far apart, the significands
599 do not overlap, which makes the subtraction a noop. */
600 if (dexp >= SIGNIFICAND_BITS)
601 {
602 *r = *a;
603 r->sign = sign;
604 return true;
605 }
606
607 inexact |= sticky_rshift_significand (&t, b, dexp);
608 b = &t;
609 }
610
611 if (subtract_p)
612 {
613 if (sub_significands (r, a, b, inexact))
614 {
615 /* We got a borrow out of the subtraction. That means that
616 A and B had the same exponent, and B had the larger
617 significand. We need to swap the sign and negate the
618 significand. */
619 sign ^= 1;
620 neg_significand (r, r);
621 }
622 }
623 else
624 {
625 if (add_significands (r, a, b))
626 {
627 /* We got carry out of the addition. This means we need to
628 shift the significand back down one bit and increase the
629 exponent. */
630 inexact |= sticky_rshift_significand (r, r, 1);
631 r->sig[SIGSZ-1] |= SIG_MSB;
632 if (++exp > MAX_EXP)
633 {
634 get_inf (r, sign);
635 return true;
636 }
637 }
638 }
639
640 r->class = rvc_normal;
641 r->sign = sign;
642 r->exp = exp;
643 /* Zero out the remaining fields. */
644 r->signalling = 0;
645 r->canonical = 0;
646
647 /* Re-normalize the result. */
648 normalize (r);
649
650 /* Special case: if the subtraction results in zero, the result
651 is positive. */
652 if (r->class == rvc_zero)
653 r->sign = 0;
654 else
655 r->sig[0] |= inexact;
656
657 return inexact;
658 }

This operation requires the operands have the same exponent part. Above lines from 585 to 591, we swap b as the value has bigger exponent part, and shift it to get the agreement from a, during which if lossing some bits, inexact will be 1, and at line 655 sets the lowest bit if result as 1 (like rounding to the nearest whole number).
As the decimal floating point’s exponent part is of power 10, thus at line 1926 in real_from_string, times_pten multiples this exponent with the significand.

2068 static void
2069 times_pten (REAL_VALUE_TYPE *r, int exp) in real.c
2070 {
2071 REAL_VALUE_TYPE pten, *rr;
2072 bool negative = (exp < 0);
2073 int i;
2074
2075 if (negative)
2076 {
2077 exp = -exp;
2078 pten = *real_digit (1);
2079 rr = &pten;
2080 }
2081 else
2082 rr = r;
2083
2084 for (i = 0; exp > 0; ++i, exp >>= 1)
2085 if (exp & 1)
2086 do_multiply (rr, rr, ten_to_ptwo (i));
2087
2088 if (negative)
2089 do_divide (r, r, &pten);
2090 }

When exponent part is negative, we put it into division, but here we don’t step into do_divide, it does binary division. When back from real_from_string, despite we can gives floating point constants in 2 forms, at last they come in Hex in form of REAL_VALUE_TYPE.式。Then at line 612 in interpret_float, real_convert will expand or trim the value according to its type. A perfect procedure.