tests/test-log1p.h

   1 /* Test of log1p*() function family.
   2    Copyright (C) 2012-2024 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 = LOG1P (L_(-1.0));
  35     ASSERT (z == - HUGEVAL);
  36   }
  37
  38   /* Randomized tests.  */
  39   {
  40     /* Error bound, in ulps.  */
  41     const DOUBLE err_bound =
  42       (sizeof (DOUBLE) > sizeof (double) ?
  43 #if defined __i386__ && defined __FreeBSD__
  44        /* On FreeBSD/x86 6.4, the 'long double' type really has only 53 bits of
  45           precision in the compiler but 64 bits of precision at runtime.  See
  46           <https://lists.gnu.org/r/bug-gnulib/2008-07/msg00063.html>.
  47           The compiler has truncated all 'long double' literals in log1pl.c to
  48           53 bits of precision.  */
  49        L_(900.0)
  50 #else
  51        L_(26.0)
  52 #endif
  53        : L_(13.0));
  54
  55     for (i = 0; i < SIZEOF (RANDOM); i++)
  56       {
  57         DOUBLE x = L_(16.0) * RANDOM[i]; /* 0.0 <= x <= 16.0 */
  58         DOUBLE y = LOG1P (x);
  59         DOUBLE z = LOG1P (- x / (L_(1.0) + x));
  60         DOUBLE err = y + z;
  61         ASSERT (y >= L_(0.0));
  62         ASSERT (z <= L_(0.0));
  63         ASSERT (err > - err_bound / TWO_MANT_DIG
  64                 && err < err_bound / TWO_MANT_DIG);
  65       }
  66   }
  67
  68   {
  69     /* Error bound, in ulps.  */
  70     const DOUBLE err_bound =
  71       (sizeof (DOUBLE) > sizeof (double) ?
  72 #if defined __i386__ && defined __FreeBSD__
  73        /* On FreeBSD/x86 6.4, the 'long double' type really has only 53 bits of
  74           precision in the compiler but 64 bits of precision at runtime.  See
  75           <https://lists.gnu.org/r/bug-gnulib/2008-07/msg00063.html>.
  76           The compiler has truncated all 'long double' literals in log1pl.c to
  77           53 bits of precision.  */
  78        L_(1020.0)
  79 #else
  80        L_(65.0)
  81 #endif
  82        : L_(61.0));
  83
  84     for (i = 0; i < SIZEOF (RANDOM) / 5; i++)
  85       for (j = 0; j < SIZEOF (RANDOM) / 5; j++)
  86         {
  87           DOUBLE x = L_(17.0) / (L_(16.0) - L_(15.0) * RANDOM[i]) - L_(2.0);
  88           DOUBLE y = L_(17.0) / (L_(16.0) - L_(15.0) * RANDOM[j]) - L_(2.0);
  89           /* -15/16 <= x,y <= 15 */
  90           DOUBLE z = L_(1.0) / ((L_(1.0) + x) * (L_(1.0) + y)) - L_(1.0);
  91           /* Approximately  (1+x) * (1+y) * (1+z) = 1.  */
  92           DOUBLE err = LOG1P (x) + LOG1P (y) + LOG1P (z);
  93           ASSERT (err > - err_bound / TWO_MANT_DIG
  94                   && err < err_bound / TWO_MANT_DIG);
  95         }
  96   }
  97 }
  98
  99 volatile DOUBLE x;
 100 DOUBLE y;