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[n].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[n], VOIDmode, t, 0, 1);

2023     }

2024     else

2025     {

2026        const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);

2027        do_multiply (&tens[n], t, t);

2028     }

2029   }

2030

2031   return &tens[n];

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[n].class == rvc_zero)

2061     real_from_integer (&num[n], VOIDmode, n, 0, 1);

2062

2063   return &num[n];

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.