tests/mpz/t-root.c

   1 /* Test mpz_root, mpz_rootrem, and mpz_perfect_power_p.
   2
   3 Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2009 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22
  23 #include "gmp.h"
  24 #include "gmp-impl.h"
  25 #include "tests.h"
  26
  27 void debug_mp __GMP_PROTO ((mpz_t, int));
  28
  29 void
  30 check_one (mpz_t root1, mpz_t x2, unsigned long nth, int i)
  31 {
  32   mpz_t temp, temp2;
  33   mpz_t root2, rem2;
  34
  35   mpz_init (root2);
  36   mpz_init (rem2);
  37   mpz_init (temp);
  38   mpz_init (temp2);
  39
  40   MPZ_CHECK_FORMAT (root1);
  41
  42   mpz_rootrem (root2, rem2, x2, nth);
  43   MPZ_CHECK_FORMAT (root2);
  44   MPZ_CHECK_FORMAT (rem2);
  45
  46   mpz_pow_ui (temp, root1, nth);
  47   MPZ_CHECK_FORMAT (temp);
  48
  49   mpz_add (temp2, temp, rem2);
  50
  51   /* Is power of result > argument?  */
  52   if (mpz_cmp (root1, root2) != 0 || mpz_cmp (x2, temp2) != 0 || mpz_cmp (temp, x2) > 0)
  53     {
  54       fprintf (stderr, "ERROR after test %d\n", i);
  55       debug_mp (x2, 10);
  56       debug_mp (root1, 10);
  57       debug_mp (root2, 10);
  58       fprintf (stderr, "nth: %lu\n", nth);
  59       abort ();
  60     }
  61
  62   if (nth > 1 && mpz_cmp_ui (temp, 1L) > 0 && ! mpz_perfect_power_p (temp))
  63     {
  64       fprintf (stderr, "ERROR in mpz_perfect_power_p after test %d\n", i);
  65       debug_mp (temp, 10);
  66       debug_mp (root1, 10);
  67       fprintf (stderr, "nth: %lu\n", nth);
  68       abort ();
  69     }
  70
  71   if (nth <= 10000)             /* skip too expensive test */
  72     {
  73       mpz_add_ui (temp2, root1, 1L);
  74       mpz_pow_ui (temp2, temp2, nth);
  75       MPZ_CHECK_FORMAT (temp2);
  76
  77       /* Is square of (result + 1) <= argument?  */
  78       if (mpz_cmp (temp2, x2) <= 0)
  79         {
  80           fprintf (stderr, "ERROR after test %d\n", i);
  81           debug_mp (x2, 10);
  82           debug_mp (root1, 10);
  83           fprintf (stderr, "nth: %lu\n", nth);
  84           abort ();
  85         }
  86     }
  87
  88   mpz_clear (root2);
  89   mpz_clear (rem2);
  90   mpz_clear (temp);
  91   mpz_clear (temp2);
  92 }
  93
  94 int
  95 main (int argc, char **argv)
  96 {
  97   mpz_t x2;
  98   mpz_t root1;
  99   mp_size_t x2_size;
 100   int i;
 101   int reps = 500;
 102   unsigned long nth;
 103   gmp_randstate_ptr rands;
 104   mpz_t bs;
 105   unsigned long bsi, size_range;
 106
 107   tests_start ();
 108   TESTS_REPS (reps, argv, argc);
 109
 110   rands = RANDS;
 111
 112   mpz_init (bs);
 113
 114   mpz_init (x2);
 115   mpz_init (root1);
 116
 117   /* This triggers a gcc 4.3.2 bug */
 118   mpz_set_str (x2, "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff80000000000000000000000000000000000000000000000000000000000000002", 16);
 119   mpz_root (root1, x2, 2);
 120   check_one (root1, x2, 2, -1);
 121
 122   for (i = 0; i < reps; i++)
 123     {
 124       mpz_urandomb (bs, rands, 32);
 125       size_range = mpz_get_ui (bs) % 17 + 2;
 126
 127       mpz_urandomb (bs, rands, size_range);
 128       x2_size = mpz_get_ui (bs) + 10;
 129       mpz_rrandomb (x2, rands, x2_size);
 130
 131       mpz_urandomb (bs, rands, 15);
 132       nth = mpz_getlimbn (bs, 0) % mpz_sizeinbase (x2, 2) + 2;
 133
 134       mpz_root (root1, x2, nth);
 135
 136       mpz_urandomb (bs, rands, 4);
 137       bsi = mpz_get_ui (bs);
 138       if ((bsi & 1) != 0)
 139         {
 140           /* With 50% probability, set x2 near a perfect power.  */
 141           mpz_pow_ui (x2, root1, nth);
 142           if ((bsi & 2) != 0)
 143             {
 144               mpz_sub_ui (x2, x2, bsi >> 2);
 145               mpz_abs (x2, x2);
 146             }
 147           else
 148             mpz_add_ui (x2, x2, bsi >> 2);
 149           mpz_root (root1, x2, nth);
 150         }
 151
 152       check_one (root1, x2, nth, i);
 153     }
 154
 155   mpz_clear (bs);
 156   mpz_clear (x2);
 157   mpz_clear (root1);
 158
 159   tests_end ();
 160   exit (0);
 161 }
 162
 163 void
 164 debug_mp (mpz_t x, int base)
 165 {
 166   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
 167 }