1 /*
2     Copyright (C) 2021 Daniel Schultz
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 <https://www.gnu.org/licenses/>.
10 */
11 
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <gmp.h>
15 #include "flint.h"
16 #include "fmpz.h"
17 #include "fmpz_vec.h"
18 #include "fmpz_mat.h"
19 #include "ulong_extras.h"
20 
21 
_fmpz_mat_full(fmpz_mat_t A,flint_bitcnt_t bits)22 void _fmpz_mat_full(fmpz_mat_t A, flint_bitcnt_t bits)
23 {
24     slong i, j;
25 
26     for (i = 0; i < A->r; i++ )
27     {
28         for (j = 0; j < A->c; j++)
29         {
30             fmpz_one(fmpz_mat_entry(A, i, j));
31             fmpz_mul_2exp(fmpz_mat_entry(A, i, j), fmpz_mat_entry(A, i, j), bits);
32             fmpz_sub_ui(fmpz_mat_entry(A, i, j), fmpz_mat_entry(A, i, j), 1);
33         }
34     }
35 }
36 
37 
main(void)38 int main(void)
39 {
40     fmpz_mat_t A, B, C, D;
41     slong i;
42     FLINT_TEST_INIT(state);
43 
44     flint_printf("mul_double_word....");
45     fflush(stdout);
46 
47     for (i = 0; i < 1000 * flint_test_multiplier(); i++)
48     {
49         int sign;
50         slong m, n, k;
51         flint_bitcnt_t abits, bbits;
52 
53         sign = n_randint(state, 2);
54         m = n_randint(state, 50);
55         k = n_randint(state, 50);
56         n = n_randint(state, 50);
57         abits = n_randint(state, 2*FLINT_BITS - sign) + 1;
58         bbits = n_randint(state, 2*FLINT_BITS - sign) + 1;
59 
60         fmpz_mat_init(A, m, k);
61         fmpz_mat_init(B, k, n);
62         fmpz_mat_init(C, m, n);
63         fmpz_mat_init(D, m, n);
64 
65         if (sign)
66         {
67             if (n_randint(state, 2))
68             {
69                 _fmpz_mat_full(A, abits);
70                 _fmpz_mat_full(B, bbits);
71 
72                 if (n_randint(state, 2))
73                     fmpz_mat_neg(A, A);
74 
75                 if (n_randint(state, 2))
76                     fmpz_mat_neg(B, B);
77             }
78             else
79             {
80                 fmpz_mat_randtest(A, state, abits);
81                 fmpz_mat_randtest(B, state, bbits);
82             }
83         }
84         else
85         {
86             if (n_randint(state, 2))
87             {
88                 _fmpz_mat_full(A, abits);
89                 _fmpz_mat_full(B, bbits);
90             }
91             else
92             {
93                 fmpz_mat_randtest_unsigned(A, state, abits);
94                 fmpz_mat_randtest_unsigned(B, state, bbits);
95             }
96         }
97 
98         fmpz_mat_randtest(C, state, n_randint(state, 200) + 1);
99 
100         _fmpz_mat_mul_double_word(C, A, B);
101         fmpz_mat_mul_classical_inline(D, A, B);
102 
103         if (!fmpz_mat_equal(C, D))
104         {
105             flint_printf("FAIL: results not equal\n\n");
106             fmpz_mat_print(A); flint_printf("\n\n");
107             fmpz_mat_print(B); flint_printf("\n\n");
108             fmpz_mat_print(C); flint_printf("\n\n");
109             fmpz_mat_print(D); flint_printf("\n\n");
110             flint_abort();
111         }
112 
113         fmpz_mat_clear(A);
114         fmpz_mat_clear(B);
115         fmpz_mat_clear(C);
116         fmpz_mat_clear(D);
117     }
118 
119     FLINT_TEST_CLEANUP(state);
120 
121     flint_printf("PASS\n");
122     return 0;
123 }
124