uppgb.c (3198B)
1 #include "common.h" 2 3 static void systemsolver(fmpq_t res, pt_t *p, pt_t *d, size_t n) { 4 fmpq_mat_t A, x, b; 5 fmpq_mat_init(A, 6, 6); // should be enough unless linear dep 6 fmpq_mat_init(x, 6, 1); 7 fmpq_mat_init(b, 6, 1); 8 9 fmpq_t t; 10 fmpq_init(t); 11 12 size_t i, j; 13 size_t ii = 0; 14 for (i = 0; i < n; i++) { 15 for (j = i + 1; j < n; j++) { 16 // A[ii,0] = d[j,1] - d[i,1] 17 fmpq_sub(fmpq_mat_entry(A, ii, 0), d[j].c[1], d[i].c[1]); 18 // A[ii,1] = d[i,0] - d[j,0] 19 fmpq_sub(fmpq_mat_entry(A, ii, 1), d[i].c[0], d[j].c[0]); 20 // A[ii,2] = 0 21 fmpq_set_ui(fmpq_mat_entry(A, ii, 2), 0, 1); 22 // A[ii,3] = p[i,1] - p[j,1] 23 fmpq_sub(fmpq_mat_entry(A, ii, 3), p[i].c[1], p[j].c[1]); 24 // A[ii,4] = p[j,0] - p[i,0] 25 fmpq_sub(fmpq_mat_entry(A, ii, 4), p[j].c[0], p[i].c[0]); 26 // A[ii,5] = 0 27 fmpq_set_ui(fmpq_mat_entry(A, ii, 5), 0, 1); 28 // b[ii] = -d[i,1]p[i,0]+d[i,0]p[i,1]+d[j,1]p[j,0]-d[j,0]p[j,1] 29 fmpq_set_ui(fmpq_mat_entry(b, ii, 0), 0, 1); 30 fmpq_mul(t, d[i].c[1], p[i].c[0]); 31 fmpq_sub(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 32 fmpq_mul(t, d[i].c[0], p[i].c[1]); 33 fmpq_add(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 34 fmpq_mul(t, d[j].c[1], p[j].c[0]); 35 fmpq_add(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 36 fmpq_mul(t, d[j].c[0], p[j].c[1]); 37 fmpq_sub(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 38 39 ii++; 40 if (ii == 6) 41 break; 42 43 // A[ii,0] = d[j,2] - d[i,2] 44 fmpq_sub(fmpq_mat_entry(A, ii, 0), d[j].c[2], d[i].c[2]); 45 // A[ii,2] = d[i,0] - d[j,0] 46 fmpq_sub(fmpq_mat_entry(A, ii, 2), d[i].c[0], d[j].c[0]); 47 // A[ii,1] = 0 48 fmpq_set_ui(fmpq_mat_entry(A, ii, 1), 0, 1); 49 // A[ii,3] = p[i,2] - p[j,2] 50 fmpq_sub(fmpq_mat_entry(A, ii, 3), p[i].c[2], p[j].c[2]); 51 // A[ii,5] = p[j,0] - p[i,0] 52 fmpq_sub(fmpq_mat_entry(A, ii, 5), p[j].c[0], p[i].c[0]); 53 // A[ii,4] = 0 54 fmpq_set_ui(fmpq_mat_entry(A, ii, 4), 0, 1); 55 // b[ii] = -d[i,2]p[i,0]+d[i,0]p[i,2]+d[j,2]p[j,0]-d[j,0]p[j,2] 56 fmpq_set_ui(fmpq_mat_entry(b, ii, 0), 0, 1); 57 fmpq_mul(t, d[i].c[2], p[i].c[0]); 58 fmpq_sub(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 59 fmpq_mul(t, d[i].c[0], p[i].c[2]); 60 fmpq_add(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 61 fmpq_mul(t, d[j].c[2], p[j].c[0]); 62 fmpq_add(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 63 fmpq_mul(t, d[j].c[0], p[j].c[2]); 64 fmpq_sub(fmpq_mat_entry(b, ii, 0), fmpq_mat_entry(b, ii, 0), t); 65 66 ii++; 67 if (ii == 6) 68 break; 69 } 70 71 if (ii == 6) 72 break; 73 } 74 75 fmpq_mat_solve(x, A, b); 76 77 fmpq_set(res, fmpq_mat_entry(x, 0, 0)); 78 fmpq_add(res, res, fmpq_mat_entry(x, 1, 0)); 79 fmpq_add(res, res, fmpq_mat_entry(x, 2, 0)); 80 81 fmpq_mat_clear(x); 82 fmpq_clear(t); 83 fmpq_mat_clear(A); 84 } 85 86 int main(int argc, char **argv) { 87 char **lines; 88 size_t nlines = readlines(&lines, "input"); 89 90 pt_t p[nlines], d[nlines]; 91 read_pts(p, d, lines, nlines); 92 93 fmpq_t res; 94 fmpq_init(res); 95 96 systemsolver(res, p, d, nlines); 97 98 fmpq_print(res); 99 printf("\n"); 100 }