/* This is a C library for multiplication, division, inversion and GCD computation of univariate polynomials over an algebraic number field. Seyed Mohammad Mahdi Javadi 2008~2009 */ #include #include #define ULONG unsigned long long int #define LONG long long int #define poly LONG * #define array int * #define longarray LONG * /* min and max functions */ #define min(X,Y) ((X) < (Y) ? (X) : (Y)) #define max(X,Y) ((X) > (Y) ? (X) : (Y)) void ip_rem_monic_b(poly A, poly B, int N, array Si, poly E, LONG M, longarray WS); void ip_rem_no_ext(poly A, poly B, LONG M); /* compute the amount of space needed to store a polynomial of degree d in x over a number field with N extensions. Si[i] is S_i in the paper. S_i is the amount of space needed to store one element in R_i. Si[1] = 1. */ int space_d_N(int d, int N, array Si) { return (d+1)*Si[N] + 1; } /* Copy the polynomial B into the space A. */ void copy_into(poly A, poly B, int N, array Si) { size_t LEN = space_d_N(B[0], N, Si)*sizeof(LONG); memcpy(A, B, LEN); } /* A := A + B. M here is the modulus. */ void ip_add(poly A, poly B, int N, array Si, LONG M) { int da, db, i, block_size, sb, dr; LONG t, s; da = A[0]; db = B[0]; if (db == -1) return; if (da == -1) { copy_into(A,B,N,Si); return; } if (N == 0) { if (da > db) { dr = da; for (i = 1; i <= db+1; i++) { t = A[i] + B[i]; s = t - M; if (s < 0) A[i] = t; else A[i] = s; } } else { dr = db; for (i = 1; i <= da+1; i++) { t = A[i] + B[i]; s = t - M; if (s < 0) A[i] = t; else A[i] = s; } for (; i <= db + 1; i++) A[i] = B[i]; } for (; (dr >= 0) && (A[dr+1] == 0); dr--); A[0] = dr; return; } block_size = Si[N]; if (da > db) { sb = 1; dr = da; for (i = 0; i <= db; i++, sb += block_size) ip_add(A+sb, B+sb, N-1,Si, M); } else { dr = db; sb = 1; for (i = 0; i <= da; i++, sb += block_size) ip_add(A+sb, B+sb, N-1,Si, M); if (da < db) { for (; i <= db; i++, sb += block_size) copy_into(A+sb,B+sb,N-1, Si); } } for (sb = dr*block_size + 1; (dr >= 0) && (A[sb] == -1); dr--, sb -= block_size); A[0] = dr; return; } /* A := -A; * Destroys A */ void ip_neg(poly A, int N, array Si, LONG M) { int i, block_size, s; if (N == 0) { for (i = 1; i <= A[0] + 1; i++) if (A[i] > 0) A[i] = M-A[i]; return; } block_size = Si[N]; for (i = 1, s = 1; i <= A[0] + 1; i++, s += block_size) ip_neg(A+s, N-1, Si, M); return; } /* A := A - B; */ void ip_sub(poly A, poly B, int N, array Si, LONG M) { int da, db, i, block_size, sb, dr; da = A[0]; db = B[0]; if (db == -1) return; if (da == -1) { ip_neg(B, N, Si, M); copy_into(A,B,N,Si); return; } if (N == 0) { if (da > db) { dr = da; for (i = 1; i <= db+1; i++) { A[i] = A[i] - B[i]; if (A[i] < 0LL) A[i] = M + A[i]; } } else { dr = db; for (i = 1; i <= da+1; i++) { A[i] = A[i] - B[i]; if (A[i] < 0LL) A[i] = M + A[i]; } for (; i <= db + 1; i++) { if (B[i] == 0) A[i] = 0; else A[i] = M - B[i]; } } for (; (dr >= 0) && (A[dr+1] == 0); dr--); A[0] = dr; return; } block_size = Si[N]; if (da > db) { dr = da; sb = 1; for (i = 0; i <= db; i++, sb += block_size) ip_sub(A+sb, B+sb, N-1,Si, M); } else { dr = db; sb = 1; for (i = 0; i <= da; i++, sb += block_size) ip_sub(A+sb, B+sb, N-1,Si, M); if (da < db) { for (; i <= db; i++, sb += block_size) { ip_neg(B+sb, N-1, Si, M); copy_into(A+sb,B+sb,N-1, Si); } } } for (sb = dr*block_size + 1; (dr >= 0) && (A[sb] == -1); dr--, sb -= block_size); A[0] = dr; return; } /* C := A * B degree(A) < 2 and degree(B) < 2 No field extension. */ void ip_mul_no_ext_quadratic(poly A, poly B, poly C, LONG M) { ULONG t,s; C[1] = (A[1]*B[1]) % M; if (A[0] == 0) { if (B[0] == 0) { if (C[1] == 0) C[0] = -1; else C[0] = 0; return; } C[2] = (A[1]*B[2]) % M; if (C[2] > 0) C[0] = 1; else if (C[1] > 0) C[0] = 0; else C[0] = -1; return; } if (B[0] == 0) { C[2] = (A[2]*B[1]) % M; if (C[2] > 0) C[0] = 1; else if (C[1] > 0) C[0] = 0; else C[0] = -1; return; } t = A[2]*B[1]; s = A[1]*B[2]; t = (s + t) % (ULONG)M; C[2] = t; C[3] = (A[2]*B[2]) % M; if (C[3] > 0) C[0] = 2; else if (C[2] > 0) C[0] = 1; else if(C[1] > 0) C[0] = 0; else C[0] = -1; return; } /* C := A * B. No field extension. */ void ip_mul_no_ext(poly A, poly B, poly C, LONG M) { int dc, k, da, db, i, mm; LONG t1, X; da = A[0]; db = B[0]; if ((da == -1) || (db == -1)) { C[0] = -1; return; } if ((da < 2) && (db < 2)) return ip_mul_no_ext_quadratic(A, B, C, M); X = M*M; dc = A[0] + B[0]; for (k = 0; k <= dc; k++) { t1 = 0; mm = min(k, da); for (i = max(0,k - db); i <= mm; i++) { if (t1 > 0) t1 -= X; t1 += A[i + 1]*B[k - i + 1]; } t1 %= M; if (t1 < 0) t1 += M; C[k+1] = t1; } for(; (dc >= 0) && (C[dc+1] == 0); dc--); C[0] = dc; return; } /* Reduce the coefficients of A modulu M. */ void reduce_int_coeff(poly A, LONG M) { int da, i; da = A[0]; for (i = 1; i <= da + 1; i++) { A[i] %= M; if (A[i] < 0) A[i] += M; } for( i=da+1; i>0 && A[i]==0; i-- ); A[0] = i-1; return; } /* C := C + A*B. No field extensions. MS = M^2. */ void ip_muladd_no_ext(poly A, poly B, poly C, LONG M, LONG MS) { int da, db, i, j, k; LONG t; da = A[0]; db = B[0]; if ((da == -1) || (db == -1)) return; for (i = 0; i <= da; i++) { t = A[i+1]; for (j=1; j<=db+1; j++) { k = i+j; if( C[k] > 0 ) C[k] -= MS; C[k] += t*B[j]; } } return; } /* C(z) = C(z) - A(z)*B(z) */ void ip_mulsub_no_ext(poly A, poly B, poly C, LONG M, LONG MS) { int da, db, i, j, k; LONG t; da = A[0]; db = B[0]; if ((da == -1) || (db == -1)) return; for (i = 0; i <= da; i++) { t = A[i+1]; for (j=1; j<=db+1; j++) { k = i+j; if( C[k] < 0LL ) C[k] += MS; C[k] -= t*B[j]; } } return; } /* C := A * B Special multiplication routine for R_1[x] MS = M^2 */ void ip_mul_one_ext(poly A, poly B, poly C, array Si, poly E, LONG M, longarray WS, LONG MS) { int da, db, dm, dc, block_size, k, i, sa, sb, sc, j; poly minpoly; poly t1; poly t2; da = A[0]; db = B[0]; if ( (da == -1) || (db == -1)) { C[0] = -1; return; } /* A = 0 or B = 0 ==> C = 0 */ minpoly = E; dm = minpoly[0]; t1 = WS; WS += 2*dm; t2 = WS; WS += 2*dm; dc = da + db; block_size = Si[1]; sc = 1; for (k = 0; k <= dc; k++) { i = max(0,k-db); sa = 1+block_size*i; sb = 1+block_size*(k-i); t1[0] = 2*dm-2; for( j=1; j<=t1[0]+1; j++ ) t1[j] = 0; for(; i <= min(k, da); i++, sa += block_size, sb -= block_size) { ip_muladd_no_ext(A+sa,B+sb,t1,M,MS); } reduce_int_coeff(t1, M); ip_rem_no_ext(t1, minpoly, M); copy_into(C+sc, t1, 0, Si); sc += block_size; } for(i = dc, sc -= block_size; (i >= 0) && (C[sc] == -1) ; i--, sc -= block_size); C[0] = i; return; } /* Main multilication routine. C := A * B WS: Working storage */ void ip_mul(poly A, poly B, poly C, int N, array Si, poly E, LONG M, longarray WS) { int da, db, dm, dc, block_size, k, i, sa, sb, sc; poly minpoly; poly F; poly t1; poly t2; LONG MS; da = A[0]; db = B[0]; if ( (da == -1) || (db == -1)) { C[0] = -1; return; } /* A = 0 or B = 0 ==> C = 0 */ if (N == 0) { ip_mul_no_ext(A, B, C, M); return; } if (N == 1) { MS = M*M; return ip_mul_one_ext(A, B, C, Si, E, M, WS, MS); } minpoly = E; dm = minpoly[0]; t1 = WS; WS += space_d_N(2*dm-2, N-1, Si); t2 = WS; WS += space_d_N(2*dm-2, N-1, Si); dc = da + db; block_size = Si[N]; F = E + block_size + Si[N-1]; sc = 1; for (k = 0; k <= dc; k++) { i = max(0,k-db); sa = 1+block_size*i; sb = 1+block_size*(k-i); ip_mul(A+sa,B+sb,t2,N-1,Si,F,M,WS); copy_into(t1, t2, N - 1, Si); for(i++, sa += block_size, sb -= block_size; i <= min(k, da); i++, sa += block_size, sb -= block_size) { ip_mul(A+sa,B+sb,t2,N-1,Si,F,M,WS); ip_add(t1, t2, N - 1, Si, M); } ip_rem_monic_b(t1, minpoly, N - 1, Si, F, M, WS); copy_into(C+sc, t1, N - 1, Si); sc += block_size; } for(i = dc, sc -= block_size; (i >= 0) && (C[sc] == -1) ; i--, sc -= block_size); C[0] = i; return; } /* C := A(z) * B(z) mod m(z) where E contains m(z) and degree(A) < 2 and degree(B) < 2. */ void ip_mul_mod_one_ext_quadratic(poly A, poly B, poly C, poly E, LONG M) { LONG LC; if ((A[0] == -1) || (B[0] == -1)) { C[0] = -1; return; } if ((A[0] == 0) || (B[0] == 0)) return ip_mul_no_ext_quadratic(A, B, C, M); LC = (A[2]*B[2]) % M; C[1] = ((A[1]*B[1]) % M - (LC*E[1]) ) % M; C[1] = ((A[1]*B[1]) % M - (LC*E[1] % M) ) % M; if (C[1] < 0LL) C[1] += M; C[2] = ((A[2]*B[1]) % M + (A[1]*B[2]) % M - (LC*E[2])) % M; C[2] = ((A[2]*B[1]) % M + (A[1]*B[2]) % M - (LC*E[2] % M)) % M; if (C[2] < 0LL) C[2] += M; if (C[2] > 0) C[0] = 1; else if (C[1] > 0) C[0] = 0; else C[0] = -1; return; } /* C := A(z) * B(z) mod m(z) where E contains m(z) */ void ip_mul_mod(poly A, poly B, poly C, int N, array Si, poly E, LONG M, longarray WS) { int bs; if (N == 0) bs = Si[N]; else bs = Si[N] + Si[N-1]; if ((N == 1) && (E[0] == 2)) return ip_mul_mod_one_ext_quadratic(A, B, C,E, M); ip_mul(A, B, C, N - 1, Si, E+bs, M, WS); ip_rem_monic_b(C, E, N-1, Si, E+bs,M, WS); } /* Store the remainder and quotient of dividing A by B in A. No field extension. */ void ip_rem_no_ext(poly A, poly B, LONG M) { int da, db, dq, dr, k, j, sa, sb, mm; LONG t1, X; da = A[0]; db = B[0]; if (da < db) return; dq = da - db; dr = db - 1; X = M*M; for (k = da; k >= 0; k--) { t1 = A[k+1]; j = max(0, k - dq); sb = j + 1; sa = k-j+db + 1; mm = min(dr, k); for(; j <= mm; j++, sb++,sa--) { if (t1 < 0) t1 += X; t1 -= B[sb]*A[sa]; } t1 = t1 % M; if (t1 < 0) t1 += M; A[k+1] = t1; } for(; (dr >= 0) && (A[dr+1] == 0); dr--); A[0] = dr; return; } /* Special subtraction routine for the case where there is no field extension. Does not reduce the coefficients mod M. MS = M^2. */ void ip_sub_no_ext_special(poly A, poly B, array Si, LONG M, LONG MS) { int da, db, i, dr; da = A[0]; db = B[0]; if (db == -1) return; if (da == -1) { A[0] = B[0]; for (i = 1; i <= B[0] + 1; i++) A[i] = (-B[i]) - MS; return; } if (da > db) { dr = da; for (i = 1; i <= db+1; i++) { A[i] = A[i] - B[i]; if (A[i] > 0) A[i] -= MS; } } else { dr = db; for (i = 1; i <= da+1; i++) { A[i] = A[i] - B[i]; if (A[i] > 0) A[i] -= MS; } for (; i <= db + 1; i++) { if (B[i] == 0) A[i] = 0; else A[i] = (-B[i])-MS; } } for (; (dr >= 0) && (A[dr+1] == 0); dr--); A[0] = dr; return; } /* Special multiplication routine. Does not reduce the coefficients modulo M. There are no extensions. */ void ip_mul_no_ext_special(poly A, poly B, poly C, LONG M, LONG MS) { int dc, k, da, db, i, mm; LONG t1; da = A[0]; db = B[0]; if ((da == -1) || (db == -1)) { C[0] = -1; return; } dc = A[0] + B[0]; for (k = 0; k <= dc; k++) { t1 = 0; mm = min(k, da); for (i = max(0,k - db); i <= mm; i++) { t1 += A[i + 1]*B[k - i + 1]; if (t1 > 0) t1 -= MS; } C[k+1] = t1; } for(; (dc >= 0) && (C[dc+1] == 0); dc--); C[0] = dc; return; } /* Store the remainder and quotient of dividing A by B in A. This is a special routine for the case where there is only one extension. */ void ip_rem_one_ext(poly A, poly B, array Si, poly E, LONG M, longarray WS, LONG MS) { int da, db, block_size, dm, i, dr, dq, k, sa, sb, sc, j; poly minpoly; poly t1; poly t2; da = A[0]; db = B[0]; if (da < db) return; dq = da - db; dr = db - 1; minpoly = E; dm = minpoly[0]; t1 = WS; WS += 2*dm; t2 = WS; WS += 2*dm; block_size = Si[1]; sc = 1 + da*block_size; for (k = da; k >= 0; k--) { t1[0] = 2*dm-2; for( j=1; j<=t1[0]+1; j++ ) t1[j] = 0; copy_into(t1, A+sc, 0, Si); t1[0] = 2*dm-2; i = max(0, k-dq); sb = 1+i*block_size; sa = 1+(k-i+db)*block_size; for (; i <= min(dr, k); i++, sb += block_size, sa -= block_size) { ip_mulsub_no_ext(A+sa,B+sb,t1,M,MS); } reduce_int_coeff(t1, M); ip_rem_no_ext(t1, minpoly, M); copy_into(A+sc,t1, 0, Si); sc -= block_size; } for(i = dr, sc = 1+dr*block_size; (i >= 0) && (A[sc] == -1) ; i--, sc -= block_size); A[0] = i; return; } /* Main division routine: stroe the remainder and the quotient of dividing A by B in A. Note: There is not enough space to store the degree of the quotient. */ void ip_rem_monic_b(poly A, poly B, int N, array Si, poly E, LONG M, longarray WS) { int da, db, block_size, dm, i, dr, dq, k, sa, sb, sc; poly minpoly; poly t1; poly t2; poly F; LONG MS; da = A[0]; db = B[0]; if (da < db) return; if (N == 0) { ip_rem_no_ext(A, B, M); return; } if (N == 1) { MS = M*M; return ip_rem_one_ext(A, B, Si, E, M, WS, MS); } dq = da - db; dr = db - 1; minpoly = E; dm = minpoly[0]; t1 = WS; WS += space_d_N(2*dm-2, N-1, Si); t2 = WS; WS += space_d_N(2*dm-2, N-1, Si); block_size = Si[N]; F = E + block_size + Si[N - 1]; sc = 1 + da*block_size; for (k = da; k >= 0; k--) { copy_into(t1, A+sc, N - 1, Si); i = max(0, k-dq); sb = 1+i*block_size; sa = 1+(k-i+db)*block_size; for (; i <= min(dr, k); i++, sb += block_size, sa -= block_size) { ip_mul(A+sa,B+sb,t2,N-1,Si,F,M,WS); ip_sub(t1,t2,N-1,Si,M); } ip_rem_monic_b(t1, minpoly, N - 1, Si, F, M, WS); copy_into(A+sc,t1, N - 1, Si); sc -= block_size; } for(i = dr, sc = 1+dr*block_size; (i >= 0) && (A[sc] == -1) ; i--, sc -= block_size); A[0] = i; return; } /* A := A * c. Multiply A by an scalar c. There is no field extension. */ void ip_scale_mul_no_ext(poly A, LONG c, LONG M) { int da, i; da = A[0]; if (da == -1) return; if (c == 0) { A[0] = -1; return; } if (c == 1) return; for (i = 1; i <= (da+1); i++) { A[i] = (A[i] * c) % M; } for( ; (da >= 0) && (A[da+1] == -1); da --); A[0] = da; return; } /* A := A * c. Multiply A by an scalar c (an element of the number field) */ void ip_scale_mul(poly A, poly c, int N, array Si, poly E, LONG M, longarray WS) { int da, i, sa, block_size, dm; poly temp; da = A[0]; if (da == -1) return; if (c[0] == -1) { A[0] = -1; return; } if (N == 0) return ip_scale_mul_no_ext(A, c[1], M); dm = E[0]; block_size = Si[N]; temp = WS; WS += space_d_N(2*dm - 2,N - 1, Si); for (i = 0, sa = 1; i <= da; i++, sa += block_size) { ip_mul_mod(A+sa, c, temp, N, Si, E, M, WS); copy_into(A+sa, temp, N - 1, Si); } for (sa -= block_size; (sa >= 1) && (A[sa] == -1); sa -= block_size, da--); A[0] = da; return; } /* Check if A is 1 or not. */ char is_one(poly A, int N) { if (N == 0) { if ((A[0] == 0) && (A[1] == 1)) return 1; return 0; } return (A[0] == 0) && (is_one(A+1, N - 1)); } /* Computing the inverse of an integer a modulu m. */ LONG ip_int_inv(LONG a, LONG m) { LONG r1, r2, s1, s2, t, q; if (a == 1) return 1; for (r1 = a, r2 = m, s1 = 1, s2 = 0; r2 != 0; ) { q = r1 / r2; t = r1 % r2; if (t < 0LL) t += m; r1 = r2; r2 = t; t = s1 - q*s2; s1 = s2; s2 = t; } s1 = s1 % m; if (s1 < 0LL) s1 += m; return s1; } /* R := A^(-1) mod m(z) where E contains m(z) and degree(m) = 2. */ int ip_inv_no_ext_quadratic(poly A, poly R, poly E, long M) { LONG in, denom, A2s, A1A2; if (A[0] == 0) { in = ip_int_inv(A[1], M); R[0] = 0; R[1] = in; return -1; } A2s = (A[2]*A[2]) % M; A1A2 = (A[1]*A[2]) % M; denom = ( (A2s*E[1]) % M - (A1A2*E[2]) % M + (A[1]*A[1]) ) % M; if (denom == 0) { R[0] = A[0]; R[1] = A[1]; R[2] = A[2]; return 0; } if (denom < 0LL) denom += M; denom = ip_int_inv(denom, M); R[0] = 1; A2s = (-denom*A[2]) % M; if (A2s < 0LL) A2s += M; R[2] = A2s; A2s = (-E[2]*A[2]) % M; A2s = (((A2s + A[1]) % M)*denom) % M; if (A2s < 0LL) A2s += M; R[1] = A2s; return -1; } /* R := A^(-1) mod m(z) where E contains m(z). */ int ip_inv_no_ext(poly A, poly R, array Si, poly E, long M, longarray WS) { poly r1; poly r2; poly s1; poly s2; LONG in, tmp; poly minpoly; poly temp; poly q; poly ttemp; int da, dm, dq, dr2; if ( is_one(A, 0) == 1 ) { R[0] = 0; R[1] = 1; return -1; } if (E[0] == 2) return ip_inv_no_ext_quadratic(A, R, E, M); minpoly = E; dm = minpoly[0]; da = A[0]; r1 = WS; WS += dm+2; r2 = WS; WS += dm+2; s1 = R; s2 = WS; WS += dm+2; ttemp = WS; WS += 2*dm; copy_into(r1, A, 0, Si); copy_into(r2, minpoly, 0, Si); s1[0] = 0; s2[0] = -1; if ( is_one(A, 0) == 1 ) { R[0] = 0; R[1] = 1; return -1; } if (E[0] == 2) return ip_inv_no_ext_quadratic(A, R, E, M); minpoly = E; dm = minpoly[0]; da = A[0]; r1 = WS; WS += dm+2; r2 = WS; WS += dm+2; s1 = R; s2 = WS; WS += dm+2; ttemp = WS; WS += 2*dm; copy_into(r1, A, 0, Si); copy_into(r2, minpoly, 0, Si); s1[0] = 0; s2[0] = -1; in = ip_int_inv(r1[ r1[0] + 1 ], M); ip_scale_mul_no_ext(r1, in, M); s1[1] = in; dr2 = r2[0]; for (; dr2 != -1; ) { in = ip_int_inv(r2[ r2[0] + 1], M); ip_scale_mul_no_ext(r2, in, M); ip_scale_mul_no_ext(s2, in, M); dq = r1[0] - dr2; if (dq < 0) dq = -1; ip_rem_no_ext(r1, r2, M); q = r1 + dr2; temp = r1; r1 = r2; r2 = temp; tmp = q[0]; q[0] = dq; ip_mul_no_ext(q, s2, ttemp,M); ip_rem_no_ext(ttemp, E, M); ip_sub(s1, ttemp, 0 , Si, M); q[0] = tmp; temp = s1; s1 = s2; s2 = temp; dr2 = r2[0]; } if (is_one(r1, 0)) { copy_into(R, s1, 0, Si); return -1; } copy_into(R, r1, 0, Si); return 0; } int ip_inv(poly A, poly R, int N, array Si, poly E, LONG M, longarray WS) { poly r1; poly r2; poly s1; poly s2; poly minpoly; poly in; poly q; poly F; poly temp; LONG tmp; poly ttemp; int i,dm, dr2, block_size, NZ, dq; if (N == 0) return ip_inv_no_ext(A, R, Si, E, M, WS); if ( is_one(A, N) == 1 ) { copy_into(R, A, N, Si); return -1; } minpoly = E; dm = minpoly[0]; block_size = Si[N]; F = E + block_size + Si[N + 1]; r1 = WS; WS += space_d_N(dm, N,Si); r2 = WS; WS += space_d_N(dm, N,Si); s1 = R; s2 = WS; WS += space_d_N(dm, N,Si); in = WS; WS += space_d_N(dm, N,Si); ttemp = WS; WS += space_d_N(2*dm - 2, N, Si); copy_into(r1, A, N, Si); copy_into(r2, minpoly, N, Si); for (i = 0; i <= N; i++) s1[i] = 0; s1[N+1] = 1; s2[0] = -1; NZ = ip_inv(r1 + r1[0]*block_size + 1 , in, N - 1, Si, F, M, WS); if (NZ >= 0) { /* Zero Divisor in R_NZ */ copy_into(R, in, NZ, Si); return NZ; } ip_scale_mul(r1, in, N, Si, F, M, WS); ip_scale_mul(s1, in, N, Si, F, M, WS); dr2 = r2[0]; for (; dr2 != -1; ) { NZ = ip_inv(r2 + dr2*block_size + 1, in, N - 1, Si, F, M, WS); if (NZ >= 0) { /* Zero Divisor in R_NZ */ copy_into(R, in, NZ, Si); return NZ; } ip_scale_mul(r2, in, N, Si, F, M, WS); ip_scale_mul(s2, in, N, Si, F, M, WS); dq = r1[0] - dr2; if (dq < 0) dq = -1; ip_rem_monic_b(r1, r2, N, Si, F, M, WS); q = r1 + dr2*block_size; temp = r1; r1 = r2; r2 = temp; tmp = q[0]; q[0] = dq; ip_mul(q, s2, ttemp, N , Si, F, M, WS); ip_rem_monic_b(ttemp, E, N, Si, F, M, WS); ip_sub(s1, ttemp, N , Si, M); q[0] = tmp; temp = s1; s1 = s2; s2 = temp; dr2 = r2[0]; } if (is_one(r1, N)) { copy_into(R, s1, N, Si); return -1; } copy_into(R, r1, N, Si); return N; } /* A := gcd(A, B). A and B are two univeriate polynomial in Z_M[x]. */ void ip_gcd_no_ext(poly A, poly B, LONG M) { LONG inv; poly temp; poly r1; poly r2; int da, db; if (A[0] < B[0]) { r1 = B; r2 = A; } else { r1 = A; r2 = B; } da = r1[0]; db = r2[0]; if (db == -1) return; inv = ip_int_inv( r1[da + 1] , M); ip_scale_mul_no_ext(r1, inv, M); for (; db != -1; ) { inv = ip_int_inv(r2[db + 1], M); ip_scale_mul_no_ext(r2, inv, M); ip_rem_no_ext(r1,r2,M); temp = r1; r1 = r2; r2 = temp; db = r2[0]; } for (da = 0; da <= r1[0] + 1; da++) A[da] = r1[da]; return; } /* Main GCD routine. A := gcd(A, B). A and B are two univariate polynomial over an algebraic number field with multiple field extensions. N is the number of field extensions. Si[i] is S_i in the paper. S_i is the amount of space needed to store one element in R_i. Si[1] = 1. E contains the N minimal polynomials. M is the integer modulus. WS is the working storage. */ int ip_gcd(poly A, poly B, int N, array Si, poly E, LONG M, longarray WS) { int da, db, NZ, block_size; poly inv; poly temp; poly r1; poly r2; if (N == 0) { ip_gcd_no_ext(A, B, M); return -1; } if (A[0] < B[0]) { r1 = B; r2 = A; } else { r1 = A; r2 = B; } da = r1[0]; db = r2[0]; if (db == -1) return -1; inv = WS; if (N == 0) WS += 1; else WS += Si[N] + Si[N - 1]; block_size = Si[N]; NZ = ip_inv(r1 + da*block_size + 1 , inv, N - 1, Si, E, M, WS); if (NZ >= 0) { copy_into(A, inv, NZ, Si); return NZ; } ip_scale_mul(r1, inv, N, Si, E, M, WS); for( ;db != -1;) { NZ = ip_inv(r2 + db*block_size + 1 , inv, N - 1, Si, E, M, WS); if (NZ >= 0) { copy_into(A, inv, NZ, Si); return NZ; } ip_scale_mul(r2, inv, N, Si, E, M, WS); ip_rem_monic_b(r1, r2, N, Si, E, M, WS); temp = r1; r1 = r2; r2 = temp; db = r2[0]; } copy_into(A, r1, N, Si); return -1; } EXT_DECL ALGEB M_DECL maple_gcd_wrapper(const MKernelVector k, const ALGEB args) { int r; ALGEB _A = ((ALGEB *)args)[ 1 ]; poly A; ALGEB _B = ((ALGEB *)args)[ 2 ]; poly B; ALGEB _N = ((ALGEB *)args)[ 3 ]; int N; ALGEB _Si = ((ALGEB *)args)[ 4 ]; array Si; ALGEB _E = ((ALGEB *)args)[ 5 ]; poly E; ALGEB _M = ((ALGEB *)args)[ 6 ]; LONG M; ALGEB _WS = ((ALGEB *)args)[ 7 ]; longarray WS; A = (poly) RTableDataBlock(k,_A); B = (poly) RTableDataBlock(k,_B); N = MapleToInteger32(k, _N); Si = (array) RTableDataBlock(k,_Si); E = (poly) RTableDataBlock(k,_E); M = (LONG) MapleToM_INT(k, _M); WS = (longarray) RTableDataBlock(k,_WS); r = ip_gcd(A,B,N,Si,E,M,WS); return ToMapleInteger(k, r); } EXT_DECL ALGEB M_DECL maple_inv_wrapper(const MKernelVector k, const ALGEB args) { int r; ALGEB _A = ((ALGEB *)args)[ 1 ]; poly A; ALGEB _B = ((ALGEB *)args)[ 2 ]; poly R; ALGEB _N = ((ALGEB *)args)[ 3 ]; int N; ALGEB _Si = ((ALGEB *)args)[ 4 ]; array Si; ALGEB _E = ((ALGEB *)args)[ 5 ]; poly E; ALGEB _M = ((ALGEB *)args)[ 6 ]; LONG M; ALGEB _WS = ((ALGEB *)args)[ 7 ]; longarray WS; A = (poly) RTableDataBlock(k,_A); R = (poly) RTableDataBlock(k,_B); N = MapleToInteger32(k, _N); Si = (array) RTableDataBlock(k,_Si); E = (poly) RTableDataBlock(k,_E); M = (LONG) MapleToM_INT(k, _M); WS = (longarray) RTableDataBlock(k,_WS); r = ip_inv(A,R,N,Si,E,M,WS); return ToMapleInteger(k, r); } EXT_DECL ALGEB M_DECL maple_mul_wrapper(const MKernelVector k, const ALGEB args) { ALGEB _A = ((ALGEB *)args)[ 1 ]; poly A; ALGEB _B = ((ALGEB *)args)[ 2 ]; poly B; ALGEB _C = ((ALGEB *)args)[ 3 ]; poly C; ALGEB _N = ((ALGEB *)args)[ 4 ]; int N; ALGEB _Si = ((ALGEB *)args)[ 5 ]; array Si; ALGEB _E = ((ALGEB *)args)[ 6 ]; poly E; ALGEB _M = ((ALGEB *)args)[ 7 ]; LONG M; ALGEB _WS = ((ALGEB *)args)[ 8 ]; longarray WS; A = (poly) RTableDataBlock(k,_A); B = (poly) RTableDataBlock(k,_B); C = (poly) RTableDataBlock(k,_C); N = MapleToInteger32(k, _N); Si = (array) RTableDataBlock(k,_Si); E = (poly) RTableDataBlock(k,_E); M = (LONG) MapleToM_INT(k, _M); WS = (longarray) RTableDataBlock(k,_WS); ip_mul(A,B,C,N,Si,E,M,WS); return ToMapleInteger(k, 0); } EXT_DECL ALGEB M_DECL maple_div_wrapper(const MKernelVector k, const ALGEB args) { ALGEB _A = ((ALGEB *)args)[ 1 ]; poly A; ALGEB _B = ((ALGEB *)args)[ 2 ]; poly B; ALGEB _N = ((ALGEB *)args)[ 3 ]; int N; ALGEB _Si = ((ALGEB *)args)[ 4 ]; array Si; ALGEB _E = ((ALGEB *)args)[ 5 ]; poly E; ALGEB _M = ((ALGEB *)args)[ 6 ]; LONG M; ALGEB _WS = ((ALGEB *)args)[ 7 ]; longarray WS; A = (poly) RTableDataBlock(k,_A); B = (poly) RTableDataBlock(k,_B); N = MapleToInteger32(k, _N); Si = (array) RTableDataBlock(k,_Si); E = (poly) RTableDataBlock(k,_E); M = (LONG) MapleToM_INT(k, _M); WS = (longarray) RTableDataBlock(k,_WS); ip_rem_monic_b(A,B,N,Si,E,M,WS); return ToMapleInteger(k, 0); }