tests/test-fmod.h

   1 /* Test of fmod*() function family.
   2    Copyright (C) 2012 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 <http://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   /* Randomized tests.  */
  33   for (i = 0; i < SIZEOF (RANDOM); i++)
  34     for (j = 0; j < SIZEOF (RANDOM); j++)
  35       {
  36         DOUBLE x = L_(16.0) * RANDOM[i]; /* 0.0 <= x <= 16.0 */
  37         DOUBLE y = RANDOM[j]; /* 0.0 <= y < 1.0 */
  38         if (y > L_(0.0))
  39           {
  40             DOUBLE z = FMOD (x, y);
  41             ASSERT (z >= L_(0.0));
  42             ASSERT (z < y);
  43             z -= x - (int) (x / y) * y;
  44             ASSERT (/* The common case.  */
  45                     (z > - L_(16.0) / TWO_MANT_DIG
  46                      && z < L_(16.0) / TWO_MANT_DIG)
  47                     || /* rounding error: x / y computed too large */
  48                        (z > y - L_(16.0) / TWO_MANT_DIG
  49                         && z < y + L_(16.0) / TWO_MANT_DIG)
  50                     || /* rounding error: x / y computed too small */
  51                        (z > - y - L_(16.0) / TWO_MANT_DIG
  52                         && z < - y + L_(16.0) / TWO_MANT_DIG));
  53           }
  54       }
  55
  56   for (i = 0; i < SIZEOF (RANDOM); i++)
  57     for (j = 0; j < SIZEOF (RANDOM); j++)
  58       {
  59         DOUBLE x = L_(1.0e9) * RANDOM[i]; /* 0.0 <= x <= 10^9 */
  60         DOUBLE y = RANDOM[j]; /* 0.0 <= y < 1.0 */
  61         if (y > L_(0.0))
  62           {
  63             DOUBLE z = FMOD (x, y);
  64             DOUBLE r;
  65             ASSERT (z >= L_(0.0));
  66             ASSERT (z < y);
  67             {
  68               /* Determine the quotient x / y in two steps, because it
  69                  may be > 2^31.  */
  70               int q1 = (int) (x / y / L_(65536.0));
  71               int q2 = (int) ((x - q1 * L_(65536.0) * y) / y);
  72               DOUBLE q = (DOUBLE) q1 * L_(65536.0) + (DOUBLE) q2;
  73               r = x - q * y;
  74             }
  75             /* The absolute error of z can be up to 1e9/2^MANT_DIG.
  76                The absolute error of r can also be up to 1e9/2^MANT_DIG.
  77                Therefore the error of z - r can be twice as large.  */
  78             z -= r;
  79             ASSERT (/* The common case.  */
  80                     (z > - L_(2.0) * L_(1.0e9) / TWO_MANT_DIG
  81                      && z < L_(2.0) * L_(1.0e9) / TWO_MANT_DIG)
  82                     || /* rounding error: x / y computed too large */
  83                        (z > y - L_(2.0) * L_(1.0e9) / TWO_MANT_DIG
  84                         && z < y + L_(2.0) * L_(1.0e9) / TWO_MANT_DIG)
  85                     || /* rounding error: x / y computed too small */
  86                        (z > - y - L_(2.0) * L_(1.0e9) / TWO_MANT_DIG
  87                         && z < - y + L_(2.0) * L_(1.0e9) / TWO_MANT_DIG));
  88           }
  89       }
  90 }
  91
  92 volatile DOUBLE x;
  93 volatile DOUBLE y;
  94 DOUBLE z;