#include <iostream> using namespace std; #include <math.h> #ifdef _M_IX86 #define umulrem(z, x, y, m) \ __asm mov eax, x \ __asm mul y \ __asm div m \ __asm mov z, edx #define umuladdrem(z, x, y, a, m) \ __asm mov eax, x \ __asm mul y \ __asm add eax, a \ __asm adc edx, 0 \ __asm div m \ __asm mov z, edx #else #ifdef _MSC_VER typedef unsigned __int64 Tu64; #else typedef unsigned long long Tu64; #endif #define umulrem(z, x, y, m) \ { \ z = (unsigned int)(x * (Tu64)y % m); \ } #define umuladdrem(z, x, y, a, m) \ { \ z = (unsigned int)((x * (Tu64)y + a) % m); \ } #endif static bool IsPrime(unsigned int n) { if (n < 2) return false; if (n < 4) return true; if (n % 2 == 0) return false; const unsigned int iMax = (int)sqrt(n) + 1; unsigned int i; for (i = 3; i <= iMax; i += 2) if (n % i == 0) return false; return true; } static unsigned int LargestPrimeFactor(unsigned int n) { if (n < 2) return 1; unsigned int r = n, p; if (r % 2 == 0) { p = 2; do { r /= 2; } while (r % 2 == 0); } unsigned int i; for (i = 3; i <= r; i += 2) { if (r % i == 0) { p = i; do { r /= i; } while (r % i == 0); } } return p; } static unsigned int Powm(unsigned int n, unsigned int e, unsigned int m) { unsigned int r = 1; unsigned int t = n % m; unsigned int i; for (i = e; i != 0; i /= 2) { if (i % 2 != 0) { umulrem(r, r, t, m); } umulrem(t, t, t, m); } return r; } static unsigned int Inv(unsigned int n, unsigned int m) { unsigned int a = n, b = m; int u = 1, v = 0; do { const unsigned int q = a / b; const unsigned int t1 = a - q*b; a = b; b = t1; const int t2 = u - (int)q*v; u = v; v = t2; } while (b != 0); if (a != 1) u = 0; if (u < 0) u += m; return u; } class CPolyMod { protected: // (mod x^r - 1, n) const unsigned int m_r; const unsigned int m_n; unsigned int m_deg; unsigned int * mp_a; private: CPolyMod():m_r(0), m_n(0) { mp_a = NULL; }; public: // default value is x CPolyMod(unsigned int r, unsigned int n) : m_r(r), m_n(n) { m_deg = 1; mp_a = new unsigned int [2]; mp_a[0] = 0; mp_a[1] = 1; } CPolyMod(const CPolyMod & p) : m_r(p.m_r), m_n(p.m_n) { m_deg = p.m_deg; mp_a = new unsigned int [p.m_deg + 1]; unsigned int i; for (i = 0; i <= p.m_deg; ++i) mp_a[i] = p.mp_a[i]; } virtual ~CPolyMod() { if (mp_a != NULL) delete [] mp_a; } private: void _polySquare() { const unsigned int deg = m_deg; const unsigned int n = m_n; const unsigned int * const p_a = mp_a; const unsigned int degr = deg + deg; unsigned int * const p_ar = new unsigned int [degr + 1]; unsigned int k; for (k = 0; k <= degr; ++k) p_ar[k] = 0; unsigned int j; for (j = 1; j <= deg; ++j) { const unsigned int x = p_a[j]; if (x != 0) { unsigned int i; for (i = 0; i < j; ++i) { const unsigned int y = 2 * p_a[i]; unsigned int t = p_ar[j + i]; umuladdrem(t, x, y, t, n); p_ar[j + i] = t; } } } unsigned int i; for (i = 0; i <= deg; ++i) { const unsigned int x = p_a[i]; unsigned int t = p_ar[2 * i]; umuladdrem(t, x, x, t, n); p_ar[2 * i] = t; } m_deg = degr; delete [] mp_a; mp_a = p_ar; } void _polyMul(const CPolyMod & p) { const unsigned int deg = m_deg; const unsigned int n = m_n; const unsigned int * const p_a = mp_a; const unsigned int degr = deg + p.m_deg; unsigned int * const p_ar = new unsigned int [degr + 1]; unsigned int k; for (k = 0; k <= degr; ++k) p_ar[k] = 0; unsigned int j; for (j = 0; j <= p.m_deg; ++j) { const unsigned int x = p.mp_a[j]; if (x != 0) { unsigned int i; for (i = 0; i <= deg; ++i) { const unsigned int y = p_a[i]; unsigned int t = p_ar[j + i]; umuladdrem(t, x, y, t, n); p_ar[j + i] = t; } } } m_deg = degr; delete [] mp_a; mp_a = p_ar; } void _Mod() { unsigned int deg = m_deg; unsigned int * const p_a = mp_a; while (deg >= m_r) { p_a[deg - m_r] += p_a[deg]; if (p_a[deg - m_r] >= m_n) p_a[deg - m_r] -= m_n; --deg; while (p_a[deg] == 0) --deg; } m_deg = deg; } void _Norm() { const unsigned int deg = m_deg; const unsigned int n = m_n; unsigned int * const p_a = mp_a; if (p_a[deg] != 1) { const unsigned int y = Inv(p_a[deg], m_n); unsigned int i; for (i = 0; i <= deg; ++i) { unsigned int t = p_a[i]; umulrem(t, t, y, n); p_a[i] = t; } } } public: CPolyMod & operator = (const CPolyMod & p) { if (&p == this) return *this; m_deg = p.m_deg; delete [] mp_a; mp_a = new unsigned int [p.m_deg + 1]; unsigned int i; for (i = 0; i <= p.m_deg; ++i) mp_a[i] = p.mp_a[i]; return *this; } int operator != (const CPolyMod & p) const { if (m_deg != p.m_deg) return true; unsigned int i; for (i = 0; i <= m_deg; ++i) if (mp_a[i] != p.mp_a[i]) return true; return false; } CPolyMod & operator += (unsigned int i) { const unsigned int t = i % m_n; mp_a[0] += t; if (mp_a[0] >= m_n) mp_a[0] -= m_n; return *this; } CPolyMod & operator -= (unsigned int i) { const unsigned int t = m_n - i % m_n; mp_a[0] += t; if (mp_a[0] >= m_n) mp_a[0] -= m_n; return *this; } CPolyMod Pow(unsigned int e) const { unsigned int er = 1; unsigned int j; for (j = e; j != 1; j /= 2) { er = 2 * er + (j % 2); } CPolyMod t(*this); unsigned int i; for (i = er; i != 1; i /= 2) { t._polySquare(); t._Mod(); if (i % 2 != 0) { t._polyMul(*this); t._Mod(); } } t._Norm(); return t; } }; int main(int argc, char * argv[]) { // AKS // n=11701, r=11699, q=5849, s=2923 // n=1000000007, r=57287, q=28643, s=14311 // // Bernstein // n=349, r=347, q=173, s=140 // n=1000000007, r=3623, q=1811, s=1785 unsigned int n0; cout << "n ? "; cin >> n0; unsigned int n; for (n = n0; true; n += 2) { bool b = false; unsigned int q, r, s; for (r = 3; r < n; r += 2) { if (IsPrime(r)) { const unsigned int m = n % r; if (m == 0) break; q = LargestPrimeFactor(r - 1); if (Powm(m, (r - 1) / q, r) <= 1) continue; /* // AKS if (q >= 4 * sqrt(r) * n.Log()/log(2)) { s = (unsigned int)(2 * sqrt(r) * n.Log()/log(2)); b = true; break; } */ // Bernstein const double cMin = 2 * floor(sqrt(r)) * log(n); double c = 0; for (s = 1; s <= q; s++) { c += log(q + s - 1) - log(s); if (c >= cMin) { b = true; break; } } if (b) break; } } if (b) { cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << "\r" << flush;; bool b = true; CPolyMod rPoly(r, n); // x try { rPoly = rPoly.Pow(n); // x^n } catch (...) { b = false; } if (b) { unsigned int a; for (a = 1; a <= s; ++a) { cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << " " << a << "\r" << flush;; CPolyMod lPoly(r, n); // x lPoly -= a; // x - a try { lPoly = lPoly.Pow(n); // (x - a)^n } catch (...) { b = false; break; } lPoly += a; // (x - a)^n + a if (lPoly != rPoly) { b = false; break; } } } if (b) cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << ", power of a prime." << endl; else cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << ", composite." << endl; } } return 0; }