1 2This file explains the functions of the library routines. 3 4Note: 5 Type FiniteField is defined as unsigned long 6 Type Double is defined as double 7 8 9extern long nullspaceLong(const long n, 10 const long m, 11 const long *A, 12 mpz_t * *mp_N_pass); 13 14/* 15 * Calling Sequence: 16 * nullspaceLong(n, m, A, mp_N_pass) 17 * 18 * Summary: 19 * Compute the right nullspace of A. 20 * 21 * Input: n: long, row dimension of A 22 * m: long, column dimension of A 23 * A: 1-dim signed long array length n*m, representing n x m matrix 24 * in row major order 25 * 26 * Output: 27 * - *mp_N_pass: points to a 1-dim mpz_t array of length m*s, where s is the 28 * dimension of the right nullspace of A 29 * - the dimension s of the nullspace is returned 30 * 31 * Notes: 32 * - The matrix A is represented by one-dimension array in row major order. 33 * - Space for what mp_N_points to is allocated by this procedure: if the 34 * nullspace is empty, mp_N_pass is set to NULL. 35 */ 36 37 38 39extern long nullspaceMP(const long n, 40 const long m, 41 const mpz_t *A, 42 mpz_t * *mp_N_pass); 43 44/* 45 * Calling Sequence: 46 * nullspaceMP(n, m, A, mp_N_pass) 47 * 48 * Summary: Compute the right nullspace of A. In this function A is a 49 * 1-dimensional mpz_t array. 50 * 51 * Input: n: long, row dimension of A 52 * m: long, column dimension of A 53 * A: 1-dim mpz_t array length n*m, representing n x m matrix 54 * in row major order 55 * 56 * Output: 57 * - *mp_N_pass: points to a 1-dim mpz_t array of length m*s, where s is the 58 * dimension of the right nullspace of A 59 * - the dimension s of the nullspace is returned 60 * 61 * Notes: 62 * - The matrix A is represented by one-dimension array in row major order. 63 * - Space for what mp_N_points to is allocated by this procedure: if the 64 * nullspace is empty, mp_N_pass is set to NULL. 65 */ 66 67 68 69extern void nonsingSolvMM (const enum SOLU_POS solupos, const long n, 70 const long m, const long *A, mpz_t *mp_B, 71 mpz_t *mp_N, mpz_t mp_D); 72 73/* 74 * Calling Sequence: 75 * nonsingSolvMM(solupos, n, m, A, mp_B, mp_N, mp_D) 76 * 77 * Summary: 78 * Solve nonsingular system of linear equations, where the left hand side 79 * input matrix is a signed long matrix. 80 * 81 * Description: 82 * Given a n x n nonsingular signed long matrix A, a n x m or m x n mpz_t 83 * matrix mp_B, this function will compute the solution of the system 84 * 1. AX = mp_B 85 * 2. XA = mp_B. 86 * The parameter solupos controls whether the system is in the type of 1 87 * or 2. 88 * 89 * Since the unique solution X is a rational matrix, the function will 90 * output the numerator matrix mp_N and denominator mp_D respectively, 91 * such that A(mp_N) = mp_D*mp_B or (mp_N)A = mp_D*mp_B. 92 * 93 * Input: 94 * solupos: enumerate, flag to indicate the system to be solved 95 * - solupos = LeftSolu: solve XA = mp_B 96 * - solupos = RightSolu: solve AX = mp_B 97 * n: long, dimension of A 98 * m: long, column or row dimension of mp_B depending on solupos 99 * A: 1-dim signed long array length n*n, representing the n x n 100 * left hand side input matrix 101 * mp_B: 1-dim mpz_t array length n*m, representing the right hand side 102 * matrix of the system 103 * - solupos = LeftSolu: mp_B a m x n matrix 104 * - solupos = RightSolu: mp_B a n x m matrix 105 * 106 * Output: 107 * mp_N: 1-dim mpz_t array length n*m, representing the numerator matrix 108 * of the solution 109 * - solupos = LeftSolu: mp_N a m x n matrix 110 * - solupos = RightSolu: mp_N a n x m matrix 111 * mp_D: mpz_t, denominator of the solution 112 * 113 * Precondition: 114 * A must be a nonsingular matrix. 115 * 116 * Note: 117 * - It is necessary to make sure the input parameters are correct, 118 * expecially the dimension, since there is no parameter checks in the 119 * function. 120 * - Input and output matrices are row majored and represented by 121 * one-dimension array. 122 * - It is needed to preallocate the memory space of mp_N and mp_D. 123 * 124 */ 125 126 127 128extern void nonsingSolvLlhsMM (const enum SOLU_POS solupos, const long n, 129 const long m, mpz_t *mp_A, mpz_t *mp_B, 130 mpz_t *mp_N, mpz_t mp_D); 131 132/* 133 * Calling Sequence: 134 * nonsingSolvLlhsMM(solupos, n, m, mp_A, mp_B, mp_N, mp_D) 135 * 136 * Summary: 137 * Solve nonsingular system of linear equations, where the left hand side 138 * input matrix is a mpz_t matrix. 139 * 140 * Description: 141 * Given a n x n nonsingular mpz_t matrix A, a n x m or m x n mpz_t 142 * matrix mp_B, this function will compute the solution of the system 143 * 1. (mp_A)X = mp_B 144 * 2. X(mp_A) = mp_B. 145 * The parameter solupos controls whether the system is in the type of 1 146 * or 2. 147 * 148 * Since the unique solution X is a rational matrix, the function will 149 * output the numerator matrix mp_N and denominator mp_D respectively, 150 * such that mp_Amp_N = mp_D*mp_B or mp_Nmp_A = mp_D*mp_B. 151 * 152 * Input: 153 * solupos: enumerate, flag to indicate the system to be solved 154 * - solupos = LeftSolu: solve XA = mp_B 155 * - solupos = RightSolu: solve AX = mp_B 156 * n: long, dimension of A 157 * m: long, column or row dimension of mp_B depending on solupos 158 * mp_A: 1-dim mpz_t array length n*n, representing the n x n left hand 159 * side input matrix 160 * mp_B: 1-dim mpz_t array length n*m, representing the right hand side 161 * matrix of the system 162 * - solupos = LeftSolu: mp_B a m x n matrix 163 * - solupos = RightSolu: mp_B a n x m matrix 164 * 165 * Output: 166 * mp_N: 1-dim mpz_t array length n*m, representing the numerator matrix 167 * of the solution 168 * - solupos = LeftSolu: mp_N a m x n matrix 169 * - solupos = RightSolu: mp_N a n x m matrix 170 * mp_D: mpz_t, denominator of the solution 171 * 172 * Precondition: 173 * mp_A must be a nonsingular matrix. 174 * 175 * Note: 176 * - It is necessary to make sure the input parameters are correct, 177 * expecially the dimension, since there is no parameter checks in the 178 * function. 179 * - Input and output matrices are row majored and represented by 180 * one-dimension array. 181 * - It is needed to preallocate the memory space of mp_N and mp_D. 182 * 183 */ 184 185 186 187extern void nonsingSolvRNSMM (const enum SOLU_POS solupos, const long n, 188 const long m, const long basislen, 189 const FiniteField *basis, Double **ARNS, 190 mpz_t *mp_B, mpz_t *mp_N, mpz_t mp_D); 191 192/* 193 * Calling Sequence: 194 * nonsingSolvRNSMM(solupos, basislen, n, m, basis, ARNS, mp_B, mp_N, mp_D) 195 * 196 * Summary: 197 * Solve nonsingular system of linear equations, where the left hand side 198 * input matrix is represented in a RNS. 199 * 200 * Description: 201 * Given a n x n nonsingular matrix A represented in a RNS, a n x m or m x n 202 * mpz_t matrix mp_B, this function will compute the solution of the system 203 * 1. AX = mp_B 204 * 2. XA = mp_B. 205 * The parameter solupos controls whether the system is in the type of 1 206 * or 2. 207 * 208 * Since the unique solution X is a rational matrix, the function will 209 * output the numerator matrix mp_N and denominator mp_D respectively, 210 * such that A(mp_N) = mp_D*mp_B or (mp_N)A = mp_D*mp_B. 211 * 212 * Input: 213 * solupos: enumerate, flag to indicate the system to be solved 214 * - solupos = LeftSolu: solve XA = mp_B 215 * - solupos = RightSolu: solve AX = mp_B 216 * basislen: long, dimension of RNS basis 217 * n: long, dimension of A 218 * m: long, column or row dimension of mp_B depending on solupos 219 * basis: 1-dim FiniteField array length basislen, RNS basis 220 * ARNS: 2-dim Double array, dimension basislen x n*n, representation of 221 * n x n input matrix A in RNS, where ARNS[i] = A mod basis[i] 222 * mp_B: 1-dim mpz_t array length n*m, representing the right hand side 223 * matrix of the system 224 * - solupos = LeftSolu: mp_B a m x n matrix 225 * - solupos = RightSolu: mp_B a n x m matrix 226 * 227 * Output: 228 * mp_N: 1-dim mpz_t array length n*m, representing the numerator matrix 229 * of the solution 230 * - solupos = LeftSolu: mp_N a m x n matrix 231 * - solupos = RightSolu: mp_N a n x m matrix 232 * mp_D: mpz_t, denominator of the solution 233 * 234 * Precondition: 235 * - A must be a nonsingular matrix. 236 * - Any element p in RNS basis must satisfy 2*(p-1)^2 <= 2^53-1. 237 * 238 * Note: 239 * - It is necessary to make sure the input parameters are correct, 240 * expecially the dimension, since there is no parameter checks in the 241 * function. 242 * - Input and output matrices are row majored and represented by 243 * one-dimension array. 244 * - It is needed to preallocate the memory space of mp_N and mp_D. 245 * 246 */ 247 248 249 250extern long certSolveLong (const long certflag, const long n, const long m, 251 const long *A, mpz_t *mp_b, mpz_t *mp_N, 252 mpz_t mp_D, mpz_t *mp_NZ, mpz_t mp_DZ); 253 254/* 255 * 256 * Calling Sequence: 257 * 1/2/3 <-- certSolveLong(certflag, n, m, A, mp_b, mp_N, mp_D, 258 * mp_NZ, mp_DZ) 259 * 260 * Summary: 261 * Certified solve a system of linear equations without reducing the 262 * solution size, where the left hand side input matrix is represented 263 * by signed long integers 264 * 265 * Description: 266 * Let the system of linear equations be Av = b, where A is a n x m matrix, 267 * and b is a n x 1 vector. There are three possibilities: 268 * 269 * 1. The system has more than one rational solution 270 * 2. The system has a unique rational solution 271 * 3. The system has no solution 272 * 273 * In the first case, there exist a solution vector v with minimal 274 * denominator and a rational certificate vector z to certify that the 275 * denominator of solution v is really the minimal denominator. 276 * 277 * The 1 x n certificate vector z satisfies that z.A is an integer vector 278 * and z.b has the same denominator as the solution vector v. 279 * In this case, the function will output the solution with minimal 280 * denominator and optional certificate vector z (if certflag = 1). 281 * 282 * Note: if choose not to compute the certificate vector z, the solution 283 * will not garantee, but with high probability, to be the minimal 284 * denominator solution, and the function will run faster. 285 * 286 * In the second case, the function will only compute the unique solution 287 * and the contents in the space for certificate vector make no sense. 288 * 289 * In the third case, there exists a certificate vector q to certify that 290 * the system has no solution. The 1 x n vector q satisfies q.A = 0 but 291 * q.b <> 0. In this case, the function will output this certificate vector 292 * q and store it into the same space for certificate z. The value of 293 * certflag also controls whether to output q or not. 294 * 295 * Note: if the function returns 3, then the system determinately does not 296 * exist solution, no matter whether to output certificate q or not. 297 * 298 * Input: 299 * certflag: 1/0, flag to indicate whether or not to compute the certificate 300 * vector z or q. 301 * - If certflag = 1, compute the certificate. 302 * - If certflag = 0, not compute the certificate. 303 * n: long, row dimension of the system 304 * m: long, column dimension of the system 305 * A: 1-dim signed long array length n*m, representation of n x m 306 * matrix A 307 * mp_b: 1-dim mpz_t array length n, representation of n x 1 vector b 308 * 309 * Return: 310 * 1: the first case, system has more than one solution 311 * 2: the second case, system has a unique solution 312 * 3: the third case, system has no solution 313 * 314 * Output: 315 * mp_N: 1-dim mpz_t array length m, 316 * - numerator vector of the solution with minimal denominator 317 * in the first case 318 * - numerator vector of the unique solution in the second case 319 * - make no sense in the third case 320 * mp_D: mpz_t, 321 * - minimal denominator of the solutions in the first case 322 * - denominator of the unique solution in the second case 323 * - make no sense in the third case 324 * 325 * The following will only be computed when certflag = 1 326 * mp_NZ: 1-dim mpz_t array length n, 327 * - numerator vector of the certificate z in the first case 328 * - make no sense in the second case 329 * - numerator vector of the certificate q in the third case 330 * mp_DZ: mpz_t, 331 * - denominator of the certificate z if in the first case 332 * - make no sense in the second case 333 * - denominator of the certificate q in the third case 334 * 335 * Note: 336 * - The space of (mp_N, mp_D) is needed to be preallocated, and entries in 337 * mp_N and integer mp_D are needed to be initiated as any integer values. 338 * - If certflag is specified to be 1, then also needs to preallocate space 339 * for (mp_NZ, mp_DZ), and initiate integer mp_DZ and entries in mp_NZ to 340 * be any integer values. 341 * Otherwise, set mp_NZ = NULL, and mp_DZ = any integer 342 * 343 */ 344 345 346 347extern long certSolveRedLong (const long certflag, const long nullcol, 348 const long n, const long m, const long *A, 349 mpz_t *mp_b, mpz_t *mp_N, mpz_t mp_D, 350 mpz_t *mp_NZ, mpz_t mp_DZ); 351 352/* 353 * 354 * Calling Sequence: 355 * 1/2/3 <-- certSolveRedLong(certflag, nullcol, n, m, A, mp_b, mp_N, mp_D, 356 * mp_NZ, mp_DZ) 357 * 358 * Summary: 359 * Certified solve a system of linear equations and reduce the solution 360 * size, where the left hand side input matrix is represented by signed 361 * long integers 362 * 363 * Description: 364 * Let the system of linear equations be Av = b, where A is a n x m matrix, 365 * and b is a n x 1 vector. There are three possibilities: 366 * 367 * 1. The system has more than one rational solution 368 * 2. The system has a unique rational solution 369 * 3. The system has no solution 370 * 371 * In the first case, there exist a solution vector v with minimal 372 * denominator and a rational certificate vector z to certify that the 373 * denominator of solution v is really the minimal denominator. 374 * 375 * The 1 x n certificate vector z satisfies that z.A is an integer vector 376 * and z.b has the same denominator as the solution vector v. 377 * In this case, the function will output the solution with minimal 378 * denominator and optional certificate vector z (if certflag = 1). 379 * 380 * Note: if choose not to compute the certificate vector z, the solution 381 * will not garantee, but with high probability, to be the minimal 382 * denominator solution, and the function will run faster. 383 * 384 * Lattice reduction will be used to reduce the solution size. Parameter 385 * nullcol designates the dimension of kernal basis we use to reduce the 386 * solution size as well as the dimension of nullspace we use to compute 387 * the minimal denominator. The heuristic results show that the solution 388 * size will be reduced by factor 1/nullcol. 389 * 390 * To find the minimum denominator as fast as possible, nullcol cannot be 391 * too small. We use NULLSPACE_COLUMN as the minimal value of nullcol. That 392 * is, if the input nullcol is less than NULLSPACE_COLUMN, NULLSPACE_COLUMN 393 * will be used instead. However, if the input nullcol becomes larger, the 394 * function will be slower. Meanwhile, it does not make sense to make 395 * nullcol greater than the dimension of nullspace of the input system. 396 * 397 * As a result, the parameter nullcol will not take effect unless 398 * NULLSPACE_COLUMN < nullcol < dimnullspace is satisfied, where 399 * dimnullspace is the dimension of nullspace of the input system. If the 400 * above condition is not satisfied, the boundary value NULLSPACE_COLUMN or 401 * dimnullspace will be used. 402 * 403 * In the second case, the function will only compute the unique solution 404 * and the contents in the space for certificate vector make no sense. 405 * 406 * In the third case, there exists a certificate vector q to certify that 407 * the system has no solution. The 1 x n vector q satisfies q.A = 0 but 408 * q.b <> 0. In this case, the function will output this certificate vector 409 * q and store it into the same space for certificate z. The value of 410 * certflag also controls whether to output q or not. 411 * 412 * Note: if the function returns 3, then the system determinately does not 413 * exist solution, no matter whether to output certificate q or not. 414 * 415 * Input: 416 * certflag: 1/0, flag to indicate whether or not to compute the certificate 417 * vector z or q. 418 * - If certflag = 1, compute the certificate. 419 * - If certflag = 0, not compute the certificate. 420 * nullcol: long, dimension of nullspace and kernel basis of conditioned 421 * system, 422 * if nullcol < NULLSPACE_COLUMN, use NULLSPACE_COLUMN instead 423 * n: long, row dimension of the system 424 * m: long, column dimension of the system 425 * A: 1-dim signed long array length n*m, representation of n x m 426 * matrix A 427 * mp_b: 1-dim mpz_t array length n, representation of n x 1 vector b 428 * 429 * Return: 430 * 1: the first case, system has more than one solution 431 * 2: the second case, system has a unique solution 432 * 3: the third case, system has no solution 433 * 434 * Output: 435 * mp_N: 1-dim mpz_t array length m, 436 * - numerator vector of the solution with minimal denominator 437 * in the first case 438 * - numerator vector of the unique solution in the second case 439 * - make no sense in the third case 440 * mp_D: mpz_t, 441 * - minimal denominator of the solutions in the first case 442 * - denominator of the unique solution in the second case 443 * - make no sense in the third case 444 * 445 * The following will only be computed when certflag = 1 446 * mp_NZ: 1-dim mpz_t array length n, 447 * - numerator vector of the certificate z in the first case 448 * - make no sense in the second case 449 * - numerator vector of the certificate q in the third case 450 * mp_DZ: mpz_t, 451 * - denominator of the certificate z if in the first case 452 * - make no sense in the second case 453 * - denominator of the certificate q in the third case 454 * 455 * Note: 456 * - The space of (mp_N, mp_D) is needed to be preallocated, and entries in 457 * mp_N and integer mp_D are needed to be initiated as any integer values. 458 * - If certflag is specified to be 1, then also needs to preallocate space 459 * for (mp_NZ, mp_DZ), and initiate integer mp_DZ and entries in mp_NZ to 460 * be any integer values. 461 * Otherwise, set mp_NZ = NULL, and mp_DZ = any integer 462 * 463 */ 464 465 466 467extern long certSolveMP (const long certflag, const long n, const long m, 468 mpz_t *mp_A, mpz_t *mp_b, mpz_t *mp_N, 469 mpz_t mp_D, mpz_t *mp_NZ, mpz_t mp_DZ); 470 471/* 472 * 473 * Calling Sequence: 474 * 1/2/3 <-- certSolveMP(certflag, n, m, mp_A, mp_b, mp_N, mp_D, 475 * mp_NZ, mp_DZ) 476 * 477 * Summary: 478 * Certified solve a system of linear equations without reducing the 479 * solution size, where the left hand side input matrix is represented 480 * by mpz_t integers 481 * 482 * Description: 483 * Let the system of linear equations be Av = b, where A is a n x m matrix, 484 * and b is a n x 1 vector. There are three possibilities: 485 * 486 * 1. The system has more than one rational solution 487 * 2. The system has a unique rational solution 488 * 3. The system has no solution 489 * 490 * In the first case, there exist a solution vector v with minimal 491 * denominator and a rational certificate vector z to certify that the 492 * denominator of solution v is really the minimal denominator. 493 * 494 * The 1 x n certificate vector z satisfies that z.A is an integer vector 495 * and z.b has the same denominator as the solution vector v. 496 * In this case, the function will output the solution with minimal 497 * denominator and optional certificate vector z (if certflag = 1). 498 * 499 * Note: if choose not to compute the certificate vector z, the solution 500 * will not garantee, but with high probability, to be the minimal 501 * denominator solution, and the function will run faster. 502 * 503 * In the second case, the function will only compute the unique solution 504 * and the contents in the space for certificate vector make no sense. 505 * 506 * In the third case, there exists a certificate vector q to certify that 507 * the system has no solution. The 1 x n vector q satisfies q.A = 0 but 508 * q.b <> 0. In this case, the function will output this certificate vector 509 * q and store it into the same space for certificate z. The value of 510 * certflag also controls whether to output q or not. 511 * 512 * Note: if the function returns 3, then the system determinately does not 513 * exist solution, no matter whether to output certificate q or not. 514 * In the first case, there exist a solution vector v with minimal 515 * denominator and a rational certificate vector z to certify that the 516 * denominator of solution v is the minimal denominator. 517 * 518 * Input: 519 * certflag: 1/0, flag to indicate whether or not to compute the certificate 520 * vector z or q. 521 * - If certflag = 1, compute the certificate. 522 * - If certflag = 0, not compute the certificate. 523 * n: long, row dimension of the system 524 * m: long, column dimension of the system 525 * mp_A: 1-dim mpz_t array length n*m, representation of n x m matrix A 526 * mp_b: 1-dim mpz_t array length n, representation of n x 1 vector b 527 * 528 * Return: 529 * 1: the first case, system has more than one solution 530 * 2: the second case, system has a unique solution 531 * 3: the third case, system has no solution 532 * 533 * Output: 534 * mp_N: 1-dim mpz_t array length m, 535 * - numerator vector of the solution with minimal denominator 536 * in the first case 537 * - numerator vector of the unique solution in the second case 538 * - make no sense in the third case 539 * mp_D: mpz_t, 540 * - minimal denominator of the solutions in the first case 541 * - denominator of the unique solution in the second case 542 * - make no sense in the third case 543 * 544 * The following will only be computed when certflag = 1 545 * mp_NZ: 1-dim mpz_t array length n, 546 * - numerator vector of the certificate z in the first case 547 * - make no sense in the second case 548 * - numerator vector of the certificate q in the third case 549 * mp_DZ: mpz_t, 550 * - denominator of the certificate z if in the first case 551 * - make no sense in the second case 552 * - denominator of the certificate q in the third case 553 * 554 * Note: 555 * - The space of (mp_N, mp_D) is needed to be preallocated, and entries in 556 * mp_N and integer mp_D are needed to be initiated as any integer values. 557 * - If certflag is specified to be 1, then also needs to preallocate space 558 * for (mp_NZ, mp_DZ), and initiate integer mp_DZ and entries in mp_NZ to 559 * be any integer values. 560 * Otherwise, set mp_NZ = NULL, and mp_DZ = any integer 561 * 562 */ 563 564 565 566extern long certSolveRedMP (const long certflag, const long nullcol, 567 const long n, const long m, mpz_t *mp_A, 568 mpz_t *mp_b, mpz_t *mp_N, mpz_t mp_D, 569 mpz_t *mp_NZ, mpz_t mp_DZ); 570 571/* 572 * 573 * Calling Sequence: 574 * 1/2/3 <-- certSolveRedMP(certflag, nullcol, n, m, mp_A, mp_b, mp_N, mp_D, 575 * mp_NZ, mp_DZ) 576 * 577 * Summary: 578 * Certified solve a system of linear equations and reduce the solution 579 * size, where the left hand side input matrix is represented by signed 580 * mpz_t integers 581 * 582 * Description: 583 * Let the system of linear equations be Av = b, where A is a n x m matrix, 584 * and b is a n x 1 vector. There are three possibilities: 585 * 586 * 1. The system has more than one rational solution 587 * 2. The system has a unique rational solution 588 * 3. The system has no solution 589 * 590 * In the first case, there exist a solution vector v with minimal 591 * denominator and a rational certificate vector z to certify that the 592 * denominator of solution v is really the minimal denominator. 593 * 594 * The 1 x n certificate vector z satisfies that z.A is an integer vector 595 * and z.b has the same denominator as the solution vector v. 596 * In this case, the function will output the solution with minimal 597 * denominator and optional certificate vector z (if certflag = 1). 598 * 599 * Note: if choose not to compute the certificate vector z, the solution 600 * will not garantee, but with high probability, to be the minimal 601 * denominator solution, and the function will run faster. 602 * 603 * Lattice reduction will be used to reduce the solution size. Parameter 604 * nullcol designates the dimension of kernal basis we use to reduce the 605 * solution size as well as the dimension of nullspace we use to compute 606 * the minimal denominator. The heuristic results show that the solution 607 * size will be reduced by factor 1/nullcol. 608 * 609 * To find the minimum denominator as fast as possible, nullcol cannot be 610 * too small. We use NULLSPACE_COLUMN as the minimal value of nullcol. That 611 * is, if the input nullcol is less than NULLSPACE_COLUMN, NULLSPACE_COLUMN 612 * will be used instead. However, if the input nullcol becomes larger, the 613 * function will be slower. Meanwhile, it does not make sense to make 614 * nullcol greater than the dimension of nullspace of the input system. 615 * 616 * As a result, the parameter nullcol will not take effect unless 617 * NULLSPACE_COLUMN < nullcol < dimnullspace is satisfied, where 618 * dimnullspace is the dimension of nullspace of the input system. If the 619 * above condition is not satisfied, the boundary value NULLSPACE_COLUMN or 620 * dimnullspace will be used. 621 * 622 * In the second case, the function will only compute the unique solution 623 * and the contents in the space for certificate vector make no sense. 624 * 625 * In the third case, there exists a certificate vector q to certify that 626 * the system has no solution. The 1 x n vector q satisfies q.A = 0 but 627 * q.b <> 0. In this case, the function will output this certificate vector 628 * q and store it into the same space for certificate z. The value of 629 * certflag also controls whether to output q or not. 630 * 631 * Note: if the function returns 3, then the system determinately does not 632 * exist solution, no matter whether to output certificate q or not. 633 * In the first case, there exist a solution vector v with minimal 634 * denominator and a rational certificate vector z to certify that the 635 * denominator of solution v is the minimal denominator. 636 * 637 * Input: 638 * certflag: 1/0, flag to indicate whether or not to compute the certificate 639 * vector z or q. 640 * - If certflag = 1, compute the certificate. 641 * - If certflag = 0, not compute the certificate. 642 * nullcol: long, dimension of nullspace and kernel basis of conditioned 643 * system, 644 * if nullcol < NULLSPACE_COLUMN, use NULLSPACE_COLUMN instead 645 * n: long, row dimension of the system 646 * m: long, column dimension of the system 647 * mp_A: 1-dim mpz_t array length n*m, representation of n x m matrix A 648 * mp_b: 1-dim mpz_t array length n, representation of n x 1 vector b 649 * 650 * Return: 651 * 1: the first case, system has more than one solution 652 * 2: the second case, system has a unique solution 653 * 3: the third case, system has no solution 654 * 655 * Output: 656 * mp_N: 1-dim mpz_t array length m, 657 * - numerator vector of the solution with minimal denominator 658 * in the first case 659 * - numerator vector of the unique solution in the second case 660 * - make no sense in the third case 661 * mp_D: mpz_t, 662 * - minimal denominator of the solutions in the first case 663 * - denominator of the unique solution in the second case 664 * - make no sense in the third case 665 * 666 * The following will only be computed when certflag = 1 667 * mp_NZ: 1-dim mpz_t array length n, 668 * - numerator vector of the certificate z in the first case 669 * - make no sense in the second case 670 * - numerator vector of the certificate q in the third case 671 * mp_DZ: mpz_t, 672 * - denominator of the certificate z if in the first case 673 * - make no sense in the second case 674 * - denominator of the certificate q in the third case 675 * 676 * Note: 677 * - The space of (mp_N, mp_D) is needed to be preallocated, and entries in 678 * mp_N and integer mp_D are needed to be initiated as any integer values. 679 * - If certflag is specified to be 1, then also needs to preallocate space 680 * for (mp_NZ, mp_DZ), and initiate integer mp_DZ and entries in mp_NZ to 681 * be any integer values. 682 * Otherwise, set mp_NZ = NULL, and mp_DZ = any integer 683 * 684 */ 685 686 687 688extern void RowEchelonTransform (const FiniteField p, Double *A, const long n, 689 const long m, const long frows, 690 const long lrows, const long redflag, 691 const long eterm, long *Q, long *rp, 692 FiniteField *d); 693 694/* 695 * Calling Sequence: 696 * RowEchelonTransform(p, A, n, m, frows, lrows, redflag, eterm, Q, rp, d) 697 * 698 * Summary: 699 * Compute a mod p row-echelon transform of a mod p input matrix 700 * 701 * Description: 702 * Given a n x m mod p matrix A, a row-echelon transform of A is a 4-tuple 703 * (U,P,rp,d) with rp the rank profile of A (the unique and strictly 704 * increasing list [j1,j2,...jr] of column indices of the row-echelon form 705 * which contain the pivots), P a permutation matrix such that all r leading 706 * submatrices of (PA)[0..r-1,rp] are nonsingular, U a nonsingular matrix 707 * such that UPA is in row-echelon form, and d the determinant of 708 * (PA)[0..r-1,rp]. 709 * 710 * Generally, it is required that p be a prime, as inverses are needed, but 711 * in some cases it is possible to obtain an echelon transform when p is 712 * composite. For the cases where the echelon transform cannot be obtained 713 * for p composite, the function returns an error indicating that p is 714 * composite. 715 * 716 * The matrix U is structured, and has last n-r columns equal to the last n-r 717 * columns of the identity matrix, n the row dimension of A. 718 * 719 * The first r rows of UPA comprise a basis in echelon form for the row 720 * space of A, while the last n-r rows of U comprise a basis for the left 721 * nullspace of PA. 722 * 723 * For efficiency, this function does not output an echelon transform 724 * (U,P,rp,d) directly, but rather the expression sequence (Q,rp,d). 725 * Q, rp, d are the form of arrays and pointers in order to operate inplace, 726 * which require to preallocate spaces and initialize them. Initially, 727 * Q[i] = i (i=0..n), rp[i] = 0 (i=0..n), and *d = 1. Upon completion, rp[0] 728 * stores the rank r, rp[1..r] stores the rank profile. i<=Q[i]<=n for 729 * i=1..r. The input Matrix A is modified inplace and used to store U. 730 * Let A' denote the state of A on completion. Then U is obtained from the 731 * identity matrix by replacing the first r columns with those of A', and P 732 * is obtained from the identity matrix by swapping row i with row Q[i], for 733 * i=1..r in succession. 734 * 735 * Parameters flrows, lrows, redflag, eterm control the specific operations 736 * this function will perform. Let (U,P,rp,d) be as constructed above. If 737 * frows=0, the first r rows of U will not be correct. If lrows=0, the last 738 * n-r rows of U will not be correct. The computation can be up to four 739 * times faster if these flags are set to 0. 740 * 741 * If redflag=1, the row-echelon form is reduced, that is (UPA)[0..r-1,rp] 742 * will be the identity matrix. If redflag=0, the row-echelon form will not 743 * be reduced, that is (UPA)[1..r,rp] will be upper triangular and U is unit 744 * lower triangular. If frows=0 then redflag has no effect. 745 * 746 * If eterm=1, then early termination is triggered if a column of the 747 * input matrix is discovered that is linearly dependant on the previous 748 * columns. In case of early termination, the third return value d will be 0 749 * and the remaining components of the echelon transform will not be correct. 750 * 751 * Input: 752 * p: FiniteField, modulus 753 * A: 1-dim Double array length n*m, representation of a n x m input 754 * matrix 755 * n: long, row dimension of A 756 * m: long, column dimension of A 757 * frows: 1/0, 758 * - if frows = 1, the first r rows of U will be correct 759 * - if frows = 0, the first r rows of U will not be correct 760 * lrows: 1/0, 761 * - if lrows = 1, the last n-r rows of U will be correct 762 * - if lrows = 0, the last n-r rows of U will not be correct 763 * redflag: 1/0, 764 * - if redflag = 1, compute row-echelon form 765 * - if redflag = 0, not compute reow-echelon form 766 * eterm: 1/0, 767 * - if eterm = 1, terminate early if not in full rank 768 * - if eterm = 0, not terminate early 769 * Q: 1-dim long array length n+1, compact representation of 770 * permutation vector, initially Q[i] = i, 0 <= i <= n 771 * rp: 1-dim long array length n+1, representation of rank profile, 772 * initially rp[i] = 0, 0 <= i <= n 773 * d: pointer to FiniteField, storing determinant of the matrix, 774 * initially *d = 1 775 * 776 * Precondition: 777 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2) 778 * 779 */ 780 781 782 783 784extern Double * mAdjoint (const FiniteField p, Double *A, const long n); 785 786/* 787 * Calling Sequence: 788 * Adj <-- mAdjoint(p, A, n) 789 * 790 * Summary: 791 * Compute the adjoint of a mod p square matrix 792 * 793 * Description: 794 * Given a n x n mod p matrix A, the function computes adjoint of A. Input 795 * A is not modified upon completion. 796 * 797 * Input: 798 * p: FiniteField, prime modulus 799 * if p is a composite number, the routine will still work if no error 800 * message is returned 801 * A: 1-dim Double array length n*n, representation of a n x n mod p matrix. 802 * The entries of A are casted from integers 803 * n: long, dimension of A 804 * 805 * Return: 806 * 1-dim Double matrix length n*n, repesentation of a n x n mod p matrix, 807 * adjoint of A 808 * 809 * Precondition: 810 * n*(p-1)^2 <= 2^53-1 = 9007199254740991 811 * 812 */ 813 814 815extern long mBasis (const FiniteField p, Double *A, const long n, 816 const long m, const long basis, const long nullsp, 817 Double **B, Double **N); 818 819/* 820 * Calling Sequence: 821 * r/-1 <-- mBasis(p, A, n, m, basis, nullsp, B, N) 822 * 823 * Summary: 824 * Compute a basis for the rowspace and/or a basis for the left nullspace 825 * of a mod p matrix 826 * 827 * Description: 828 * Given a n x m mod p matrix A, the function computes a basis for the 829 * rowspace B and/or a basis for the left nullspace N of A. Row vectors in 830 * the r x m matrix B consist of basis of A, where r is the rank of A in 831 * Z/pZ. If r is zero, then B will be NULL. Row vectors in the n-r x n 832 * matrix N consist of the left nullspace of A. N will be NULL if A is full 833 * rank. 834 * 835 * The pointers are passed into argument lists to store the computed basis 836 * and nullspace. Upon completion, the rank r will be returned. The 837 * parameters basis and nullsp control whether to compute basis and/or 838 * nullspace. If set basis and nullsp in the way that both basis and 839 * nullspace will not be computed, an error message will be printed and 840 * instead of rank r, -1 will be returned. 841 * 842 * Input: 843 * p: FiniteField, prime modulus 844 * if p is a composite number, the routine will still work if no 845 * error message is returned 846 * A: 1-dim Double array length n*m, representation of a n x m mod p 847 * matrix. The entries of A are casted from integers 848 * n: long, row dimension of A 849 * m: long, column dimension of A 850 * basis: 1/0, flag to indicate whether to compute basis for rowspace or 851 * not 852 * - basis = 1, compute the basis 853 * - basis = 0, not compute the basis 854 * nullsp: 1/0, flag to indicate whether to compute basis for left nullspace 855 * or not 856 * - nullsp = 1, compute the nullspace 857 * - nullsp = 0, not compute the nullspace 858 * 859 * Output: 860 * B: pointer to (Double *), if basis = 1, *B will be a 1-dim r*m Double 861 * array, representing the r x m basis matrix. If basis = 1 and r = 0, 862 * *B = NULL 863 * 864 * N: pointer to (Double *), if nullsp = 1, *N will be a 1-dim (n-r)*n Double 865 * array, representing the n-r x n nullspace matrix. If nullsp = 1 and 866 * r = n, *N = NULL. 867 * 868 * Return: 869 * - if basis and/or nullsp are set to be 1, then return the rank r of A 870 * - if both basis and nullsp are set to be 0, then return -1 871 * 872 * Precondition: 873 * n*(p-1)^2 <= 2^53-1 = 9007199254740991 874 * 875 * Note: 876 * - In case basis = 0, nullsp = 1, A will be destroyed inplace. Otherwise, 877 * A will not be changed. 878 * - Space of B and/or N will be allocated in the function 879 * 880 */ 881 882 883 884extern long mDeterminant (const FiniteField p, Double *A, const long n); 885 886/* 887 * Calling Sequence: 888 * det <-- mDeterminant(p, A, n) 889 * 890 * Summary: 891 * Compute the determinant of a square mod p matrix 892 * 893 * Input: 894 * p: FiniteField, prime modulus 895 * if p is a composite number, the routine will still work if no error 896 * message is returned 897 * A: 1-dim Double array length n*n, representation of a n x n mod p matrix. 898 * The entries of A are casted from integers 899 * n: long, dimension of A 900 * 901 * Output: 902 * det(A) mod p, the determinant of square matrix A 903 * 904 * Precondition: 905 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2) 906 * 907 * Note: 908 * A is destroyed inplace 909 * 910 */ 911 912 913extern long mInverse (const FiniteField p, Double *A, const long n); 914 915/* 916 * Calling Sequence: 917 * 1/0 <-- mInverse(p, A, n) 918 * 919 * Summary: 920 * Certified compute the inverse of a mod p matrix inplace 921 * 922 * Description: 923 * Given a n x n mod p matrix A, the function computes A^(-1) mod p 924 * inplace in case A is a nonsingular matrix in Z/Zp. If the inverse does 925 * not exist, the function returns 0. 926 * 927 * A will be destroyed at the end in both cases. If the inverse exists, A is 928 * inplaced by its inverse. Otherwise, the inplaced A is not the inverse. 929 * 930 * Input: 931 * p: FiniteField, prime modulus 932 * if p is a composite number, the routine will still work if no error 933 * message is returned 934 * A: 1-dim Double array length n*n, representation of a n x n mod p matrix. 935 * The entries of A are casted from integers 936 * n: long, dimension of A 937 * 938 * Return: 939 * - 1, if A^(-1) mod p exists 940 * - 0, if A^(-1) mod p does not exist 941 * 942 * Precondition: 943 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2) 944 * 945 * Note: 946 * A is destroyed inplace 947 * 948 */ 949 950 951 952extern long mRank (const FiniteField p, Double *A, const long n, const long m); 953 954/* 955 * Calling Sequence: 956 * r <-- mRank(p, A, n, m) 957 * 958 * Summary: 959 * Compute the rank of a mod p matrix 960 * 961 * Input: 962 * p: FiniteField, prime modulus 963 * if p is a composite number, the routine will still work if no 964 * error message is returned 965 * A: 1-dim Double array length n*m, representation of a n x m mod p 966 * matrix. The entries of A are casted from integers 967 * n: long, row dimension of A 968 * m: long, column dimension of A 969 * 970 * Return: 971 * r: long, rank of matrix A 972 * 973 * Precondition: 974 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2) 975 * 976 * Note: 977 * A is destroyed inplace 978 * 979 */ 980 981 982 983extern long * mRankProfile (const FiniteField p, Double *A, 984 const long n, const long m); 985 986/* 987 * Calling Sequence: 988 * rp <-- mRankProfile(p, A, n, m) 989 * 990 * Summary: 991 * Compute the rank profile of a mod p matrix 992 * 993 * Input: 994 * p: FiniteField, prime modulus 995 * if p is a composite number, the routine will still work if no 996 * error message is returned 997 * A: 1-dim Double array length n*m, representation of a n x m mod p 998 * matrix. The entries of A are casted from integers 999 * n: long, row dimension of A 1000 * m: long, column dimension of A 1001 * 1002 * Return: 1003 * rp: 1-dim long array length n+1, where 1004 * - rp[0] is the rank of matrix A 1005 * - rp[1..r] is the rank profile of matrix A 1006 * 1007 * Precondition: 1008 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2) 1009 * 1010 * Note: 1011 * A is destroyed inplace 1012 * 1013 */ 1014