1 /* Tests matrix22_mul.
2
3 Copyright 2008 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library test suite.
6
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 Public License for more details.
16
17 You should have received a copy of the GNU General Public License along with
18 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
19
20 #include <stdio.h>
21 #include <stdlib.h>
22
23 #include "gmp-impl.h"
24 #include "tests.h"
25
26 struct matrix {
27 mp_size_t alloc;
28 mp_size_t n;
29 mp_ptr e00, e01, e10, e11;
30 };
31
32 static void
33 matrix_init (struct matrix *M, mp_size_t n)
34 {
35 mp_ptr p = refmpn_malloc_limbs (4*(n+1));
36 M->e00 = p; p += n+1;
37 M->e01 = p; p += n+1;
38 M->e10 = p; p += n+1;
39 M->e11 = p;
40 M->alloc = n + 1;
41 M->n = 0;
42 }
43
44 static void
45 matrix_clear (struct matrix *M)
46 {
47 refmpn_free_limbs (M->e00);
48 }
49
50 static void
51 matrix_copy (struct matrix *R, const struct matrix *M)
52 {
53 R->n = M->n;
54 MPN_COPY (R->e00, M->e00, M->n);
55 MPN_COPY (R->e01, M->e01, M->n);
56 MPN_COPY (R->e10, M->e10, M->n);
57 MPN_COPY (R->e11, M->e11, M->n);
58 }
59
60 /* Used with same size, so no need for normalization. */
61 static int
62 matrix_equal_p (const struct matrix *A, const struct matrix *B)
63 {
64 return (A->n == B->n
65 && mpn_cmp (A->e00, B->e00, A->n) == 0
66 && mpn_cmp (A->e01, B->e01, A->n) == 0
67 && mpn_cmp (A->e10, B->e10, A->n) == 0
68 && mpn_cmp (A->e11, B->e11, A->n) == 0);
69 }
70
71 static void
72 matrix_random(struct matrix *M, mp_size_t n, gmp_randstate_ptr rands)
73 {
74 M->n = n;
75 mpn_random (M->e00, n);
76 mpn_random (M->e01, n);
77 mpn_random (M->e10, n);
78 mpn_random (M->e11, n);
79 }
80
81 #define MUL(rp, ap, an, bp, bn) do { \
82 if (an > bn) \
83 mpn_mul (rp, ap, an, bp, bn); \
84 else \
85 mpn_mul (rp, bp, bn, ap, an); \
86 } while(0)
87
88 static void
89 ref_matrix22_mul (struct matrix *R,
90 const struct matrix *A,
91 const struct matrix *B, mp_ptr tp)
92 {
93 mp_size_t an, bn, n;
94 mp_ptr r00, r01, r10, r11, a00, a01, a10, a11, b00, b01, b10, b11;
95
96 if (A->n >= B->n)
97 {
98 r00 = R->e00; a00 = A->e00; b00 = B->e00;
99 r01 = R->e01; a01 = A->e01; b01 = B->e01;
100 r10 = R->e10; a10 = A->e10; b10 = B->e10;
101 r11 = R->e11; a11 = A->e11; b11 = B->e11;
102 an = A->n, bn = B->n;
103 }
104 else
105 {
106 /* Transpose */
107 r00 = R->e00; a00 = B->e00; b00 = A->e00;
108 r01 = R->e10; a01 = B->e10; b01 = A->e10;
109 r10 = R->e01; a10 = B->e01; b10 = A->e01;
110 r11 = R->e11; a11 = B->e11; b11 = A->e11;
111 an = B->n, bn = A->n;
112 }
113 n = an + bn;
114 R->n = n + 1;
115
116 mpn_mul (r00, a00, an, b00, bn);
117 mpn_mul (tp, a01, an, b10, bn);
118 r00[n] = mpn_add_n (r00, r00, tp, n);
119
120 mpn_mul (r01, a00, an, b01, bn);
121 mpn_mul (tp, a01, an, b11, bn);
122 r01[n] = mpn_add_n (r01, r01, tp, n);
123
124 mpn_mul (r10, a10, an, b00, bn);
125 mpn_mul (tp, a11, an, b10, bn);
126 r10[n] = mpn_add_n (r10, r10, tp, n);
127
128 mpn_mul (r11, a10, an, b01, bn);
129 mpn_mul (tp, a11, an, b11, bn);
130 r11[n] = mpn_add_n (r11, r11, tp, n);
131 }
132
133 static void
134 one_test (const struct matrix *A, const struct matrix *B, int i)
135 {
136 struct matrix R;
137 struct matrix P;
138 mp_ptr tp;
139
140 matrix_init (&R, A->n + B->n + 1);
141 matrix_init (&P, A->n + B->n + 1);
142
143 tp = refmpn_malloc_limbs (mpn_matrix22_mul_itch (A->n, B->n));
144
145 ref_matrix22_mul (&R, A, B, tp);
146 matrix_copy (&P, A);
147 mpn_matrix22_mul (P.e00, P.e01, P.e10, P.e11, A->n,
148 B->e00, B->e01, B->e10, B->e11, B->n, tp);
149 P.n = A->n + B->n + 1;
150 if (!matrix_equal_p (&R, &P))
151 {
152 fprintf (stderr, "ERROR in test %d\n", i);
153 gmp_fprintf (stderr, "A = (%Nx, %Nx\n %Nx, %Nx)\n"
154 "B = (%Nx, %Nx\n %Nx, %Nx)\n"
155 "R = (%Nx, %Nx (expected)\n %Nx, %Nx)\n"
156 "P = (%Nx, %Nx (incorrect)\n %Nx, %Nx)\n",
157 A->e00, A->n, A->e01, A->n, A->e10, A->n, A->e11, A->n,
158 B->e00, B->n, B->e01, B->n, B->e10, B->n, B->e11, B->n,
159 R.e00, R.n, R.e01, R.n, R.e10, R.n, R.e11, R.n,
160 P.e00, P.n, P.e01, P.n, P.e10, P.n, P.e11, P.n);
161 abort();
162 }
163 refmpn_free_limbs (tp);
164 matrix_clear (&R);
165 matrix_clear (&P);
166 }
167
168 #define MAX_SIZE (2+2*MATRIX22_STRASSEN_THRESHOLD)
169
170 int
171 main (int argc, char **argv)
172 {
173 struct matrix A;
174 struct matrix B;
175
176 gmp_randstate_ptr rands;
177 mpz_t bs;
178 int i;
179
180 tests_start ();
181 rands = RANDS;
182
183 matrix_init (&A, MAX_SIZE);
184 matrix_init (&B, MAX_SIZE);
185 mpz_init (bs);
186
187 for (i = 0; i < 1000; i++)
188 {
189 mp_size_t an, bn;
190 mpz_urandomb (bs, rands, 32);
191 an = 1 + mpz_get_ui (bs) % MAX_SIZE;
192 mpz_urandomb (bs, rands, 32);
193 bn = 1 + mpz_get_ui (bs) % MAX_SIZE;
194
195 matrix_random (&A, an, rands);
196 matrix_random (&B, bn, rands);
197
198 one_test (&A, &B, i);
199 }
200 mpz_clear (bs);
201 matrix_clear (&A);
202 matrix_clear (&B);
203
204 tests_end ();
205 return 0;
206 }