+ 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);
+ }
+ }
+}
+
+static void factor_modular(const umodpoly& p, upvec& upv)
+{
+ upvec factors;
+ vector<int> mult;
+ modsqrfree(p, factors, mult);
+
+#define USE_SAME_DEGREE_FACTOR
+#ifdef USE_SAME_DEGREE_FACTOR
+ for ( size_t i=0; i<factors.size(); ++i ) {
+ upvec upvbuf;
+ same_degree_factor(factors[i], upvbuf);
+ for ( int j=mult[i]; j>0; --j ) {
+ upv.insert(upv.end(), upvbuf.begin(), upvbuf.end());
+ }
+ }
+#else
+ for ( size_t i=0; i<factors.size(); ++i ) {
+ upvec upvbuf;
+ berlekamp(factors[i], upvbuf);
+ if ( upvbuf.size() ) {
+ for ( size_t j=0; j<upvbuf.size(); ++j ) {
+ upv.insert(upv.end(), mult[i], upvbuf[j]);
+ }
+ }
+ else {
+ for ( int j=mult[i]; j>0; --j ) {
+ upv.push_back(factors[i]);
+ }
+ }
+ }
+#endif
+}
+
+/** Calculates polynomials s and t such that a*s+b*t==1.
+ * Assertion: a and b are relatively prime and not zero.
+ *
+ * @param[in] a polynomial
+ * @param[in] b polynomial
+ * @param[out] s polynomial
+ * @param[out] t polynomial
+ */
+static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpoly& t)
+{
+ if ( degree(a) < degree(b) ) {
+ exteuclid(b, a, t, s);