+ if ( ++r == k ) {
+ r = 1;
+ ++u;
+ }
+ }
+}
+
+static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
+{
+ size_t newdeg = degree(a)/prime;
+ ap.resize(newdeg+1);
+ ap[0] = a[0];
+ for ( size_t i=1; i<=newdeg; ++i ) {
+ ap[i] = a[i*prime];
+ }
+}
+
+static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
+{
+ const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus);
+ int i = 1;
+ umodpoly b;
+ deriv(a, b);
+ if ( b.size() ) {
+ umodpoly c;
+ gcd(a, b, c);
+ umodpoly w;
+ div(a, c, w);
+ while ( unequal_one(w) ) {
+ umodpoly y;
+ gcd(w, c, y);
+ umodpoly z;
+ div(w, y, z);
+ factors.push_back(z);
+ mult.push_back(i);
+ ++i;
+ w = y;
+ umodpoly buf;
+ div(c, y, buf);
+ c = buf;
+ }
+ if ( unequal_one(c) ) {
+ umodpoly cp;
+ expt_1_over_p(c, prime, cp);
+ size_t previ = mult.size();
+ modsqrfree(cp, factors, mult);
+ for ( size_t i=previ; i<mult.size(); ++i ) {
+ mult[i] *= prime;
+ }
+ }
+ }
+ else {
+ umodpoly ap;
+ expt_1_over_p(a, prime, ap);
+ size_t previ = mult.size();
+ modsqrfree(ap, factors, mult);
+ for ( size_t i=previ; i<mult.size(); ++i ) {
+ mult[i] *= prime;
+ }
+ }
+}
+
+static void distinct_degree_factor(const umodpoly& a_, vector<int>& degrees, upvec& ddfactors)
+{
+ umodpoly a = a_;
+
+ cl_modint_ring R = a[0].ring();
+ int q = cl_I_to_int(R->modulus);
+ int nhalf = degree(a)/2;
+
+ int i = 1;
+ umodpoly w(2);
+ w[0] = R->zero();
+ w[1] = R->one();
+ umodpoly x = w;
+
+ bool nontrivial = false;
+ while ( i <= nhalf ) {
+ expt_pos(w, q);
+ umodpoly buf;
+ rem(w, a, buf);
+ w = buf;
+ umodpoly wx = w - x;
+ gcd(a, wx, buf);
+ if ( unequal_one(buf) ) {
+ degrees.push_back(i);
+ ddfactors.push_back(buf);
+ }
+ if ( unequal_one(buf) ) {
+ umodpoly buf2;
+ div(a, buf, buf2);
+ a = buf2;
+ nhalf = degree(a)/2;
+ rem(w, a, buf);
+ w = buf;
+ }
+ ++i;
+ }
+ if ( unequal_one(a) ) {
+ degrees.push_back(degree(a));
+ ddfactors.push_back(a);
+ }
+}
+
+static void same_degree_factor(const umodpoly& a, upvec& upv)
+{
+ cl_modint_ring R = a[0].ring();
+ int deg = degree(a);
+
+ vector<int> degrees;
+ upvec ddfactors;
+ distinct_degree_factor(a, degrees, ddfactors);
+
+ for ( size_t i=0; i<degrees.size(); ++i ) {
+ if ( degrees[i] == degree(ddfactors[i]) ) {
+ upv.push_back(ddfactors[i]);
+ }
+ else {
+ berlekamp(ddfactors[i], upv);