1  /* Test of log10*() function family.
       2     Copyright (C) 2012-2023 Free Software Foundation, Inc.
       3  
       4     This program is free software: you can redistribute it and/or modify
       5     it under the terms of the GNU General Public License as published by
       6     the Free Software Foundation, either version 3 of the License, or
       7     (at your option) any later version.
       8  
       9     This program is distributed in the hope that it will be useful,
      10     but WITHOUT ANY WARRANTY; without even the implied warranty of
      11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      12     GNU General Public License for more details.
      13  
      14     You should have received a copy of the GNU General Public License
      15     along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
      16  
      17  static void
      18  test_function (void)
      19  {
      20    int i;
      21    int j;
      22    const DOUBLE TWO_MANT_DIG =
      23      /* Assume MANT_DIG <= 5 * 31.
      24         Use the identity
      25           n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5).  */
      26      (DOUBLE) (1U << ((MANT_DIG - 1) / 5))
      27      * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5))
      28      * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5))
      29      * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5))
      30      * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5));
      31  
      32    /* Pole.  */
      33    {
      34      DOUBLE z = LOG10 (L_(0.0));
      35      ASSERT (z == - HUGEVAL);
      36    }
      37    {
      38      DOUBLE z = LOG10 (MINUS_ZERO);
      39      ASSERT (z == - HUGEVAL);
      40    }
      41  
      42    /* Randomized tests.  */
      43    {
      44      /* Error bound, in ulps.  */
      45      const DOUBLE err_bound =
      46        (sizeof (DOUBLE) > sizeof (double) ?
      47  #if defined __i386__ && defined __FreeBSD__
      48         /* On FreeBSD/x86 6.4, the 'long double' type really has only 53 bits of
      49            precision in the compiler but 64 bits of precision at runtime.  See
      50            <https://lists.gnu.org/r/bug-gnulib/2008-07/msg00063.html>.
      51            The compiler has truncated all 'long double' literals in log10l.c to
      52            53 bits of precision.  */
      53         L_(18.0)
      54  #else
      55         L_(3.0)
      56  #endif
      57         : L_(3.0));
      58  
      59      for (i = 0; i < SIZEOF (RANDOM); i++)
      60        {
      61          DOUBLE x = L_(16.0) * RANDOM[i] + L_(1.0); /* 1.0 <= x <= 17.0 */
      62          DOUBLE y = LOG10 (x);
      63          DOUBLE z = LOG10 (L_(1.0) / x);
      64          DOUBLE err = y + z;
      65          ASSERT (y >= L_(0.0));
      66          ASSERT (z <= L_(0.0));
      67          ASSERT (err > - err_bound / TWO_MANT_DIG
      68                  && err < err_bound / TWO_MANT_DIG);
      69        }
      70    }
      71  
      72    {
      73      /* Error bound, in ulps.  */
      74      const DOUBLE err_bound =
      75        (sizeof (DOUBLE) > sizeof (double) ?
      76  #if defined __i386__ && defined __FreeBSD__
      77         /* On FreeBSD/x86 6.4, the 'long double' type really has only 53 bits of
      78            precision in the compiler but 64 bits of precision at runtime.  See
      79            <https://lists.gnu.org/r/bug-gnulib/2008-07/msg00063.html>.
      80            The compiler has truncated all 'long double' literals in log10l.c to
      81            53 bits of precision.  */
      82         L_(38.0)
      83  #else
      84         L_(5.0)
      85  #endif
      86         : L_(5.0));
      87  
      88      for (i = 0; i < SIZEOF (RANDOM) / 5; i++)
      89        for (j = 0; j < SIZEOF (RANDOM) / 5; j++)
      90          {
      91            DOUBLE x = L_(17.0) / (L_(16.0) - L_(15.0) * RANDOM[i]) - L_(1.0);
      92            DOUBLE y = L_(17.0) / (L_(16.0) - L_(15.0) * RANDOM[j]) - L_(1.0);
      93            /* 1/16 <= x,y <= 16 */
      94            DOUBLE z = L_(1.0) / (x * y);
      95            /* Approximately  x * y * z = 1.  */
      96            DOUBLE err = LOG10 (x) + LOG10 (y) + LOG10 (z);
      97            ASSERT (err > - err_bound / TWO_MANT_DIG
      98                    && err < err_bound / TWO_MANT_DIG);
      99          }
     100    }
     101  }
     102  
     103  volatile DOUBLE x;
     104  DOUBLE y;