1 /*
2     Copyright (C) 2009, 2015 William Hart
3 
4     This file is part of FLINT.
5 
6     FLINT is free software: you can redistribute it and/or modify it under
7     the terms of the GNU Lesser General Public License (LGPL) as published
8     by the Free Software Foundation; either version 2.1 of the License, or
9     (at your option) any later version.  See <http://www.gnu.org/licenses/>.
10 */
11 
12 #include <gmp.h>
13 #include "flint.h"
14 #include "ulong_extras.h"
15 #include "fmpz.h"
16 
17 int
main(void)18 main(void)
19 {
20     int i, result;
21     fmpz_t n;
22     FLINT_TEST_INIT(state);
23 
24     flint_printf("ll_mod_preinv....");
25     fflush(stdout);
26 
27     /* check (m*d + r1) % d == r1 */
28     for (i = 0; i < 100000 * flint_test_multiplier(); i++)
29     {
30         ulong d, dinv, nh, nl, r1, r2, m;
31 
32         d = n_randtest_not_zero(state);
33         m = n_randtest(state);
34         r1 = n_randtest(state) % d;
35 
36         /* <nh, nl> = m*d + r1 < (m + 1)*d <= 2^FLINT_BITS * d */
37 
38         umul_ppmm(nh, nl, m, d);
39         add_ssaaaa(nh, nl, nh, nl, UWORD(0), r1);
40 
41         dinv = n_preinvert_limb(d);
42 
43         r2 = n_ll_mod_preinv(nh, nl, d, dinv);
44 
45         result = (r1 == r2);
46         if (!result)
47         {
48             flint_printf("FAIL:\n");
49             flint_printf("<nh, nl> = (m*d + r1) but <nh, nl> % d != r1\n");
50             flint_printf("nh = %wu, nl = %wu, d = %wu, dinv = %wu\n", nh, nl,
51                          d, dinv);
52             flint_printf("r1 = %wu, r2 = %wu\n", r1, r2);
53             abort();
54         }
55     }
56 
57     fmpz_init(n);
58 
59     /* compare with fmpz_fdiv_ui */
60     for (i = 0; i < 100000 * flint_test_multiplier(); i++)
61     {
62         ulong d, dinv, nh, nl, r1, r2;
63 
64         d = n_randtest_not_zero(state);
65         nh = n_randtest(state);
66         nl = n_randtest(state);
67 
68         dinv = n_preinvert_limb(d);
69 
70         /* n = <nh, nl> */
71         fmpz_set_ui(n, nh);
72         fmpz_mul_2exp(n, n, FLINT_BITS);
73         fmpz_add_ui(n, n, nl);
74 
75         r1 = n_ll_mod_preinv(nh, nl, d, dinv);
76 
77         r2 = fmpz_fdiv_ui(n, d);
78 
79         result = (r1 == r2);
80         if (!result)
81         {
82             flint_printf("FAIL:\n");
83             flint_printf
84                 ("n = <nh, nl> but n % d does not agree with fmpz_fdiv_ui\n");
85             flint_printf("nh = %wu, nl = %wu, d = %wu, dinv = %wu\n", nh, nl,
86                          d, dinv);
87             flint_printf("r1 = %wu, r2 = %wu\n", r1, r2);
88             fmpz_clear(n);
89             abort();
90         }
91     }
92 
93     fmpz_clear(n);
94 
95     FLINT_TEST_CLEANUP(state);
96 
97     flint_printf("PASS\n");
98     return 0;
99 }
100