1 /*
2 Copyright (C) 2011 Fredrik Johansson
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 "fmpz_mat.h"
13 #include "longlong.h"
14
15 void
fmpz_mat_mul_classical_inline(fmpz_mat_t C,const fmpz_mat_t A,const fmpz_mat_t B)16 fmpz_mat_mul_classical_inline(fmpz_mat_t C, const fmpz_mat_t A,
17 const fmpz_mat_t B)
18 {
19 slong ar, bc, br;
20 slong i, j, k;
21
22 fmpz a, b;
23 mpz_t t;
24
25 mp_limb_t au, bu;
26 mp_limb_t pos[3];
27 mp_limb_t neg[3];
28
29 ar = A->r;
30 br = B->r;
31 bc = B->c;
32
33 mpz_init(t);
34
35 for (i = 0; i < ar; i++)
36 {
37 for (j = 0; j < bc; j++)
38 {
39 flint_mpz_set_ui(t, UWORD(0));
40
41 pos[2] = pos[1] = pos[0] = neg[2] = neg[1] = neg[0] = UWORD(0);
42
43 for (k = 0; k < br; k++)
44 {
45 a = A->rows[i][k];
46 b = B->rows[k][j];
47
48 if (a == 0 || b == 0)
49 continue;
50
51 if (!COEFF_IS_MPZ(a)) /* a is small */
52 {
53 if (!COEFF_IS_MPZ(b)) /* both are small */
54 {
55 au = FLINT_ABS(a);
56 bu = FLINT_ABS(b);
57
58 umul_ppmm(au, bu, au, bu);
59
60 if ((a ^ b) >= WORD(0))
61 add_sssaaaaaa(pos[2], pos[1], pos[0],
62 pos[2], pos[1], pos[0], 0, au, bu);
63 else
64 add_sssaaaaaa(neg[2], neg[1], neg[0],
65 neg[2], neg[1], neg[0], 0, au, bu);
66 }
67 else
68 {
69 if (a >= 0)
70 flint_mpz_addmul_ui(t, COEFF_TO_PTR(b), a);
71 else
72 flint_mpz_submul_ui(t, COEFF_TO_PTR(b), -a);
73 }
74 }
75 else if (!COEFF_IS_MPZ(b)) /* b is small */
76 {
77 if (b >= 0)
78 flint_mpz_addmul_ui(t, COEFF_TO_PTR(a), b);
79 else
80 flint_mpz_submul_ui(t, COEFF_TO_PTR(a), -b);
81 }
82 else
83 {
84 mpz_addmul(t, COEFF_TO_PTR(a), COEFF_TO_PTR(b));
85 }
86 }
87
88 if (mpz_sgn(t) != 0 || pos[2] || neg[2] || pos[1] || neg[1])
89 {
90 __mpz_struct r;
91
92 r._mp_size = pos[2] ? 3 : (pos[1] ? 2 : pos[0] != 0);
93 r._mp_alloc = r._mp_size;
94 r._mp_d = pos;
95
96 mpz_add(t, t, &r);
97
98 r._mp_size = neg[2] ? 3 : (neg[1] ? 2 : neg[0] != 0);
99 r._mp_alloc = r._mp_size;
100 r._mp_d = neg;
101
102 mpz_sub(t, t, &r);
103
104 fmpz_set_mpz(fmpz_mat_entry(C, i, j), t);
105 }
106 else
107 {
108 if (neg[0] > pos[0])
109 fmpz_neg_ui(fmpz_mat_entry(C, i, j), neg[0] - pos[0]);
110 else
111 fmpz_set_ui(fmpz_mat_entry(C, i, j), pos[0] - neg[0]);
112 }
113 }
114 }
115
116 mpz_clear(t);
117 }
118