// Proof-of-concept implementation of factoring into coprimes // 2004.04.07, Daniel J. Bernstein // // To quote the paper: ``Warning: The algorithms in this paper are _not_ // designed to be as fast as possible. They are designed to be as _simple_ // as possible, under the constraints of (1) running in essentially linear // time and (2) using a limited set of arithmetic operations. There are // many ways to save more time; these speedups are important but are beyond // the scope of this paper.'' // // Sample use: echo '[91 119 221 1547 6898073]' | ./dcba-concept #include #include #include using namespace LiDIA; using namespace std; /* a^2^n */ bigint twopower(bigint a,unsigned int n) { /* Algorithm 10.1 */ for (;;) { if (!n) return a; a = a * a; n = n - 1; } } /* gout = gcd{a,c}; x = ppi(a,c); y = ppo(a,c) */ void gcdppippo(bigint &gout,bigint &x,bigint &y,bigint a,bigint c) { /* Algorithm 11.3 */ bigint g; x = gcd(a,c); gout = x; y = a / x; for (;;) { g = gcd(x,y); if (g == 1) return; x = x * g; y = y / g; } } /* x = ppi(a,c); y = ppo(a,c) */ void ppippo(bigint &x,bigint &y,bigint a,bigint c) { bigint g; gcdppippo(g,x,y,a,c); } /* ppi(a,c) */ bigint ppi(bigint a,bigint c) { bigint g; bigint x; bigint y; gcdppippo(g,x,y,a,c); return x; } /* gout = gcd{a,b}; x = ppg(a,b); y = pple(a,b) */ void gcdppgpple(bigint &gout,bigint &x,bigint &y,bigint a,bigint b) { /* Algorithm 11.4 */ bigint g; y = gcd(a,b); gout = y; x = a / y; for (;;) { g = gcd(x,y); if (g == 1) return; x = x * g; y = y / g; } } /* out = out + cb{a,b} */ void append_cb(base_vector &out,bigint a,bigint b) { /* Algorithm 13.2 */ if (b == 1) { if (a != 1) out[out.size()] = a; return; } bigint r; ppippo(a,r,a,b); if (r != 1) out[out.size()] = r; bigint g; bigint h; bigint c; gcdppgpple(g,h,c,a,b); bigint c0 = c; bigint x = c0; unsigned int n = 1; for (;;) { gcdppgpple(g,h,c,h,g * g); bigint d = gcd(c,b); x = x * d; bigint y = twopower(d,n - 1); append_cb(out,c / y,d); if (h == 1) break; n = n + 1; } append_cb(out,b / x,c0); } /* prod S */ bigint prod(base_vector S) { /* Algorithm 14.1 */ unsigned int n = S.size(); if (!n) return 1; if (n == 1) return S[0]; base_vector T(0,EXPAND); T.assign(0,S,0,n/2 - 1); bigint X = prod(T); T.assign(0,S,n/2,n - 1); bigint Y = prod(T); return X * Y; } /* append each ppi(a,P[i]) to out */ void splitappend(base_vector &out,bigint a,base_vector P) { /* Algorithm 15.3 */ unsigned int n = P.size(); if (!n) return; bigint b = ppi(a,prod(P)); if (n == 1) { out[out.size()] = b; return; } base_vector Q(0,EXPAND); Q.assign(0,P,0,n/2 - 1); splitappend(out,b,Q); Q.assign(0,P,n/2,n - 1); splitappend(out,b,Q); } /* split(a,P) */ base_vector split(bigint a,base_vector P) { base_vector result(0,EXPAND); splitappend(result,a,P); return result; } /* cb(P union {b}) */ base_vector cbextend(base_vector P,bigint b) { /* Algorithm 16.2 */ base_vector result(0,EXPAND); if (!P.size()) { if (b != 1) result[0] = b; return result; } bigint x = prod(P); bigint a; bigint r; ppippo(a,r,b,x); if (r != 1) result[0] = r; base_vector S(0,EXPAND); S = split(a,P); for (unsigned int i = 0;i < P.size();++i) append_cb(result,P[i],S[i]); return result; } /* bit_i k */ int bit(unsigned int i,unsigned int k) { if (k & (1 << i)) return 1; return 0; } /* cb(P union Q) */ base_vector cbmerge(base_vector P,base_vector Q) { /* Algorithm 17.3 */ base_vector S(0,EXPAND); base_vector T(0,EXPAND); base_vector R(0,EXPAND); unsigned int n = Q.size(); unsigned int b = 1; while ((1 << b) < n) ++b; S = P; unsigned int i = 0; for (;;) { if (i == b) return S; R.set_size(0); for (unsigned int k = 0;k < n;++k) if (bit(i,k) == 0) R[R.size()] = Q[k]; T = cbextend(S,prod(R)); R.set_size(0); for (unsigned int k = 0;k < n;++k) if (bit(i,k) == 1) R[R.size()] = Q[k]; S = cbextend(T,prod(R)); ++i; } } /* cb S */ base_vector cb(base_vector S) { /* Algorithm 18.1 */ base_vector result(0,EXPAND); unsigned int n = S.size(); if (n == 0) return result; if (n == 1) { if (S[0] != 1) result[0] = S[0]; return result; } base_vector T(0,EXPAND); T.assign(0,S,0,n/2 - 1); base_vector P = cb(T); T.assign(0,S,n/2,n - 1); base_vector Q = cb(T); return cbmerge(P,Q); } /* (i,pai) = reduce(p,a) */ void reduce(bigint &i,bigint &pai,bigint p,bigint a) { /* Algorithm 19.2 */ if (!!(a % p)) { i = 0; pai = a; return; } bigint j; bigint b; reduce(j,b,p * p,a / p); if (!(b % p)) { i = 2 * j + 2; pai = b / p; return; } i = 2 * j + 1; pai = b; } int printfactors(bigint a,base_vector P) { /* Algorithm 20.1 */ unsigned int n = P.size(); if (!n) { if (a != 1) { cout << "fail"; return 0; } return 1; } if (n == 1) { bigint m; bigint c; reduce(m,c,P[0],a); if (c != 1) { cout << "fail"; return 0; } cout << P[0] << "^" << m << " "; return 1; } base_vector Q(0,EXPAND); Q.assign(0,P,0,n/2 - 1); bigint y = prod(Q); bigint b; bigint c; ppippo(b,c,a,y); if (!printfactors(b,Q)) return 0; Q.assign(0,P,n/2,n - 1); if (!printfactors(c,Q)) return 0; return 1; } void printfactors(base_vector S,base_vector P) { /* Algorithm 21.2 */ unsigned int n = S.size(); if (!n) return; bigint x = prod(P); bigint y = prod(S); bigint z = ppi(x,y); base_vector D = split(z,P); base_vector Q(0,EXPAND); for (unsigned int i = 0;i < P.size();++i) if (D[i] == P[i]) Q[Q.size()] = P[i]; if (n == 1) { cout << y << " = "; printfactors(y,Q); cout << "\n"; return; } base_vector T(0,EXPAND); T.assign(0,S,0,n/2 - 1); printfactors(T,Q); T.assign(0,S,n/2,n - 1); printfactors(T,Q); } int main() { base_vector S; base_vector P; cin >> S; P = cb(S); cout << "S = " << S << "\n"; cout << "cb S = " << P << "\n"; printfactors(S,P); return 0; }