[2024f6a] | 1 | /*********************************************** |
---|
| 2 | * Copyright (C) 2011 Sebastian Jambor * |
---|
| 3 | * sebastian@momo.math.rwth-aachen.de * |
---|
| 4 | ***********************************************/ |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | #include<cmath> |
---|
[f5b40a] | 8 | #include<kernel/mod2.h> |
---|
[41b5c0] | 9 | //#include<iomanip> |
---|
[2024f6a] | 10 | |
---|
| 11 | #include "minpoly.h" |
---|
| 12 | |
---|
| 13 | // TODO: avoid code copying, use subclassing instead! |
---|
| 14 | |
---|
[5b2745] | 15 | LinearDependencyMatrix::LinearDependencyMatrix (unsigned n, unsigned long p) |
---|
| 16 | { |
---|
| 17 | this->n = n; |
---|
| 18 | this->p = p; |
---|
| 19 | |
---|
| 20 | matrix = new unsigned long *[n]; |
---|
[abcc692] | 21 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 22 | { |
---|
| 23 | matrix[i] = new unsigned long[2 * n + 1]; |
---|
| 24 | } |
---|
| 25 | pivots = new unsigned[n]; |
---|
| 26 | tmprow = new unsigned long[2 * n + 1]; |
---|
| 27 | rows = 0; |
---|
[2024f6a] | 28 | } |
---|
| 29 | |
---|
[5b2745] | 30 | LinearDependencyMatrix::~LinearDependencyMatrix () |
---|
| 31 | { |
---|
| 32 | delete[]tmprow; |
---|
| 33 | delete[]pivots; |
---|
[2024f6a] | 34 | |
---|
[abcc692] | 35 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 36 | { |
---|
| 37 | delete[](matrix[i]); |
---|
| 38 | } |
---|
| 39 | delete[]matrix; |
---|
[2024f6a] | 40 | } |
---|
| 41 | |
---|
[5b2745] | 42 | void LinearDependencyMatrix::resetMatrix () |
---|
| 43 | { |
---|
| 44 | rows = 0; |
---|
[2024f6a] | 45 | } |
---|
| 46 | |
---|
[5b2745] | 47 | int LinearDependencyMatrix::firstNonzeroEntry (unsigned long *row) |
---|
| 48 | { |
---|
[abcc692] | 49 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 50 | if(row[i] != 0) |
---|
| 51 | return i; |
---|
[2024f6a] | 52 | |
---|
[5b2745] | 53 | return -1; |
---|
[2024f6a] | 54 | } |
---|
| 55 | |
---|
[5b2745] | 56 | void LinearDependencyMatrix::reduceTmpRow () |
---|
| 57 | { |
---|
[abcc692] | 58 | for(int i = 0; i < rows; i++) |
---|
[5b2745] | 59 | { |
---|
| 60 | unsigned piv = pivots[i]; |
---|
[abcc692] | 61 | unsigned x = tmprow[piv]; |
---|
[5b2745] | 62 | // if the corresponding entry in the row is zero, there is nothing to do |
---|
| 63 | if(x != 0) |
---|
| 64 | { |
---|
| 65 | // subtract tmprow[i] times the i-th row |
---|
[abcc692] | 66 | for(int j = piv; j < n + rows + 1; j++) |
---|
[5b2745] | 67 | { |
---|
[abcc692] | 68 | unsigned long tmp = multMod (matrix[i][j], x, p); |
---|
[5b2745] | 69 | tmp = p - tmp; |
---|
| 70 | tmprow[j] += tmp; |
---|
[abcc692] | 71 | tmprow[j] %= p; |
---|
[5b2745] | 72 | } |
---|
| 73 | } |
---|
| 74 | } |
---|
[2024f6a] | 75 | } |
---|
| 76 | |
---|
| 77 | |
---|
[5b2745] | 78 | void LinearDependencyMatrix::normalizeTmp (unsigned i) |
---|
| 79 | { |
---|
| 80 | unsigned long inv = modularInverse (tmprow[i], p); |
---|
| 81 | tmprow[i] = 1; |
---|
[abcc692] | 82 | for(int j = i + 1; j < 2 * n + 1; j++) |
---|
| 83 | tmprow[j] = multMod (tmprow[j], inv, p); |
---|
[2024f6a] | 84 | } |
---|
[ec9db9] | 85 | |
---|
[5b2745] | 86 | bool LinearDependencyMatrix::findLinearDependency (unsigned long *newRow, |
---|
| 87 | unsigned long *dep) |
---|
| 88 | { |
---|
| 89 | // Copy newRow to tmprow (we need to add RHS) |
---|
[abcc692] | 90 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 91 | { |
---|
| 92 | tmprow[i] = newRow[i]; |
---|
| 93 | tmprow[n + i] = 0; |
---|
| 94 | } |
---|
| 95 | tmprow[2 * n] = 0; |
---|
| 96 | tmprow[n + rows] = 1; |
---|
| 97 | |
---|
| 98 | reduceTmpRow (); |
---|
| 99 | |
---|
| 100 | // Is tmprow reduced to zero? Then we have found a linear dependence. |
---|
| 101 | // Otherwise, we just add tmprow to the matrix. |
---|
[abcc692] | 102 | unsigned newpivot = firstNonzeroEntry (tmprow); |
---|
[5b2745] | 103 | if(newpivot == -1) |
---|
| 104 | { |
---|
[abcc692] | 105 | for(int i = 0; i <= n; i++) |
---|
[5b2745] | 106 | { |
---|
| 107 | dep[i] = tmprow[n + i]; |
---|
| 108 | } |
---|
| 109 | |
---|
| 110 | return true; |
---|
| 111 | } |
---|
| 112 | else |
---|
| 113 | { |
---|
| 114 | normalizeTmp (newpivot); |
---|
| 115 | |
---|
[abcc692] | 116 | for(int i = 0; i < 2 * n + 1; i++) |
---|
[5b2745] | 117 | { |
---|
| 118 | matrix[rows][i] = tmprow[i]; |
---|
| 119 | } |
---|
| 120 | |
---|
| 121 | pivots[rows] = newpivot; |
---|
| 122 | rows++; |
---|
| 123 | |
---|
| 124 | return false; |
---|
| 125 | } |
---|
[2024f6a] | 126 | } |
---|
| 127 | |
---|
[ec9db9] | 128 | #if 0 |
---|
[5b2745] | 129 | std::ostream & operator<< (std::ostream & out, |
---|
| 130 | const LinearDependencyMatrix & mat) |
---|
| 131 | { |
---|
| 132 | int width = ((int) log10 (mat.p)) + 1; |
---|
[2024f6a] | 133 | |
---|
[abcc692] | 134 | out << "Pivots: " << std::endl; |
---|
| 135 | for(int j = 0; j < mat.n; j++) |
---|
| 136 | { |
---|
| 137 | out << std::setw (width) << mat.pivots[j] << " "; |
---|
| 138 | } |
---|
| 139 | out << std::endl; |
---|
| 140 | out << "matrix:" << std::endl; |
---|
[5b2745] | 141 | for(int i = 0; i < mat.rows; i++) |
---|
| 142 | { |
---|
| 143 | for(int j = 0; j < mat.n; j++) |
---|
| 144 | { |
---|
| 145 | out << std::setw (width) << mat.matrix[i][j] << " "; |
---|
[2024f6a] | 146 | } |
---|
| 147 | out << "| "; |
---|
[5b2745] | 148 | for(int j = mat.n; j <= 2 * mat.n; j++) |
---|
| 149 | { |
---|
| 150 | out << std::setw (width) << mat.matrix[i][j] << " "; |
---|
[2024f6a] | 151 | } |
---|
| 152 | out << std::endl; |
---|
[5b2745] | 153 | } |
---|
[abcc692] | 154 | out << "tmprow: " << std::endl; |
---|
[5b2745] | 155 | for(int j = 0; j < mat.n; j++) |
---|
| 156 | { |
---|
| 157 | out << std::setw (width) << mat.tmprow[j] << " "; |
---|
| 158 | } |
---|
| 159 | out << "| "; |
---|
| 160 | for(int j = mat.n; j <= 2 * mat.n; j++) |
---|
| 161 | { |
---|
| 162 | out << std::setw (width) << mat.tmprow[j] << " "; |
---|
| 163 | } |
---|
| 164 | out << std::endl; |
---|
| 165 | |
---|
| 166 | return out; |
---|
[2024f6a] | 167 | } |
---|
[ec9db9] | 168 | #endif |
---|
[2024f6a] | 169 | |
---|
| 170 | |
---|
[5b2745] | 171 | NewVectorMatrix::NewVectorMatrix (unsigned n, unsigned long p) |
---|
| 172 | { |
---|
| 173 | this->n = n; |
---|
| 174 | this->p = p; |
---|
[2024f6a] | 175 | |
---|
[5b2745] | 176 | matrix = new unsigned long *[n]; |
---|
[abcc692] | 177 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 178 | { |
---|
| 179 | matrix[i] = new unsigned long[n]; |
---|
| 180 | } |
---|
[2024f6a] | 181 | |
---|
[5b2745] | 182 | pivots = new unsigned[n]; |
---|
| 183 | rows = 0; |
---|
[2024f6a] | 184 | } |
---|
| 185 | |
---|
[5b2745] | 186 | NewVectorMatrix::~NewVectorMatrix () |
---|
| 187 | { |
---|
| 188 | delete pivots; |
---|
[2024f6a] | 189 | |
---|
[abcc692] | 190 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 191 | { |
---|
| 192 | delete[]matrix[i]; |
---|
| 193 | } |
---|
| 194 | delete matrix; |
---|
[2024f6a] | 195 | } |
---|
| 196 | |
---|
[5b2745] | 197 | int NewVectorMatrix::firstNonzeroEntry (unsigned long *row) |
---|
| 198 | { |
---|
[abcc692] | 199 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 200 | if(row[i] != 0) |
---|
| 201 | return i; |
---|
[2024f6a] | 202 | |
---|
[5b2745] | 203 | return -1; |
---|
[2024f6a] | 204 | } |
---|
| 205 | |
---|
[5b2745] | 206 | void NewVectorMatrix::normalizeRow (unsigned long *row, unsigned i) |
---|
| 207 | { |
---|
| 208 | unsigned long inv = modularInverse (row[i], p); |
---|
| 209 | row[i] = 1; |
---|
| 210 | |
---|
[abcc692] | 211 | for(int j = i + 1; j < n; j++) |
---|
[5b2745] | 212 | { |
---|
[abcc692] | 213 | row[j] = multMod (row[j], inv, p); |
---|
[5b2745] | 214 | } |
---|
[2024f6a] | 215 | } |
---|
| 216 | |
---|
[5b2745] | 217 | void NewVectorMatrix::insertRow (unsigned long *row) |
---|
| 218 | { |
---|
[abcc692] | 219 | for(int i = 0; i < rows; i++) |
---|
[5b2745] | 220 | { |
---|
| 221 | unsigned piv = pivots[i]; |
---|
[abcc692] | 222 | unsigned x = row[piv]; |
---|
[5b2745] | 223 | // if the corresponding entry in the row is zero, there is nothing to do |
---|
| 224 | if(x != 0) |
---|
| 225 | { |
---|
| 226 | // subtract row[i] times the i-th row |
---|
[abcc692] | 227 | for(int j = piv; j < n; j++) |
---|
[5b2745] | 228 | { |
---|
[abcc692] | 229 | unsigned long tmp = multMod (matrix[i][j], x, p); |
---|
[5b2745] | 230 | tmp = p - tmp; |
---|
| 231 | row[j] += tmp; |
---|
[abcc692] | 232 | row[j] %= p; |
---|
[5b2745] | 233 | } |
---|
| 234 | } |
---|
| 235 | } |
---|
| 236 | |
---|
[abcc692] | 237 | unsigned piv = firstNonzeroEntry (row); |
---|
[5b2745] | 238 | |
---|
| 239 | if(piv != -1) |
---|
| 240 | { |
---|
| 241 | // normalize and insert row into the matrix |
---|
| 242 | normalizeRow (row, piv); |
---|
[abcc692] | 243 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 244 | { |
---|
| 245 | matrix[rows][i] = row[i]; |
---|
| 246 | } |
---|
| 247 | |
---|
| 248 | pivots[rows] = piv; |
---|
| 249 | rows++; |
---|
| 250 | } |
---|
[2024f6a] | 251 | } |
---|
| 252 | |
---|
| 253 | |
---|
[5b2745] | 254 | void NewVectorMatrix::insertMatrix (LinearDependencyMatrix & mat) |
---|
| 255 | { |
---|
| 256 | // The matrix in LinearDependencyMatrix is already in reduced form. |
---|
| 257 | // Thus, if the current matrix is empty, we can simply copy the other matrix. |
---|
| 258 | if(rows == 0) |
---|
| 259 | { |
---|
[abcc692] | 260 | for(int i = 0; i < mat.rows; i++) |
---|
[5b2745] | 261 | { |
---|
[abcc692] | 262 | for(int j = 0; j < n; j++) |
---|
[5b2745] | 263 | { |
---|
| 264 | matrix[i][j] = mat.matrix[i][j]; |
---|
| 265 | |
---|
| 266 | rows = mat.rows; |
---|
[abcc692] | 267 | for(int i = 0; i < rows; i++) |
---|
[5b2745] | 268 | { |
---|
[abcc692] | 269 | pivots[i] = mat.pivots[i]; |
---|
[2024f6a] | 270 | } |
---|
[5b2745] | 271 | } |
---|
[2024f6a] | 272 | } |
---|
[5b2745] | 273 | } |
---|
| 274 | else |
---|
| 275 | { |
---|
[abcc692] | 276 | for(int i = 0; i < mat.rows; i++) |
---|
[5b2745] | 277 | { |
---|
| 278 | insertRow (mat.matrix[i]); |
---|
| 279 | } |
---|
| 280 | } |
---|
[2024f6a] | 281 | } |
---|
| 282 | |
---|
[5b2745] | 283 | int NewVectorMatrix::findSmallestNonpivot () |
---|
| 284 | { |
---|
| 285 | // This method isn't very efficient, but it is called at most a few times, so efficiency is not important. |
---|
| 286 | if(rows == n) |
---|
| 287 | return -1; |
---|
[2024f6a] | 288 | |
---|
[abcc692] | 289 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 290 | { |
---|
| 291 | bool isPivot = false; |
---|
[abcc692] | 292 | for(int j = 0; j < rows; j++) |
---|
[5b2745] | 293 | { |
---|
| 294 | if(pivots[j] == i) |
---|
| 295 | { |
---|
| 296 | isPivot = true; |
---|
| 297 | break; |
---|
| 298 | } |
---|
| 299 | } |
---|
| 300 | |
---|
| 301 | if(!isPivot) |
---|
| 302 | { |
---|
| 303 | return i; |
---|
| 304 | } |
---|
| 305 | } |
---|
[2024f6a] | 306 | } |
---|
| 307 | |
---|
| 308 | |
---|
[5b2745] | 309 | void vectorMatrixMult (unsigned long *vec, unsigned long **mat, |
---|
| 310 | unsigned long *result, unsigned n, unsigned long p) |
---|
| 311 | { |
---|
| 312 | unsigned long tmp; |
---|
[abcc692] | 313 | for(int i = 0; i < n; i++) |
---|
[5b2745] | 314 | { |
---|
| 315 | result[i] = 0; |
---|
[abcc692] | 316 | for(int j = 0; j < n; j++) |
---|
[5b2745] | 317 | { |
---|
[abcc692] | 318 | tmp = multMod (vec[j], mat[j][i], p); |
---|
[5b2745] | 319 | result[i] += tmp; |
---|
[abcc692] | 320 | result[i] %= p; |
---|
[5b2745] | 321 | } |
---|
| 322 | } |
---|
[2024f6a] | 323 | } |
---|
| 324 | |
---|
[abcc692] | 325 | |
---|
[ec9db9] | 326 | #if 0 |
---|
[5b2745] | 327 | void printVec (unsigned long *vec, int n) |
---|
| 328 | { |
---|
| 329 | for(int i = 0; i < n; i++) |
---|
| 330 | { |
---|
| 331 | std::cout << vec[i] << ", "; |
---|
| 332 | } |
---|
| 333 | |
---|
| 334 | std::cout << std::endl; |
---|
[2024f6a] | 335 | } |
---|
[ec9db9] | 336 | #endif |
---|
| 337 | |
---|
[41b5c0] | 338 | |
---|
[5b2745] | 339 | unsigned long *computeMinimalPolynomial (unsigned long **matrix, unsigned n, |
---|
| 340 | unsigned long p) |
---|
| 341 | { |
---|
| 342 | LinearDependencyMatrix lindepmat (n, p); |
---|
| 343 | NewVectorMatrix newvectormat (n, p); |
---|
[2024f6a] | 344 | |
---|
[5b2745] | 345 | unsigned long *result = new unsigned long[n + 1]; |
---|
| 346 | unsigned long *mpvec = new unsigned long[n + 1]; |
---|
| 347 | unsigned long *tmp = new unsigned long[n + 1]; |
---|
[2024f6a] | 348 | |
---|
[41b5c0] | 349 | // initialize result = 1 |
---|
[abcc692] | 350 | for(int i = 0; i <= n; i++) |
---|
[5b2745] | 351 | { |
---|
| 352 | result[i] = 0; |
---|
| 353 | } |
---|
| 354 | result[0] = 1; |
---|
[2024f6a] | 355 | |
---|
[5b2745] | 356 | int degresult = 0; |
---|
[2024f6a] | 357 | |
---|
| 358 | |
---|
[5b2745] | 359 | int i = 0; |
---|
[2024f6a] | 360 | |
---|
[5b2745] | 361 | unsigned long *vec = new unsigned long[n]; |
---|
| 362 | unsigned long *vecnew = new unsigned long[n]; |
---|
[ec9db9] | 363 | |
---|
[5b2745] | 364 | while(i != -1) |
---|
| 365 | { |
---|
[abcc692] | 366 | for(int j = 0; j < n; j++) |
---|
[5b2745] | 367 | { |
---|
| 368 | vec[j] = 0; |
---|
[2024f6a] | 369 | } |
---|
[5b2745] | 370 | vec[i] = 1; |
---|
[2024f6a] | 371 | |
---|
[5b2745] | 372 | lindepmat.resetMatrix (); |
---|
[2024f6a] | 373 | |
---|
[5b2745] | 374 | while(true) |
---|
| 375 | { |
---|
| 376 | bool ld = lindepmat.findLinearDependency (vec, mpvec); |
---|
[2024f6a] | 377 | |
---|
[5b2745] | 378 | if(ld) |
---|
| 379 | { |
---|
| 380 | break; |
---|
| 381 | } |
---|
[2024f6a] | 382 | |
---|
[5b2745] | 383 | vectorMatrixMult (vec, matrix, vecnew, n, p); |
---|
| 384 | unsigned long *swap = vec; |
---|
| 385 | vec = vecnew; |
---|
| 386 | vecnew = swap; |
---|
| 387 | } |
---|
[2024f6a] | 388 | |
---|
| 389 | |
---|
[5b2745] | 390 | unsigned degmpvec = n; |
---|
| 391 | while(mpvec[degmpvec] == 0) |
---|
| 392 | { |
---|
| 393 | degmpvec--; |
---|
[2024f6a] | 394 | } |
---|
| 395 | |
---|
[41b5c0] | 396 | // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree, |
---|
[5b2745] | 397 | // then we are done. |
---|
| 398 | if(degmpvec == n) |
---|
| 399 | { |
---|
| 400 | unsigned long *swap = result; |
---|
| 401 | result = mpvec; |
---|
| 402 | mpvec = swap; |
---|
| 403 | i = -1; |
---|
| 404 | } |
---|
| 405 | else |
---|
| 406 | { |
---|
| 407 | // otherwise, we have to compute the lcm of mpvec and prior result |
---|
[2024f6a] | 408 | |
---|
[5b2745] | 409 | // tmp = zeropol |
---|
[abcc692] | 410 | for(int j = 0; j <= n; j++) |
---|
[5b2745] | 411 | { |
---|
| 412 | tmp[j] = 0; |
---|
| 413 | } |
---|
| 414 | degresult = lcm (tmp, result, mpvec, p, degresult, degmpvec); |
---|
| 415 | unsigned long *swap = result; |
---|
| 416 | result = tmp; |
---|
| 417 | tmp = swap; |
---|
[2024f6a] | 418 | |
---|
[5b2745] | 419 | if(degresult == n) |
---|
| 420 | { |
---|
| 421 | // minimal polynomial has highest possible degree, so we are done. |
---|
| 422 | i = -1; |
---|
| 423 | } |
---|
| 424 | else |
---|
| 425 | { |
---|
| 426 | newvectormat.insertMatrix (lindepmat); |
---|
| 427 | i = newvectormat.findSmallestNonpivot (); |
---|
| 428 | } |
---|
[2024f6a] | 429 | } |
---|
[5b2745] | 430 | } |
---|
[2024f6a] | 431 | |
---|
[5b2745] | 432 | // TODO: take lcms of the different monomials! |
---|
[2024f6a] | 433 | |
---|
[5b2745] | 434 | delete[]vecnew; |
---|
| 435 | delete[]vec; |
---|
| 436 | delete[]tmp; |
---|
| 437 | delete[]mpvec; |
---|
[2024f6a] | 438 | |
---|
[5b2745] | 439 | return result; |
---|
[2024f6a] | 440 | } |
---|
| 441 | |
---|
| 442 | |
---|
[5b2745] | 443 | void rem (unsigned long *a, unsigned long *q, unsigned long p, int °a, |
---|
| 444 | int degq) |
---|
| 445 | { |
---|
| 446 | while(degq <= dega) |
---|
| 447 | { |
---|
| 448 | unsigned d = dega - degq; |
---|
[abcc692] | 449 | long factor = multMod (a[dega], modularInverse (q[degq], p), p); |
---|
[5b2745] | 450 | for(int i = degq; i >= 0; i--) |
---|
| 451 | { |
---|
[abcc692] | 452 | long tmp = p - multMod (factor, q[i], p); |
---|
[5b2745] | 453 | a[d + i] += tmp; |
---|
| 454 | a[d + i] %= p; |
---|
| 455 | } |
---|
[2024f6a] | 456 | |
---|
[5b2745] | 457 | while(a[dega] == 0 && dega >= 0) |
---|
| 458 | { |
---|
| 459 | dega--; |
---|
[2024f6a] | 460 | } |
---|
[5b2745] | 461 | } |
---|
[2024f6a] | 462 | } |
---|
| 463 | |
---|
| 464 | |
---|
[5b2745] | 465 | void quo (unsigned long *a, unsigned long *q, unsigned long p, int °a, |
---|
| 466 | int degq) |
---|
| 467 | { |
---|
| 468 | unsigned degres = dega - degq; |
---|
| 469 | unsigned long *result = new unsigned long[degres + 1]; |
---|
| 470 | |
---|
| 471 | while(degq <= dega) |
---|
| 472 | { |
---|
| 473 | unsigned d = dega - degq; |
---|
| 474 | long inv = modularInverse (q[degq], p); |
---|
[abcc692] | 475 | result[d] = multMod (a[dega], inv, p); |
---|
| 476 | for(int i = degq; i >= 0; i--) |
---|
[5b2745] | 477 | { |
---|
[abcc692] | 478 | long tmp = p - multMod (result[d], q[i], p); |
---|
[5b2745] | 479 | a[d + i] += tmp; |
---|
| 480 | a[d + i] %= p; |
---|
| 481 | } |
---|
| 482 | |
---|
| 483 | while(a[dega] == 0 && dega >= 0) |
---|
| 484 | { |
---|
| 485 | dega--; |
---|
| 486 | } |
---|
| 487 | } |
---|
| 488 | |
---|
| 489 | // copy result into a |
---|
[abcc692] | 490 | for(int i = 0; i <= degres; i++) |
---|
[5b2745] | 491 | { |
---|
| 492 | a[i] = result[i]; |
---|
| 493 | } |
---|
| 494 | // set all following entries of a to zero |
---|
[abcc692] | 495 | for(int i = degres + 1; i <= degq + degres; i++) |
---|
[5b2745] | 496 | { |
---|
| 497 | a[i] = 0; |
---|
| 498 | } |
---|
| 499 | |
---|
| 500 | dega = degres; |
---|
| 501 | |
---|
| 502 | delete[]result; |
---|
| 503 | } |
---|
[2024f6a] | 504 | |
---|
| 505 | |
---|
[5b2745] | 506 | void mult (unsigned long *result, unsigned long *a, unsigned long *b, |
---|
| 507 | unsigned long p, int dega, int degb) |
---|
| 508 | { |
---|
| 509 | // NOTE: we assume that every entry in result is preinitialized to zero! |
---|
[2024f6a] | 510 | |
---|
[5b2745] | 511 | for(int i = 0; i <= dega; i++) |
---|
| 512 | { |
---|
| 513 | for(int j = 0; j <= degb; j++) |
---|
| 514 | { |
---|
[abcc692] | 515 | result[i + j] += multMod (a[i], b[j], p); |
---|
[5b2745] | 516 | result[i + j] %= p; |
---|
[2024f6a] | 517 | } |
---|
[5b2745] | 518 | } |
---|
[2024f6a] | 519 | } |
---|
| 520 | |
---|
| 521 | |
---|
[5b2745] | 522 | int gcd (unsigned long *g, unsigned long *a, unsigned long *b, |
---|
| 523 | unsigned long p, int dega, int degb) |
---|
| 524 | { |
---|
| 525 | unsigned long *tmp1 = new unsigned long[dega + 1]; |
---|
| 526 | unsigned long *tmp2 = new unsigned long[degb + 1]; |
---|
| 527 | for(int i = 0; i <= dega; i++) |
---|
| 528 | { |
---|
| 529 | tmp1[i] = a[i]; |
---|
| 530 | } |
---|
| 531 | for(int i = 0; i <= degb; i++) |
---|
| 532 | { |
---|
| 533 | tmp2[i] = b[i]; |
---|
| 534 | } |
---|
| 535 | int degtmp1 = dega; |
---|
| 536 | int degtmp2 = degb; |
---|
| 537 | |
---|
| 538 | unsigned long *swappol; |
---|
| 539 | int swapdeg; |
---|
| 540 | |
---|
| 541 | while(degtmp2 >= 0) |
---|
| 542 | { |
---|
| 543 | rem (tmp1, tmp2, p, degtmp1, degtmp2); |
---|
| 544 | swappol = tmp1; |
---|
| 545 | tmp1 = tmp2; |
---|
| 546 | tmp2 = swappol; |
---|
| 547 | |
---|
| 548 | swapdeg = degtmp1; |
---|
| 549 | degtmp1 = degtmp2; |
---|
| 550 | degtmp2 = swapdeg; |
---|
| 551 | } |
---|
| 552 | |
---|
| 553 | for(int i = 0; i <= degtmp1; i++) |
---|
| 554 | { |
---|
| 555 | g[i] = tmp1[i]; |
---|
| 556 | } |
---|
| 557 | |
---|
| 558 | delete[]tmp1; |
---|
| 559 | delete[]tmp2; |
---|
| 560 | |
---|
| 561 | return degtmp1; |
---|
| 562 | } |
---|
[2024f6a] | 563 | |
---|
| 564 | |
---|
[5b2745] | 565 | int lcm (unsigned long *l, unsigned long *a, unsigned long *b, |
---|
| 566 | unsigned long p, int dega, int degb) |
---|
| 567 | { |
---|
| 568 | unsigned long *g = new unsigned long[dega + 1]; |
---|
| 569 | // initialize entries of g to zero |
---|
| 570 | for(int i = 0; i <= dega; i++) |
---|
| 571 | { |
---|
| 572 | g[i] = 0; |
---|
| 573 | } |
---|
| 574 | |
---|
| 575 | int degg = gcd (g, a, b, p, dega, degb); |
---|
| 576 | |
---|
| 577 | if(degg > 0) |
---|
| 578 | { |
---|
| 579 | // non-trivial gcd, so compute a = (a/g) |
---|
| 580 | quo (a, g, p, dega, degg); |
---|
| 581 | } |
---|
| 582 | mult (l, a, b, p, dega, degb); |
---|
| 583 | |
---|
| 584 | // normalize |
---|
| 585 | if(l[dega + degb + 1] != 1) |
---|
| 586 | { |
---|
| 587 | unsigned long inv = modularInverse (l[dega + degb], p); |
---|
| 588 | for(int i = 0; i <= dega + degb; i++) |
---|
| 589 | { |
---|
[abcc692] | 590 | l[i] = multMod (l[i], inv, p); |
---|
[5b2745] | 591 | } |
---|
| 592 | } |
---|
| 593 | |
---|
| 594 | return dega + degb; |
---|
[2024f6a] | 595 | } |
---|
| 596 | |
---|
| 597 | |
---|
[abcc692] | 598 | // computes the gcd and the cofactors of u and v |
---|
| 599 | // gcd(u,v) = u3 = u*u1 + v*u2 |
---|
| 600 | unsigned long modularInverse (long long x, long long p) |
---|
[5b2745] | 601 | { |
---|
[abcc692] | 602 | long long u1 = 1; |
---|
| 603 | long long u2 = 0; |
---|
| 604 | long long u3 = x; |
---|
| 605 | long long v1 = 0; |
---|
| 606 | long long v2 = 1; |
---|
| 607 | long long v3 = p; |
---|
| 608 | |
---|
| 609 | long long q, t1, t2, t3; |
---|
[5b2745] | 610 | while(v3 != 0) |
---|
| 611 | { |
---|
| 612 | q = u3 / v3; |
---|
| 613 | t1 = u1 - q * v1; |
---|
| 614 | t2 = u2 - q * v2; |
---|
| 615 | t3 = u3 - q * v3; |
---|
| 616 | u1 = v1; |
---|
| 617 | u2 = v2; |
---|
| 618 | u3 = v3; |
---|
| 619 | v1 = t1; |
---|
| 620 | v2 = t2; |
---|
| 621 | v3 = t3; |
---|
| 622 | } |
---|
| 623 | |
---|
| 624 | if(u1 < 0) |
---|
| 625 | { |
---|
| 626 | u1 += p; |
---|
| 627 | } |
---|
| 628 | |
---|
[abcc692] | 629 | return (unsigned long) u1; |
---|
[2024f6a] | 630 | } |
---|
[41b5c0] | 631 | |
---|