[863f53] | 1 | #include <stdio.h> |
---|
| 2 | #include <config.h> |
---|
[28f88cb] | 3 | #ifndef NOSTREAMIO |
---|
| 4 | #ifdef HAVE_IOSTREAM |
---|
| 5 | #include <iostream> |
---|
| 6 | #define OSTREAM std::ostream |
---|
| 7 | #elif defined(HAVE_IOSTREAM_H) |
---|
| 8 | #include <iostream.h> |
---|
| 9 | #define OSTREAM ostream |
---|
| 10 | #endif |
---|
| 11 | #endif /* NOSTREAMIO */ |
---|
[863f53] | 12 | |
---|
| 13 | #include "cf_defs.h" |
---|
| 14 | #include "canonicalform.h" |
---|
| 15 | #include "cf_iter.h" |
---|
| 16 | #include "cf_primes.h" |
---|
| 17 | #include "cf_algorithm.h" |
---|
| 18 | #include "algext.h" |
---|
| 19 | #include "fieldGCD.h" |
---|
| 20 | #include "cf_map.h" |
---|
| 21 | #include "cf_generator.h" |
---|
| 22 | |
---|
[27bb97f] | 23 | void out_cf(const char *s1,const CanonicalForm &f,const char *s2); |
---|
[13be01] | 24 | |
---|
| 25 | |
---|
[863f53] | 26 | CanonicalForm fieldGCD( const CanonicalForm & F, const CanonicalForm & G ); |
---|
| 27 | void CRA(const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew); |
---|
| 28 | |
---|
| 29 | |
---|
| 30 | CanonicalForm fieldGCD( const CanonicalForm & F, const CanonicalForm & G ) |
---|
| 31 | {// this is the modular method by Brown |
---|
| 32 | // assume F,G are multivariate polys over Z/p for big prime p |
---|
[13be01] | 33 | if(F.isZero()) |
---|
[863f53] | 34 | { |
---|
| 35 | if(G.isZero()) |
---|
| 36 | return G; // G is zero |
---|
| 37 | if(G.inCoeffDomain()) |
---|
| 38 | return CanonicalForm(1); |
---|
| 39 | return G/lc(G); // return monic G |
---|
| 40 | } |
---|
| 41 | if(G.isZero()) // F is non-zero |
---|
| 42 | { |
---|
| 43 | if(F.inCoeffDomain()) |
---|
| 44 | return CanonicalForm(1); |
---|
| 45 | return F/lc(F); // return monic F |
---|
| 46 | } |
---|
| 47 | if(F.inCoeffDomain() || G.inCoeffDomain()) |
---|
| 48 | return CanonicalForm(1); |
---|
[13be01] | 49 | //out_cf("F=",F,"\n"); |
---|
| 50 | //out_cf("G=",G,"\n"); |
---|
[863f53] | 51 | CFMap MM,NN; |
---|
| 52 | CFArray ps(1,2); |
---|
| 53 | ps[1] = F; |
---|
| 54 | ps[2] = G; |
---|
| 55 | compress(ps,MM,NN); // maps MM, NN are created |
---|
| 56 | CanonicalForm f=MM(F); |
---|
| 57 | CanonicalForm g=MM(G); |
---|
| 58 | // here: f, g are compressed |
---|
| 59 | // compute largest variable in f or g (least one is Variable(1)) |
---|
| 60 | int mv = f.level(); |
---|
| 61 | if(g.level() > mv) |
---|
| 62 | mv = g.level(); |
---|
| 63 | // here: mv is level of the largest variable in f, g |
---|
| 64 | if(mv == 1) // f,g univariate |
---|
| 65 | return NN(gcd( f, g )); // do not forget to map back |
---|
| 66 | // here: mv > 1 |
---|
| 67 | CanonicalForm cf = vcontent(f, Variable(2)); // cf is univariate poly in var(1) |
---|
| 68 | CanonicalForm cg = vcontent(g, Variable(2)); |
---|
| 69 | CanonicalForm c = gcd( cf, cg ); |
---|
| 70 | f/=cf; |
---|
| 71 | g/=cg; |
---|
| 72 | if(f.inCoeffDomain() || g.inCoeffDomain()) |
---|
[13be01] | 73 | { |
---|
| 74 | //printf("=============== inCoeffDomain\n"); |
---|
[863f53] | 75 | return NN(c); |
---|
[13be01] | 76 | } |
---|
[863f53] | 77 | int *L = new int[mv+1]; // L is addressed by i from 2 to mv |
---|
| 78 | int *M = new int[mv+1]; |
---|
| 79 | for(int i=2; i<=mv; i++) |
---|
| 80 | L[i] = M[i] = 0; |
---|
| 81 | L = leadDeg(f, L); |
---|
| 82 | M = leadDeg(g, M); |
---|
| 83 | CanonicalForm gamma = gcd( firstLC(f), firstLC(g) ); |
---|
| 84 | for(int i=2; i<=mv; i++) // entry at i=1 is not visited |
---|
| 85 | if(M[i] < L[i]) |
---|
| 86 | L[i] = M[i]; |
---|
| 87 | // L is now upper bound for leading degrees of gcd |
---|
| 88 | int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1 |
---|
| 89 | for(int i=2; i<=mv; i++) |
---|
| 90 | dg_im[i] = 0; // initialize |
---|
| 91 | CanonicalForm gamma_image, m=1; |
---|
| 92 | CanonicalForm gm=0; |
---|
| 93 | CanonicalForm g_image, alpha, gnew, mnew; |
---|
| 94 | FFGenerator gen = FFGenerator(); |
---|
| 95 | for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next()) |
---|
| 96 | { |
---|
| 97 | alpha = gen.item(); |
---|
| 98 | gamma_image = gamma(alpha, Variable(1)); // plug in alpha for var(1) |
---|
| 99 | if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha |
---|
| 100 | continue; |
---|
[13be01] | 101 | g_image = gcd( f(alpha, Variable(1)), g(alpha, Variable(1)) ); // recursive call with one var less |
---|
[863f53] | 102 | if(g_image.inCoeffDomain()) // early termination |
---|
[13be01] | 103 | { |
---|
| 104 | //printf("================== inCoeffDomain\n"); |
---|
[863f53] | 105 | return NN(c); |
---|
[13be01] | 106 | } |
---|
[863f53] | 107 | for(int i=2; i<=mv; i++) |
---|
| 108 | dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg) |
---|
| 109 | dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer |
---|
| 110 | if(isEqual(dg_im, L, 2, mv)) |
---|
| 111 | { |
---|
| 112 | g_image /= lc(g_image); // make g_image monic |
---|
| 113 | g_image *= gamma_image; // multiply by multiple of image lc(gcd) |
---|
| 114 | CRA( g_image, Variable(1)-alpha, gm, m, gnew, mnew ); |
---|
| 115 | // gnew = gm mod m |
---|
| 116 | // gnew = g_image mod var(1)-alpha |
---|
| 117 | // mnew = m * (var(1)-alpha) |
---|
| 118 | m = mnew; |
---|
| 119 | if(gnew == gm) // gnew did not change |
---|
| 120 | { |
---|
| 121 | g_image = gm / vcontent(gm, Variable(2)); |
---|
[28f88cb] | 122 | //out_cf("=========== try ",g_image,"\n"); |
---|
[863f53] | 123 | if(fdivides(g_image,f) && fdivides(g_image,g)) // trial division |
---|
[28f88cb] | 124 | { |
---|
| 125 | //printf("=========== okay\n"); |
---|
[863f53] | 126 | return NN(c*g_image); |
---|
[28f88cb] | 127 | } |
---|
[863f53] | 128 | } |
---|
| 129 | gm = gnew; |
---|
| 130 | continue; |
---|
| 131 | } |
---|
| 132 | if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky |
---|
| 133 | continue; |
---|
| 134 | |
---|
[ad8e1b] | 135 | // here: dg_im < L --> all previous points were unlucky |
---|
| 136 | //printf("=========== reset\n"); |
---|
| 137 | m = CanonicalForm(1); // reset |
---|
| 138 | gm = 0; // reset |
---|
| 139 | for(int i=2; i<=mv; i++) // tighten bound |
---|
| 140 | L[i] = dg_im[i]; |
---|
[863f53] | 141 | } |
---|
| 142 | // hopefully, we never reach this point |
---|
| 143 | Off( SW_USE_fieldGCD ); |
---|
| 144 | g_image = gcd(f,g); |
---|
| 145 | On( SW_USE_fieldGCD ); |
---|
| 146 | return g_image; |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | |
---|
| 150 | void CRA(const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew) |
---|
| 151 | { // this is polynomial Chinese Remaindering. q1, q2 are assumed to be coprime. |
---|
| 152 | // polynomials of level <= 1 are considered coefficients |
---|
| 153 | // xnew = x1 mod q1 (coefficientwise in the above sense) |
---|
| 154 | // xnew = x2 mod q2 |
---|
| 155 | // qnew = q1*q2 |
---|
| 156 | if(x1.level() <= 1 && x2.level() <= 1) // base case |
---|
| 157 | { |
---|
| 158 | (void) extgcd(q1,q2,xnew,qnew); |
---|
| 159 | xnew = x1 + (x2-x1) * xnew * q1; |
---|
| 160 | qnew = q1*q2; |
---|
| 161 | xnew = mod(xnew,qnew); |
---|
| 162 | return; |
---|
| 163 | } |
---|
| 164 | CanonicalForm tmp,tmp2; |
---|
| 165 | xnew = 0; |
---|
| 166 | qnew = q1 * q2; |
---|
| 167 | // here: x1.level() > 1 || x2.level() > 1 |
---|
| 168 | if(x1.level() > x2.level()) |
---|
| 169 | { |
---|
| 170 | for(CFIterator i=x1; i.hasTerms(); i++) |
---|
| 171 | { |
---|
| 172 | if(i.exp() == 0) // const. term |
---|
| 173 | { |
---|
| 174 | CRA(i.coeff(),q1,x2,q2,tmp,tmp2); |
---|
| 175 | xnew += tmp; |
---|
| 176 | } |
---|
| 177 | else |
---|
| 178 | { |
---|
| 179 | CRA(i.coeff(),q1,0,q2,tmp,tmp2); |
---|
| 180 | xnew += tmp * power(x1.mvar(),i.exp()); |
---|
| 181 | } |
---|
| 182 | } |
---|
| 183 | return; |
---|
| 184 | } |
---|
| 185 | // here: x1.level() <= x2.level() && ( x1.level() > 1 || x2.level() > 1 ) |
---|
| 186 | if(x2.level() > x1.level()) |
---|
| 187 | { |
---|
| 188 | for(CFIterator j=x2; j.hasTerms(); j++) |
---|
| 189 | { |
---|
| 190 | if(j.exp() == 0) // const. term |
---|
| 191 | { |
---|
| 192 | CRA(x1,q1,j.coeff(),q2,tmp,tmp2); |
---|
| 193 | xnew += tmp; |
---|
| 194 | } |
---|
| 195 | else |
---|
| 196 | { |
---|
| 197 | CRA(0,q1,j.coeff(),q2,tmp,tmp2); |
---|
| 198 | xnew += tmp * power(x2.mvar(),j.exp()); |
---|
| 199 | } |
---|
| 200 | } |
---|
| 201 | return; |
---|
| 202 | } |
---|
| 203 | // here: x1.level() == x2.level() && x1.level() > 1 && x2.level() > 1 |
---|
| 204 | CFIterator i = x1; |
---|
| 205 | CFIterator j = x2; |
---|
| 206 | while(i.hasTerms() || j.hasTerms()) |
---|
| 207 | { |
---|
| 208 | if(i.hasTerms()) |
---|
| 209 | { |
---|
| 210 | if(j.hasTerms()) |
---|
| 211 | { |
---|
| 212 | if(i.exp() == j.exp()) |
---|
| 213 | { |
---|
| 214 | CRA(i.coeff(),q1,j.coeff(),q2,tmp,tmp2); |
---|
| 215 | xnew += tmp * power(x1.mvar(),i.exp()); |
---|
| 216 | i++; j++; |
---|
| 217 | } |
---|
| 218 | else |
---|
| 219 | { |
---|
| 220 | if(i.exp() < j.exp()) |
---|
| 221 | { |
---|
| 222 | CRA(i.coeff(),q1,0,q2,tmp,tmp2); |
---|
| 223 | xnew += tmp * power(x1.mvar(),i.exp()); |
---|
| 224 | i++; |
---|
| 225 | } |
---|
| 226 | else // i.exp() > j.exp() |
---|
| 227 | { |
---|
| 228 | CRA(0,q1,j.coeff(),q2,tmp,tmp2); |
---|
| 229 | xnew += tmp * power(x1.mvar(),j.exp()); |
---|
| 230 | j++; |
---|
| 231 | } |
---|
| 232 | } |
---|
| 233 | } |
---|
| 234 | else // j is out of terms |
---|
| 235 | { |
---|
| 236 | CRA(i.coeff(),q1,0,q2,tmp,tmp2); |
---|
| 237 | xnew += tmp * power(x1.mvar(),i.exp()); |
---|
| 238 | i++; |
---|
| 239 | } |
---|
| 240 | } |
---|
| 241 | else // i is out of terms |
---|
| 242 | { |
---|
| 243 | CRA(0,q1,j.coeff(),q2,tmp,tmp2); |
---|
| 244 | xnew += tmp * power(x1.mvar(),j.exp()); |
---|
| 245 | j++; |
---|
| 246 | } |
---|
| 247 | } |
---|
| 248 | } |
---|