+ 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