LCOV - code coverage report
Current view: top level - tests/sieve - test_trialdiv.c (source / functions) Hit Total Coverage
Test: unnamed Lines: 0 158 0.0 %
Date: Wed Dec 13 07:00:11 2017 Functions: 0 3 0.0 %
Commit: 0a56d02

          Line data    Source code
       1             : #include "cado.h"
       2             : #include <stdio.h>
       3             : #include <stdlib.h>
       4             : #include <limits.h>
       5             : #include <sys/time.h>
       6             : #include "sieve/trialdiv.h"
       7             : #include "utils.h"
       8             : #include "test_iter.h"
       9             : #include "tests_common.h"
      10             : 
      11             : void
      12           0 : trialdiv_stdinput(const unsigned long pmax, const int verbose)
      13             : {
      14             :   unsigned long *primes;
      15             :   trialdiv_divisor_t *d;
      16             :   unsigned long *factors;
      17           0 :   int nr_primes = 0;
      18             :   unsigned long p;
      19           0 :   const size_t max_div = 32;
      20             :   mpz_t N;
      21             :   prime_info pi;
      22             : 
      23           0 :   prime_info_init (pi);
      24           0 :   for (p = getprime_mt (pi); p <= pmax; p = getprime_mt (pi))
      25           0 :     nr_primes++;
      26           0 :   primes = malloc (nr_primes * sizeof (unsigned long));
      27           0 :   prime_info_clear (pi);
      28             : 
      29           0 :   prime_info_init (pi);
      30           0 :   nr_primes = 0;
      31           0 :   for (p = getprime_mt (pi); p <= pmax; p = getprime_mt (pi))
      32           0 :     primes[nr_primes++] = p;
      33           0 :   prime_info_clear (pi);
      34             : 
      35           0 :   d = trialdiv_init (primes, nr_primes);
      36           0 :   free (primes);
      37           0 :   factors = (unsigned long *) malloc (max_div * sizeof (unsigned long));
      38             :   
      39           0 :   mpz_init(N);
      40           0 :   while (!feof(stdin)) {
      41             :     size_t t, i;
      42             :     unsigned long bit; /* old versions of GMP do not have mp_bitcnt_t */
      43             :     
      44           0 :     if (mpz_inp_str (N, stdin, 10) == 0)
      45           0 :       break;
      46             :     
      47           0 :     if (mpz_sgn(N) == 0)
      48           0 :       break;
      49             :     
      50           0 :     if (mpz_sgn(N) < 0)
      51           0 :       mpz_neg (N, N);
      52             : 
      53             :     /* Now N is positive */
      54           0 :     bit = mpz_scan1(N, 0);
      55           0 :     if (bit > 0)
      56           0 :       mpz_tdiv_q_2exp (N, N, bit);
      57             : 
      58             :     /* Now N is positive and odd */
      59           0 :     if (mpz_cmp_ui(N, 1UL) == 0) {
      60           0 :       printf ("1\n");
      61           0 :       continue;
      62             :     }
      63           0 :     t = trialdiv (factors, N, d, max_div);
      64           0 :     ASSERT_ALWAYS (t <= max_div);
      65             : 
      66           0 :     if (verbose) {
      67           0 :       for (i = 0; i < bit; i++) {
      68           0 :         printf ("2 ");
      69             :       }
      70           0 :       for (i = 0; i < t; i++)
      71           0 :         printf ("%lu ", (unsigned long) factors[i]);
      72             :     }
      73           0 :     gmp_printf ("%Zd\n", N);
      74             :   }
      75             :   
      76           0 :   mpz_clear (N);
      77           0 :   trialdiv_clear (d);
      78           0 : }
      79             : 
      80             : /* performs iter random tests with a cofactor of n limbs */
      81             : void
      82           0 : test_trialdiv (int n, unsigned long iter)
      83             : {
      84             :   mpz_t N;
      85             :   trialdiv_divisor_t *d;
      86             :   unsigned long f[1];
      87             :   unsigned long p, g[1], i, pmax;
      88             :   size_t s;
      89             :   int ret;
      90             : 
      91           0 :   mpz_init (N);
      92           0 :   pmax = trialdiv_get_max_p();
      93           0 :   for (i = 0; i < iter; i++)
      94             :     {
      95           0 :       if (i == 0)
      96           0 :         p = 3;
      97           0 :       else if (i == 1) {
      98             :         /* Find largest prime <= pmax */
      99             :         unsigned long r;
     100           0 :         for (r = pmax, p = 0; p == 0 || p > pmax; r -= 2)
     101           0 :           p = ulong_nextprime (r);
     102             :       } else {
     103           0 :           do p = ulong_nextprime (lrand48 () % pmax); while (p > pmax || p < 3);
     104             :       }
     105           0 :       f[0] = p;
     106           0 :       d = trialdiv_init (f, 1);
     107             : 
     108           0 :       mpz_urandomb (N, state, n * mp_bits_per_limb);
     109           0 :       ret = mpz_divisible_ui_p (N, p);
     110           0 :       s = trialdiv (g, N, d, 1);
     111           0 :       ASSERT_ALWAYS (s <= 2); /* s can be max_div+1, i.e., 2 */
     112           0 :       if (ret) {
     113           0 :         ASSERT_ALWAYS (s >= 1);
     114           0 :         ASSERT_ALWAYS (g[0] == p);
     115             :       } else
     116           0 :         ASSERT_ALWAYS (s == 0);
     117             : 
     118             :       /* now test a case where it should divide */
     119           0 :       mpz_sub_ui (N, N, mpz_fdiv_ui (N, p));
     120           0 :       if (mpz_sgn(N) == 0)
     121           0 :         mpz_set_ui(N, p);
     122           0 :       s = trialdiv (g, N, d, 1);
     123           0 :       ASSERT_ALWAYS (1 <= s && s <= 2); /* s can be max_div+1, i.e., 2 */
     124           0 :       ASSERT_ALWAYS (g[0] == p);
     125           0 :       trialdiv_clear(d);
     126             :     }
     127           0 :   mpz_clear (N);
     128           0 : }
     129             : 
     130           0 : int main (int argc, const char **argv)
     131             : {
     132             :   unsigned long *primes;
     133             :   trialdiv_divisor_t *d;
     134           0 :   const size_t max_div = 16;
     135             :   unsigned long *factors;
     136           0 :   int i, len = 1, nr_primes = 1000, nr_N = 100000;
     137           0 :   unsigned long expect = 0, nr_div = 0;
     138             :   mpz_t M, N, pk, t1, t2;
     139             :   double usrtime;
     140           0 :   int verbose = 0, input = 0;
     141           0 :   unsigned long iter = 10000;
     142             : 
     143           0 :   tests_common_cmdline(&argc, &argv, PARSE_SEED | PARSE_ITER);
     144           0 :   tests_common_get_iter(&iter);
     145             : 
     146             :   while (1) {
     147           0 :     if (argc > 1 && argv[1][0] == '-' && argv[1][1] == 'v')
     148             :       {
     149           0 :         verbose = 1;
     150           0 :         argc--;
     151           0 :         argv++;
     152             :       }
     153           0 :     else if (argc > 1 && argv[1][0] == '-' && argv[1][1] == 'i')
     154             :       {
     155           0 :         input = 1;
     156           0 :         argc--;
     157           0 :         argv++;
     158             :       }
     159             :     else
     160             :       break;
     161           0 :   }
     162             : 
     163           0 :   if (argc > 1)
     164           0 :     len = atoi (argv[1]);
     165             :   
     166           0 :   if (argc > 2)
     167           0 :     nr_N = atoi (argv[2]);
     168             : 
     169           0 :   if (argc > 3)
     170           0 :     nr_primes = atoi (argv[3]);
     171             :   
     172           0 :   if (input) {
     173             :     /* First parameter is pmax, and is stored in len */
     174           0 :     trialdiv_stdinput (len, verbose);
     175           0 :     exit (EXIT_SUCCESS);
     176             :   }
     177             : 
     178           0 :   if (len > TRIALDIV_MAXLEN)
     179             :     {
     180           0 :       printf ("Error, trialdiv not implemented for input sizes greater than "
     181             :               "%d words\n", TRIALDIV_MAXLEN);
     182           0 :       exit (EXIT_FAILURE);
     183             :     }
     184             : 
     185           0 :   mpz_init (M);
     186           0 :   mpz_init (N);
     187           0 :   mpz_init (pk);
     188           0 :   mpz_init (t1);
     189           0 :   mpz_init (t2);
     190           0 :   mpz_set_ui (N, 1UL);
     191           0 :   mpz_mul_2exp (N, N, 8 * sizeof(unsigned long) * len);
     192           0 :   mpz_tdiv_q_ui (N, N, 3UL);
     193           0 :   mpz_nextprime (N, N);
     194             : 
     195             :   prime_info pi;
     196           0 :   prime_info_init (pi);
     197           0 :   for (i = 0; i < 0; i++) /* To skip some primes */
     198           0 :     getprime_mt (pi);
     199           0 :   primes = malloc (nr_primes * sizeof (unsigned long));
     200             : 
     201           0 :   mpz_add_ui (M, N, nr_N - 1); /* We'll divide integers in [N, M] */
     202           0 :   mpz_sub_ui (N, N, 1UL);
     203           0 :   for (i = 0; i < nr_primes; i++)
     204             :     {
     205           0 :       primes[i] = getprime_mt (pi);
     206             :       /* Estimate the number of divisors we'll find when trial dividing by
     207             :          this p and its powers. 
     208             :          The number of times d divides in [0, n] is floor(n/d)+1,
     209             :          so the number of times it divides in [N, M] is
     210             :          floor(M/d) - floor((N-1)/d).
     211             :       */
     212           0 :       mpz_set_ui (pk, 1UL);
     213             :       do
     214             :         {
     215           0 :           mpz_mul_ui (pk, pk, primes[i]);
     216           0 :           mpz_tdiv_q (t1, N, pk);
     217           0 :           mpz_tdiv_q (t2, M, pk);
     218           0 :           mpz_sub (t2, t2, t1);
     219           0 :           ASSERT_ALWAYS (mpz_fits_ulong_p (t2));
     220           0 :           expect += mpz_get_ui (t2);
     221           0 :         } while (mpz_sgn (t2) != 0);
     222             :     }
     223           0 :   mpz_add_ui (N, N, 1UL);
     224           0 :   if (verbose)
     225           0 :     gmp_printf ("Trial dividing integers in [%Zd, %Zd] "
     226             :                 "by the %d primes in [%lu, %lu]\n", 
     227             :                 N, M, nr_primes, (unsigned long) primes[0], 
     228           0 :                 (unsigned long) primes[nr_primes - 1]);
     229           0 :   prime_info_clear (pi);
     230             : 
     231           0 :   d = trialdiv_init (primes, nr_primes);
     232           0 :   free (primes);
     233           0 :   factors = (unsigned long *) malloc (max_div * sizeof (unsigned long));
     234             :   
     235           0 :   for (i = 0; i < nr_N; i++)
     236             :     {
     237             :       size_t t;
     238           0 :       mpz_set (M, N);
     239             :       do {
     240           0 :         t = trialdiv (factors, M, d, max_div);
     241           0 :         nr_div += (t > max_div) ? max_div : t;
     242           0 :       } while (t > max_div); /* If more factors were found than fit in factor 
     243             :                                 array, continue factoring on cofactor */
     244           0 :       mpz_add_ui (N, N, 1UL);
     245             :     }
     246             :   
     247           0 :   usrtime = microseconds();
     248             : 
     249           0 :   mpz_clear (t1);
     250           0 :   mpz_clear (t2);
     251           0 :   mpz_clear (pk);
     252           0 :   mpz_clear (M);
     253           0 :   mpz_clear (N);
     254           0 :   trialdiv_clear (d);
     255           0 :   free (factors);
     256           0 :   if (verbose)
     257             :     {
     258           0 :       printf ("Found %lu divisors, expected %lu\n", nr_div, expect);
     259           0 :       printf ("Time: %f s, per N: %f mus, per prime: %f ns\n", 
     260           0 :               usrtime / 1e6, usrtime / nr_N, 
     261           0 :               usrtime / nr_N / nr_primes * 1e3);
     262             :     }
     263           0 :   if (nr_div != expect)
     264             :     {
     265           0 :       printf ("Error: did not find the expected number of divisors!\n");
     266           0 :       exit (EXIT_FAILURE);
     267             :     }
     268             : 
     269             :   /* random tests */
     270           0 :   for (int i = 1; i <= TRIALDIV_MAXLEN; i++)
     271           0 :     test_trialdiv (i, iter);
     272             : 
     273           0 :   tests_common_clear();
     274           0 :   exit (EXIT_SUCCESS);
     275             : }

Generated by: LCOV version 1.11