[f613a24] | 1 | /************************************************* |
---|
| 2 | * Author: Sebastian Jambor, 2011 * |
---|
| 3 | * GPL (e-mail from June 6, 2012, 17:00:31 MESZ) * |
---|
| 4 | * sebastian@momo.math.rwth-aachen.de * |
---|
| 5 | ************************************************/ |
---|
[2024f6a] | 6 | |
---|
| 7 | |
---|
| 8 | #include<cmath> |
---|
[056089] | 9 | #include <cstdlib> |
---|
[762407] | 10 | #include "config.h" |
---|
[f5b40a] | 11 | #include<kernel/mod2.h> |
---|
[f613a24] | 12 | |
---|
[41b5c0] | 13 | //#include<iomanip> |
---|
[2024f6a] | 14 | |
---|
| 15 | #include "minpoly.h" |
---|
| 16 | |
---|
| 17 | // TODO: avoid code copying, use subclassing instead! |
---|
| 18 | |
---|
[5b2745] | 19 | LinearDependencyMatrix::LinearDependencyMatrix (unsigned n, unsigned long p) |
---|
| 20 | { |
---|
| 21 | this->n = n; |
---|
| 22 | this->p = p; |
---|
| 23 | |
---|
| 24 | matrix = new unsigned long *[n]; |
---|
[abcc692] | 25 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 26 | { |
---|
| 27 | matrix[i] = new unsigned long[2 * n + 1]; |
---|
| 28 | } |
---|
| 29 | pivots = new unsigned[n]; |
---|
| 30 | tmprow = new unsigned long[2 * n + 1]; |
---|
| 31 | rows = 0; |
---|
[2024f6a] | 32 | } |
---|
| 33 | |
---|
[5b2745] | 34 | LinearDependencyMatrix::~LinearDependencyMatrix () |
---|
| 35 | { |
---|
| 36 | delete[]tmprow; |
---|
| 37 | delete[]pivots; |
---|
[2024f6a] | 38 | |
---|
[abcc692] | 39 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 40 | { |
---|
| 41 | delete[](matrix[i]); |
---|
| 42 | } |
---|
| 43 | delete[]matrix; |
---|
[2024f6a] | 44 | } |
---|
| 45 | |
---|
[5b2745] | 46 | void LinearDependencyMatrix::resetMatrix () |
---|
| 47 | { |
---|
| 48 | rows = 0; |
---|
[2024f6a] | 49 | } |
---|
| 50 | |
---|
[5b2745] | 51 | int LinearDependencyMatrix::firstNonzeroEntry (unsigned long *row) |
---|
| 52 | { |
---|
[abcc692] | 53 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 54 | if(row[i] != 0) |
---|
| 55 | return i; |
---|
[2024f6a] | 56 | |
---|
[5b2745] | 57 | return -1; |
---|
[2024f6a] | 58 | } |
---|
| 59 | |
---|
[5b2745] | 60 | void LinearDependencyMatrix::reduceTmpRow () |
---|
| 61 | { |
---|
[abcc692] | 62 | for(int i = 0; i < rows; i++) |
---|
[5b2745] | 63 | { |
---|
| 64 | unsigned piv = pivots[i]; |
---|
[abcc692] | 65 | unsigned x = tmprow[piv]; |
---|
[5b2745] | 66 | // if the corresponding entry in the row is zero, there is nothing to do |
---|
| 67 | if(x != 0) |
---|
| 68 | { |
---|
| 69 | // subtract tmprow[i] times the i-th row |
---|
[abcc692] | 70 | for(int j = piv; j < n + rows + 1; j++) |
---|
[5b2745] | 71 | { |
---|
[f613a24] | 72 | if (matrix[i][j] != 0) |
---|
| 73 | { |
---|
| 74 | unsigned long tmp = multMod (matrix[i][j], x, p); |
---|
| 75 | tmp = p - tmp; |
---|
| 76 | tmprow[j] += tmp; |
---|
| 77 | if (tmprow[j] >= p) |
---|
| 78 | { |
---|
| 79 | tmprow[j] -= p; |
---|
| 80 | } |
---|
| 81 | } |
---|
[5b2745] | 82 | } |
---|
| 83 | } |
---|
| 84 | } |
---|
[2024f6a] | 85 | } |
---|
| 86 | |
---|
| 87 | |
---|
[5b2745] | 88 | void LinearDependencyMatrix::normalizeTmp (unsigned i) |
---|
| 89 | { |
---|
| 90 | unsigned long inv = modularInverse (tmprow[i], p); |
---|
| 91 | tmprow[i] = 1; |
---|
[abcc692] | 92 | for(int j = i + 1; j < 2 * n + 1; j++) |
---|
| 93 | tmprow[j] = multMod (tmprow[j], inv, p); |
---|
[2024f6a] | 94 | } |
---|
[ec9db9] | 95 | |
---|
[5b2745] | 96 | bool LinearDependencyMatrix::findLinearDependency (unsigned long *newRow, |
---|
| 97 | unsigned long *dep) |
---|
| 98 | { |
---|
| 99 | // Copy newRow to tmprow (we need to add RHS) |
---|
[abcc692] | 100 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 101 | { |
---|
| 102 | tmprow[i] = newRow[i]; |
---|
| 103 | tmprow[n + i] = 0; |
---|
| 104 | } |
---|
| 105 | tmprow[2 * n] = 0; |
---|
| 106 | tmprow[n + rows] = 1; |
---|
| 107 | |
---|
| 108 | reduceTmpRow (); |
---|
| 109 | |
---|
| 110 | // Is tmprow reduced to zero? Then we have found a linear dependence. |
---|
| 111 | // Otherwise, we just add tmprow to the matrix. |
---|
[abcc692] | 112 | unsigned newpivot = firstNonzeroEntry (tmprow); |
---|
[5b2745] | 113 | if(newpivot == -1) |
---|
| 114 | { |
---|
[abcc692] | 115 | for(int i = 0; i <= n; i++) |
---|
[5b2745] | 116 | { |
---|
| 117 | dep[i] = tmprow[n + i]; |
---|
| 118 | } |
---|
| 119 | |
---|
| 120 | return true; |
---|
| 121 | } |
---|
| 122 | else |
---|
| 123 | { |
---|
| 124 | normalizeTmp (newpivot); |
---|
| 125 | |
---|
[abcc692] | 126 | for(int i = 0; i < 2 * n + 1; i++) |
---|
[5b2745] | 127 | { |
---|
| 128 | matrix[rows][i] = tmprow[i]; |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | pivots[rows] = newpivot; |
---|
| 132 | rows++; |
---|
| 133 | |
---|
| 134 | return false; |
---|
| 135 | } |
---|
[2024f6a] | 136 | } |
---|
| 137 | |
---|
[ec9db9] | 138 | #if 0 |
---|
[5b2745] | 139 | std::ostream & operator<< (std::ostream & out, |
---|
| 140 | const LinearDependencyMatrix & mat) |
---|
| 141 | { |
---|
| 142 | int width = ((int) log10 (mat.p)) + 1; |
---|
[2024f6a] | 143 | |
---|
[abcc692] | 144 | out << "Pivots: " << std::endl; |
---|
| 145 | for(int j = 0; j < mat.n; j++) |
---|
| 146 | { |
---|
| 147 | out << std::setw (width) << mat.pivots[j] << " "; |
---|
| 148 | } |
---|
| 149 | out << std::endl; |
---|
| 150 | out << "matrix:" << std::endl; |
---|
[5b2745] | 151 | for(int i = 0; i < mat.rows; i++) |
---|
| 152 | { |
---|
| 153 | for(int j = 0; j < mat.n; j++) |
---|
| 154 | { |
---|
| 155 | out << std::setw (width) << mat.matrix[i][j] << " "; |
---|
[2024f6a] | 156 | } |
---|
| 157 | out << "| "; |
---|
[5b2745] | 158 | for(int j = mat.n; j <= 2 * mat.n; j++) |
---|
| 159 | { |
---|
| 160 | out << std::setw (width) << mat.matrix[i][j] << " "; |
---|
[2024f6a] | 161 | } |
---|
| 162 | out << std::endl; |
---|
[5b2745] | 163 | } |
---|
[abcc692] | 164 | out << "tmprow: " << std::endl; |
---|
[5b2745] | 165 | for(int j = 0; j < mat.n; j++) |
---|
| 166 | { |
---|
| 167 | out << std::setw (width) << mat.tmprow[j] << " "; |
---|
| 168 | } |
---|
| 169 | out << "| "; |
---|
| 170 | for(int j = mat.n; j <= 2 * mat.n; j++) |
---|
| 171 | { |
---|
| 172 | out << std::setw (width) << mat.tmprow[j] << " "; |
---|
| 173 | } |
---|
| 174 | out << std::endl; |
---|
| 175 | |
---|
| 176 | return out; |
---|
[2024f6a] | 177 | } |
---|
[ec9db9] | 178 | #endif |
---|
[2024f6a] | 179 | |
---|
| 180 | |
---|
[5b2745] | 181 | NewVectorMatrix::NewVectorMatrix (unsigned n, unsigned long p) |
---|
| 182 | { |
---|
| 183 | this->n = n; |
---|
| 184 | this->p = p; |
---|
[2024f6a] | 185 | |
---|
[5b2745] | 186 | matrix = new unsigned long *[n]; |
---|
[abcc692] | 187 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 188 | { |
---|
| 189 | matrix[i] = new unsigned long[n]; |
---|
| 190 | } |
---|
[2024f6a] | 191 | |
---|
[5b2745] | 192 | pivots = new unsigned[n]; |
---|
[f613a24] | 193 | |
---|
| 194 | nonPivots = new unsigned[n]; |
---|
| 195 | |
---|
| 196 | for (int i = 0; i < n; i++) |
---|
| 197 | { |
---|
| 198 | nonPivots[i] = i; |
---|
| 199 | } |
---|
| 200 | |
---|
[5b2745] | 201 | rows = 0; |
---|
[2024f6a] | 202 | } |
---|
| 203 | |
---|
[5b2745] | 204 | NewVectorMatrix::~NewVectorMatrix () |
---|
| 205 | { |
---|
[f613a24] | 206 | delete nonPivots; |
---|
[5b2745] | 207 | delete pivots; |
---|
[2024f6a] | 208 | |
---|
[abcc692] | 209 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 210 | { |
---|
| 211 | delete[]matrix[i]; |
---|
| 212 | } |
---|
| 213 | delete matrix; |
---|
[2024f6a] | 214 | } |
---|
| 215 | |
---|
[5b2745] | 216 | int NewVectorMatrix::firstNonzeroEntry (unsigned long *row) |
---|
| 217 | { |
---|
[abcc692] | 218 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 219 | if(row[i] != 0) |
---|
| 220 | return i; |
---|
[2024f6a] | 221 | |
---|
[5b2745] | 222 | return -1; |
---|
[2024f6a] | 223 | } |
---|
| 224 | |
---|
[5b2745] | 225 | void NewVectorMatrix::normalizeRow (unsigned long *row, unsigned i) |
---|
| 226 | { |
---|
| 227 | unsigned long inv = modularInverse (row[i], p); |
---|
| 228 | row[i] = 1; |
---|
| 229 | |
---|
[abcc692] | 230 | for(int j = i + 1; j < n; j++) |
---|
[5b2745] | 231 | { |
---|
[abcc692] | 232 | row[j] = multMod (row[j], inv, p); |
---|
[5b2745] | 233 | } |
---|
[2024f6a] | 234 | } |
---|
| 235 | |
---|
[5b2745] | 236 | void NewVectorMatrix::insertRow (unsigned long *row) |
---|
| 237 | { |
---|
[abcc692] | 238 | for(int i = 0; i < rows; i++) |
---|
[5b2745] | 239 | { |
---|
| 240 | unsigned piv = pivots[i]; |
---|
[abcc692] | 241 | unsigned x = row[piv]; |
---|
[5b2745] | 242 | // if the corresponding entry in the row is zero, there is nothing to do |
---|
| 243 | if(x != 0) |
---|
| 244 | { |
---|
| 245 | // subtract row[i] times the i-th row |
---|
[f613a24] | 246 | // only the non-pivot entries have to be considered |
---|
| 247 | // (and also the first entry) |
---|
| 248 | row[piv] = 0; |
---|
| 249 | |
---|
| 250 | int smallestNonPivIndex = 0; |
---|
| 251 | while (nonPivots[smallestNonPivIndex] < piv) |
---|
[5b2745] | 252 | { |
---|
[f613a24] | 253 | smallestNonPivIndex++; |
---|
| 254 | } |
---|
| 255 | |
---|
| 256 | for (int j = smallestNonPivIndex; j < n-rows; j++) |
---|
| 257 | { |
---|
| 258 | unsigned ind = nonPivots[j]; |
---|
| 259 | if (matrix[i][ind] != 0) |
---|
| 260 | { |
---|
| 261 | unsigned long tmp = multMod (matrix[i][ind], x, p); |
---|
| 262 | tmp = p - tmp; |
---|
| 263 | row[ind] += tmp; |
---|
| 264 | if (row[ind] >= p) |
---|
| 265 | { |
---|
| 266 | row[ind] -= p; |
---|
| 267 | } |
---|
| 268 | } |
---|
[5b2745] | 269 | } |
---|
| 270 | } |
---|
| 271 | } |
---|
| 272 | |
---|
[abcc692] | 273 | unsigned piv = firstNonzeroEntry (row); |
---|
[5b2745] | 274 | |
---|
| 275 | if(piv != -1) |
---|
| 276 | { |
---|
[f613a24] | 277 | // Normalize and insert row into the matrix. |
---|
| 278 | // Then reduce upwards. |
---|
[5b2745] | 279 | normalizeRow (row, piv); |
---|
[abcc692] | 280 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 281 | { |
---|
| 282 | matrix[rows][i] = row[i]; |
---|
| 283 | } |
---|
| 284 | |
---|
[f613a24] | 285 | |
---|
| 286 | for (int i = 0; i < rows; i++) |
---|
| 287 | { |
---|
| 288 | unsigned x = matrix[i][piv]; |
---|
| 289 | // if the corresponding entry in the matrix is zero, |
---|
| 290 | // there is nothing to do |
---|
| 291 | if (x != 0) |
---|
| 292 | { |
---|
| 293 | for (int j = piv; j < n; j++) |
---|
| 294 | { |
---|
| 295 | if (row[j] != 0) |
---|
| 296 | { |
---|
| 297 | unsigned long tmp = multMod(row[j], x, p); |
---|
| 298 | tmp = p - tmp; |
---|
| 299 | matrix[i][j] += tmp; |
---|
| 300 | if (matrix[i][j] >= p) |
---|
| 301 | { |
---|
| 302 | matrix[i][j] -= p; |
---|
| 303 | } |
---|
| 304 | } |
---|
| 305 | } |
---|
| 306 | } |
---|
| 307 | } |
---|
| 308 | |
---|
[5b2745] | 309 | pivots[rows] = piv; |
---|
[f613a24] | 310 | |
---|
| 311 | // update nonPivots |
---|
| 312 | for (int i = 0; i < n-rows; i++) |
---|
| 313 | { |
---|
| 314 | if (nonPivots[i] == piv) |
---|
| 315 | { |
---|
| 316 | // shift everything one position to the left |
---|
| 317 | for (int j = i; j < n-rows-1; j++) |
---|
| 318 | { |
---|
| 319 | nonPivots[j] = nonPivots[j+1]; |
---|
| 320 | } |
---|
| 321 | |
---|
| 322 | break; |
---|
| 323 | } |
---|
| 324 | } |
---|
| 325 | |
---|
[5b2745] | 326 | rows++; |
---|
| 327 | } |
---|
[2024f6a] | 328 | } |
---|
| 329 | |
---|
| 330 | |
---|
[5b2745] | 331 | void NewVectorMatrix::insertMatrix (LinearDependencyMatrix & mat) |
---|
| 332 | { |
---|
[f613a24] | 333 | for(int i = 0; i < mat.rows; i++) |
---|
| 334 | { |
---|
| 335 | insertRow (mat.matrix[i]); |
---|
| 336 | } |
---|
| 337 | } |
---|
| 338 | |
---|
| 339 | int NewVectorMatrix::findSmallestNonpivot () |
---|
| 340 | { |
---|
| 341 | // This method isn't very efficient, but it is called at most a few times, |
---|
| 342 | // so efficiency is not important. |
---|
| 343 | if(rows == n) |
---|
| 344 | return -1; |
---|
| 345 | |
---|
| 346 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 347 | { |
---|
[f613a24] | 348 | bool isPivot = false; |
---|
| 349 | for(int j = 0; j < rows; j++) |
---|
[5b2745] | 350 | { |
---|
[f613a24] | 351 | if(pivots[j] == i) |
---|
[5b2745] | 352 | { |
---|
[f613a24] | 353 | isPivot = true; |
---|
| 354 | break; |
---|
[5b2745] | 355 | } |
---|
[2024f6a] | 356 | } |
---|
[f613a24] | 357 | |
---|
| 358 | if(!isPivot) |
---|
[5b2745] | 359 | { |
---|
[f613a24] | 360 | return i; |
---|
[5b2745] | 361 | } |
---|
| 362 | } |
---|
[2024f6a] | 363 | } |
---|
| 364 | |
---|
[f613a24] | 365 | int NewVectorMatrix::findLargestNonpivot () |
---|
[5b2745] | 366 | { |
---|
| 367 | // This method isn't very efficient, but it is called at most a few times, so efficiency is not important. |
---|
| 368 | if(rows == n) |
---|
| 369 | return -1; |
---|
[2024f6a] | 370 | |
---|
[f613a24] | 371 | for(int i = n-1; i >= 0; i--) |
---|
[5b2745] | 372 | { |
---|
| 373 | bool isPivot = false; |
---|
[abcc692] | 374 | for(int j = 0; j < rows; j++) |
---|
[5b2745] | 375 | { |
---|
| 376 | if(pivots[j] == i) |
---|
| 377 | { |
---|
| 378 | isPivot = true; |
---|
| 379 | break; |
---|
| 380 | } |
---|
| 381 | } |
---|
| 382 | |
---|
| 383 | if(!isPivot) |
---|
| 384 | { |
---|
| 385 | return i; |
---|
| 386 | } |
---|
| 387 | } |
---|
[056089] | 388 | abort(); |
---|
[2024f6a] | 389 | } |
---|
| 390 | |
---|
| 391 | |
---|
[5b2745] | 392 | void vectorMatrixMult (unsigned long *vec, unsigned long **mat, |
---|
[f613a24] | 393 | unsigned **nonzeroIndices, unsigned *nonzeroCounts, |
---|
[5b2745] | 394 | unsigned long *result, unsigned n, unsigned long p) |
---|
| 395 | { |
---|
| 396 | unsigned long tmp; |
---|
[f613a24] | 397 | |
---|
[abcc692] | 398 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 399 | { |
---|
| 400 | result[i] = 0; |
---|
[f613a24] | 401 | for(int j = 0; j < nonzeroCounts[i]; j++) |
---|
[5b2745] | 402 | { |
---|
[f613a24] | 403 | tmp = multMod (vec[nonzeroIndices[i][j]], mat[nonzeroIndices[i][j]][i], p); |
---|
[5b2745] | 404 | result[i] += tmp; |
---|
[f613a24] | 405 | if (result[i] >= p) |
---|
| 406 | { |
---|
| 407 | result[i] -= p; |
---|
| 408 | } |
---|
[5b2745] | 409 | } |
---|
| 410 | } |
---|
[2024f6a] | 411 | } |
---|
| 412 | |
---|
[abcc692] | 413 | |
---|
[ec9db9] | 414 | #if 0 |
---|
[5b2745] | 415 | void printVec (unsigned long *vec, int n) |
---|
| 416 | { |
---|
| 417 | for(int i = 0; i < n; i++) |
---|
| 418 | { |
---|
| 419 | std::cout << vec[i] << ", "; |
---|
| 420 | } |
---|
| 421 | |
---|
| 422 | std::cout << std::endl; |
---|
[2024f6a] | 423 | } |
---|
[ec9db9] | 424 | #endif |
---|
| 425 | |
---|
[41b5c0] | 426 | |
---|
[5b2745] | 427 | unsigned long *computeMinimalPolynomial (unsigned long **matrix, unsigned n, |
---|
| 428 | unsigned long p) |
---|
| 429 | { |
---|
| 430 | LinearDependencyMatrix lindepmat (n, p); |
---|
| 431 | NewVectorMatrix newvectormat (n, p); |
---|
[2024f6a] | 432 | |
---|
[5b2745] | 433 | unsigned long *result = new unsigned long[n + 1]; |
---|
| 434 | unsigned long *mpvec = new unsigned long[n + 1]; |
---|
| 435 | unsigned long *tmp = new unsigned long[n + 1]; |
---|
[2024f6a] | 436 | |
---|
[41b5c0] | 437 | // initialize result = 1 |
---|
[abcc692] | 438 | for(int i = 0; i <= n; i++) |
---|
[5b2745] | 439 | { |
---|
| 440 | result[i] = 0; |
---|
| 441 | } |
---|
| 442 | result[0] = 1; |
---|
[2024f6a] | 443 | |
---|
[5b2745] | 444 | int degresult = 0; |
---|
[2024f6a] | 445 | |
---|
| 446 | |
---|
[f613a24] | 447 | // Store the indices where the matrix has non-zero entries. |
---|
| 448 | // This has a huge impact on spares matrices. |
---|
| 449 | unsigned* nonzeroCounts = new unsigned[n]; |
---|
| 450 | unsigned** nonzeroIndices = new unsigned*[n]; |
---|
| 451 | for (int i = 0; i < n; i++) |
---|
| 452 | { |
---|
| 453 | nonzeroIndices[i] = new unsigned[n]; |
---|
| 454 | nonzeroCounts[i] = 0; |
---|
| 455 | for (int j = 0; j < n; j++) |
---|
| 456 | { |
---|
| 457 | if (matrix[j][i] != 0) |
---|
| 458 | { |
---|
| 459 | nonzeroIndices[i][nonzeroCounts[i]] = j; |
---|
| 460 | nonzeroCounts[i]++; |
---|
| 461 | } |
---|
| 462 | } |
---|
| 463 | } |
---|
| 464 | |
---|
| 465 | int i = n-1; |
---|
[2024f6a] | 466 | |
---|
[5b2745] | 467 | unsigned long *vec = new unsigned long[n]; |
---|
| 468 | unsigned long *vecnew = new unsigned long[n]; |
---|
[ec9db9] | 469 | |
---|
[f613a24] | 470 | unsigned loopsEven = true; |
---|
[5b2745] | 471 | while(i != -1) |
---|
| 472 | { |
---|
[abcc692] | 473 | for(int j = 0; j < n; j++) |
---|
[5b2745] | 474 | { |
---|
| 475 | vec[j] = 0; |
---|
[2024f6a] | 476 | } |
---|
[5b2745] | 477 | vec[i] = 1; |
---|
[2024f6a] | 478 | |
---|
[5b2745] | 479 | lindepmat.resetMatrix (); |
---|
[2024f6a] | 480 | |
---|
[5b2745] | 481 | while(true) |
---|
| 482 | { |
---|
| 483 | bool ld = lindepmat.findLinearDependency (vec, mpvec); |
---|
[2024f6a] | 484 | |
---|
[5b2745] | 485 | if(ld) |
---|
| 486 | { |
---|
| 487 | break; |
---|
| 488 | } |
---|
[2024f6a] | 489 | |
---|
[f613a24] | 490 | vectorMatrixMult (vec, matrix, nonzeroIndices, nonzeroCounts, vecnew, n, p); |
---|
[5b2745] | 491 | unsigned long *swap = vec; |
---|
| 492 | vec = vecnew; |
---|
| 493 | vecnew = swap; |
---|
| 494 | } |
---|
[2024f6a] | 495 | |
---|
| 496 | |
---|
[5b2745] | 497 | unsigned degmpvec = n; |
---|
| 498 | while(mpvec[degmpvec] == 0) |
---|
| 499 | { |
---|
| 500 | degmpvec--; |
---|
[2024f6a] | 501 | } |
---|
| 502 | |
---|
[41b5c0] | 503 | // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree, |
---|
[5b2745] | 504 | // then we are done. |
---|
| 505 | if(degmpvec == n) |
---|
| 506 | { |
---|
| 507 | unsigned long *swap = result; |
---|
| 508 | result = mpvec; |
---|
| 509 | mpvec = swap; |
---|
| 510 | i = -1; |
---|
| 511 | } |
---|
| 512 | else |
---|
| 513 | { |
---|
| 514 | // otherwise, we have to compute the lcm of mpvec and prior result |
---|
[2024f6a] | 515 | |
---|
[5b2745] | 516 | // tmp = zeropol |
---|
[abcc692] | 517 | for(int j = 0; j <= n; j++) |
---|
[5b2745] | 518 | { |
---|
| 519 | tmp[j] = 0; |
---|
| 520 | } |
---|
| 521 | degresult = lcm (tmp, result, mpvec, p, degresult, degmpvec); |
---|
| 522 | unsigned long *swap = result; |
---|
| 523 | result = tmp; |
---|
| 524 | tmp = swap; |
---|
[2024f6a] | 525 | |
---|
[5b2745] | 526 | if(degresult == n) |
---|
| 527 | { |
---|
| 528 | // minimal polynomial has highest possible degree, so we are done. |
---|
| 529 | i = -1; |
---|
| 530 | } |
---|
| 531 | else |
---|
| 532 | { |
---|
| 533 | newvectormat.insertMatrix (lindepmat); |
---|
[f613a24] | 534 | |
---|
| 535 | // choose new unit vector from the front or the end, alternating |
---|
| 536 | // for each round. If the matrix is the companion matrix for x^n, |
---|
| 537 | // then taking vectors from the end is best. If it is the transpose, |
---|
| 538 | // taking vectors from the front is best. |
---|
| 539 | // This tries to take the middle way |
---|
| 540 | if (loopsEven) |
---|
| 541 | { |
---|
| 542 | i = newvectormat.findSmallestNonpivot (); |
---|
| 543 | } |
---|
| 544 | else |
---|
| 545 | { |
---|
| 546 | i = newvectormat.findLargestNonpivot (); |
---|
| 547 | } |
---|
[5b2745] | 548 | } |
---|
[2024f6a] | 549 | } |
---|
[f613a24] | 550 | |
---|
| 551 | loopsEven = !loopsEven; |
---|
| 552 | } |
---|
| 553 | |
---|
| 554 | for (int i = 0; i < n; i++) |
---|
| 555 | { |
---|
| 556 | delete[] nonzeroIndices[i]; |
---|
[5b2745] | 557 | } |
---|
[f613a24] | 558 | delete[] nonzeroIndices; |
---|
| 559 | delete[] nonzeroCounts; |
---|
[2024f6a] | 560 | |
---|
| 561 | |
---|
[5b2745] | 562 | delete[]vecnew; |
---|
| 563 | delete[]vec; |
---|
| 564 | delete[]tmp; |
---|
| 565 | delete[]mpvec; |
---|
[2024f6a] | 566 | |
---|
[5b2745] | 567 | return result; |
---|
[2024f6a] | 568 | } |
---|
| 569 | |
---|
| 570 | |
---|
[5b2745] | 571 | void rem (unsigned long *a, unsigned long *q, unsigned long p, int °a, |
---|
| 572 | int degq) |
---|
| 573 | { |
---|
| 574 | while(degq <= dega) |
---|
| 575 | { |
---|
| 576 | unsigned d = dega - degq; |
---|
[abcc692] | 577 | long factor = multMod (a[dega], modularInverse (q[degq], p), p); |
---|
[5b2745] | 578 | for(int i = degq; i >= 0; i--) |
---|
| 579 | { |
---|
[abcc692] | 580 | long tmp = p - multMod (factor, q[i], p); |
---|
[5b2745] | 581 | a[d + i] += tmp; |
---|
[f613a24] | 582 | if (a[d + i] >= p) |
---|
| 583 | { |
---|
| 584 | a[d + i] -= p; |
---|
| 585 | } |
---|
[5b2745] | 586 | } |
---|
[2024f6a] | 587 | |
---|
[f613a24] | 588 | while(dega >= 0 && a[dega] == 0) |
---|
[5b2745] | 589 | { |
---|
| 590 | dega--; |
---|
[2024f6a] | 591 | } |
---|
[5b2745] | 592 | } |
---|
[2024f6a] | 593 | } |
---|
| 594 | |
---|
| 595 | |
---|
[5b2745] | 596 | void quo (unsigned long *a, unsigned long *q, unsigned long p, int °a, |
---|
| 597 | int degq) |
---|
| 598 | { |
---|
| 599 | unsigned degres = dega - degq; |
---|
| 600 | unsigned long *result = new unsigned long[degres + 1]; |
---|
| 601 | |
---|
[f613a24] | 602 | // initialize to zero |
---|
| 603 | for (int i = 0; i <= degres; i++) |
---|
| 604 | { |
---|
| 605 | result[i] = 0; |
---|
| 606 | } |
---|
| 607 | |
---|
[5b2745] | 608 | while(degq <= dega) |
---|
| 609 | { |
---|
| 610 | unsigned d = dega - degq; |
---|
[f613a24] | 611 | unsigned long inv = modularInverse (q[degq], p); |
---|
[abcc692] | 612 | result[d] = multMod (a[dega], inv, p); |
---|
| 613 | for(int i = degq; i >= 0; i--) |
---|
[5b2745] | 614 | { |
---|
[f613a24] | 615 | unsigned long tmp = p - multMod (result[d], q[i], p); |
---|
[5b2745] | 616 | a[d + i] += tmp; |
---|
[f613a24] | 617 | if (a[d + i] >= p) |
---|
| 618 | { |
---|
| 619 | a[d + i] -= p; |
---|
| 620 | } |
---|
[5b2745] | 621 | } |
---|
| 622 | |
---|
[f613a24] | 623 | while(dega >= 0 && a[dega] == 0) |
---|
[5b2745] | 624 | { |
---|
| 625 | dega--; |
---|
| 626 | } |
---|
| 627 | } |
---|
| 628 | |
---|
| 629 | // copy result into a |
---|
[abcc692] | 630 | for(int i = 0; i <= degres; i++) |
---|
[5b2745] | 631 | { |
---|
| 632 | a[i] = result[i]; |
---|
| 633 | } |
---|
| 634 | // set all following entries of a to zero |
---|
[abcc692] | 635 | for(int i = degres + 1; i <= degq + degres; i++) |
---|
[5b2745] | 636 | { |
---|
| 637 | a[i] = 0; |
---|
| 638 | } |
---|
| 639 | |
---|
| 640 | dega = degres; |
---|
| 641 | |
---|
| 642 | delete[]result; |
---|
| 643 | } |
---|
[2024f6a] | 644 | |
---|
| 645 | |
---|
[5b2745] | 646 | void mult (unsigned long *result, unsigned long *a, unsigned long *b, |
---|
| 647 | unsigned long p, int dega, int degb) |
---|
| 648 | { |
---|
| 649 | // NOTE: we assume that every entry in result is preinitialized to zero! |
---|
[2024f6a] | 650 | |
---|
[5b2745] | 651 | for(int i = 0; i <= dega; i++) |
---|
| 652 | { |
---|
| 653 | for(int j = 0; j <= degb; j++) |
---|
| 654 | { |
---|
[abcc692] | 655 | result[i + j] += multMod (a[i], b[j], p); |
---|
[f613a24] | 656 | if (result[i + j] >= p) |
---|
| 657 | { |
---|
| 658 | result[i + j] -= p; |
---|
| 659 | } |
---|
[2024f6a] | 660 | } |
---|
[5b2745] | 661 | } |
---|
[2024f6a] | 662 | } |
---|
| 663 | |
---|
| 664 | |
---|
[5b2745] | 665 | int gcd (unsigned long *g, unsigned long *a, unsigned long *b, |
---|
| 666 | unsigned long p, int dega, int degb) |
---|
| 667 | { |
---|
| 668 | unsigned long *tmp1 = new unsigned long[dega + 1]; |
---|
| 669 | unsigned long *tmp2 = new unsigned long[degb + 1]; |
---|
| 670 | for(int i = 0; i <= dega; i++) |
---|
| 671 | { |
---|
| 672 | tmp1[i] = a[i]; |
---|
| 673 | } |
---|
| 674 | for(int i = 0; i <= degb; i++) |
---|
| 675 | { |
---|
| 676 | tmp2[i] = b[i]; |
---|
| 677 | } |
---|
| 678 | int degtmp1 = dega; |
---|
| 679 | int degtmp2 = degb; |
---|
| 680 | |
---|
| 681 | unsigned long *swappol; |
---|
| 682 | int swapdeg; |
---|
| 683 | |
---|
| 684 | while(degtmp2 >= 0) |
---|
| 685 | { |
---|
| 686 | rem (tmp1, tmp2, p, degtmp1, degtmp2); |
---|
| 687 | swappol = tmp1; |
---|
| 688 | tmp1 = tmp2; |
---|
| 689 | tmp2 = swappol; |
---|
| 690 | |
---|
| 691 | swapdeg = degtmp1; |
---|
| 692 | degtmp1 = degtmp2; |
---|
| 693 | degtmp2 = swapdeg; |
---|
| 694 | } |
---|
| 695 | |
---|
| 696 | for(int i = 0; i <= degtmp1; i++) |
---|
| 697 | { |
---|
| 698 | g[i] = tmp1[i]; |
---|
| 699 | } |
---|
| 700 | |
---|
| 701 | delete[]tmp1; |
---|
| 702 | delete[]tmp2; |
---|
| 703 | |
---|
| 704 | return degtmp1; |
---|
| 705 | } |
---|
[2024f6a] | 706 | |
---|
| 707 | |
---|
[5b2745] | 708 | int lcm (unsigned long *l, unsigned long *a, unsigned long *b, |
---|
| 709 | unsigned long p, int dega, int degb) |
---|
| 710 | { |
---|
| 711 | unsigned long *g = new unsigned long[dega + 1]; |
---|
| 712 | // initialize entries of g to zero |
---|
| 713 | for(int i = 0; i <= dega; i++) |
---|
| 714 | { |
---|
| 715 | g[i] = 0; |
---|
| 716 | } |
---|
| 717 | |
---|
| 718 | int degg = gcd (g, a, b, p, dega, degb); |
---|
| 719 | |
---|
| 720 | if(degg > 0) |
---|
| 721 | { |
---|
| 722 | // non-trivial gcd, so compute a = (a/g) |
---|
| 723 | quo (a, g, p, dega, degg); |
---|
| 724 | } |
---|
| 725 | mult (l, a, b, p, dega, degb); |
---|
| 726 | |
---|
| 727 | // normalize |
---|
| 728 | if(l[dega + degb + 1] != 1) |
---|
| 729 | { |
---|
| 730 | unsigned long inv = modularInverse (l[dega + degb], p); |
---|
| 731 | for(int i = 0; i <= dega + degb; i++) |
---|
| 732 | { |
---|
[abcc692] | 733 | l[i] = multMod (l[i], inv, p); |
---|
[5b2745] | 734 | } |
---|
| 735 | } |
---|
| 736 | |
---|
| 737 | return dega + degb; |
---|
[2024f6a] | 738 | } |
---|
| 739 | |
---|
| 740 | |
---|
[abcc692] | 741 | // computes the gcd and the cofactors of u and v |
---|
| 742 | // gcd(u,v) = u3 = u*u1 + v*u2 |
---|
| 743 | unsigned long modularInverse (long long x, long long p) |
---|
[5b2745] | 744 | { |
---|
[abcc692] | 745 | long long u1 = 1; |
---|
| 746 | long long u2 = 0; |
---|
| 747 | long long u3 = x; |
---|
| 748 | long long v1 = 0; |
---|
| 749 | long long v2 = 1; |
---|
| 750 | long long v3 = p; |
---|
| 751 | |
---|
| 752 | long long q, t1, t2, t3; |
---|
[5b2745] | 753 | while(v3 != 0) |
---|
| 754 | { |
---|
| 755 | q = u3 / v3; |
---|
| 756 | t1 = u1 - q * v1; |
---|
| 757 | t2 = u2 - q * v2; |
---|
| 758 | t3 = u3 - q * v3; |
---|
| 759 | u1 = v1; |
---|
| 760 | u2 = v2; |
---|
| 761 | u3 = v3; |
---|
| 762 | v1 = t1; |
---|
| 763 | v2 = t2; |
---|
| 764 | v3 = t3; |
---|
| 765 | } |
---|
| 766 | |
---|
| 767 | if(u1 < 0) |
---|
| 768 | { |
---|
| 769 | u1 += p; |
---|
| 770 | } |
---|
| 771 | |
---|
[abcc692] | 772 | return (unsigned long) u1; |
---|
[2024f6a] | 773 | } |
---|
[41b5c0] | 774 | |
---|