factorize.cpp

#include <vector>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

typedef long long Int;

Int gcd(Int a, Int b) {
  return b != 0 ? gcd(b, a % b) : a;
}

Int extgcd(Int a, Int b, Int &x, Int &y) {
    Int g = a; x = 1; y = 0;
    if (b != 0) g = extgcd(b, a % b, y, x), y -= (a / b) * x;
    return g;
}
Int invMod(Int a, Int m) {
    Int x, y;
    if (extgcd(a, m, x, y) == 1) return (x + m) % m;
    else                         return 0; // unsolvable
}
Int powMod(Int x, Int k, Int m) {
    if (k == 0)     return 1;
    if (k % 2 == 0) return powMod(x*x % m, k/2, m);
    else            return x*powMod(x, k-1, m) % m;
}

#define NEGPOW(e) ((e) % 2 ? -1 : 1)
Int jacobi(Int a, Int m) {
    if (a == 0) return m == 1 ? 1 : 0;
    if (a % 2)  return NEGPOW((a-1)*(m-1)/4)*jacobi(m%a, a);
    else        return NEGPOW((m*m-1)/8)*jacobi(a/2, m);
}

Int sqrtMod(Int n, Int p) {
    Int S, Q, W, i, m = invMod(n, p);
    for (Q = p - 1, S = 0; Q % 2 == 0; Q /= 2, ++S);
    do { W = rand() % p; } while (W == 0 || jacobi(W, p) != -1);
    for (Int R = powMod(n, (Q+1)/2, p), V = powMod(W, Q, p); ;) {
        Int z = R * R * m % p;
        for (i = 0; i < S && z % p != 1; z *= z, ++i);
        if (i == 0) return R;
        R = (R * powMod(V, 1 << (S-i-1), p)) % p;
    }
}

#define PRIME_MAX 80
void sieve(Int n, vector<Int>& primes, vector<Int>& psqrt) {
    int s[PRIME_MAX];
    for(int i = 0; i < PRIME_MAX; i++) {
        s[i] = 1;
    }
    for(int i = 2; i < PRIME_MAX; i++) {
        if(!s[i])continue;
        if(i==2) {
            psqrt.push_back(n%2);
            primes.push_back(i);
        } else if(jacobi(n,i)==1) {
            psqrt.push_back(sqrtMod(n,i));
            primes.push_back(i);
        }
        for(int j = i; j < PRIME_MAX; j+=i) {
            s[j] = 0;
        }
    }
}

int main(int argc, char **argv)
{
    Int n = 55751;
    vector<Int> primes;
    vector<Int> psqrt;
    primes.push_back(-1);
    psqrt.push_back(0);
    sieve(n,primes,psqrt);

    //for(int i = 0; i < (int)primes.size(); i++) {
    //    printf("(%lld,%lld), ", primes[i], psqrt[i]);
    //}
    //printf("\n");

    //pexp[i]は0か1で変動する
    vector<int> pexp(primes.size());
    pexp[3] = 1;

    Int f_a = 1;
    for(int i = 0; i < (int)pexp.size(); i++) {
        if(pexp[i]) {
            f_a *= primes[i];
        }
    }
    //f_bはpexp[i]==1なるiとその剰余平方根から中国人剰余定理による乗算で算出する
    Int f_b = 5;
    int fsize = 30;
    vector<Int> kfactor_p(fsize*2, 1);
    vector<vector<int> > kfactor(fsize*2, vector<int>(primes.size()));
    {
        // process -1 factor
        Int sqrtn = sqrt(n)/f_a;
        Int m1start = -sqrtn-f_b/f_a;
        Int m1end = sqrtn-f_b/f_a;
        for(Int x=m1start;x<=m1end;x++) {
            kfactor[x+fsize][0] = 1;
            kfactor_p[x+fsize] *= -1;
        }
    }
    for(int i = 1; i < (int)primes.size(); i++) {
        Int prime = primes[i];
        vector<Int> fs;
        fs.push_back(psqrt[i]);
        fs.push_back(prime-psqrt[i]);
        Int f_a_p = pexp[i] ? f_a/prime : f_a;
        Int primediv_p = pexp[i] ? prime : 1;
        Int prime_exp = prime;
        Int prime_exp_p = pexp[i] ? 1 : prime;
        bool flag = true;

        while(flag) {
            flag = false;
            sort(fs.begin(),fs.end());
            fs.resize(unique(fs.begin(),fs.end())-fs.begin());
            vector<Int> fss;
            vector<Int> fsn;
            for(int j = 0; j < (int)fs.size(); j++) {
                fss.push_back((fs[j]-f_b+prime_exp)*invMod(f_a_p,prime_exp)%prime_exp);
                if(prime!=2) {
                    fsn.push_back((
                                (n-fs[j]*fs[j])
                                *invMod(2,prime_exp*prime)
                                *invMod(fs[j],prime_exp*prime)
                                +fs[j]
                                )%(prime_exp*prime));
                }
            }
            sort(fss.begin(),fss.end());
            fss.resize(unique(fss.begin(),fss.end())-fss.begin());
            for(int j = 0; j < (int)fss.size(); j++) {
                Int kx = fss[j];
                if(kx%primediv_p!=0)continue;
                kx/=primediv_p;
                //printf("%lld:%lld(from %lld)\n",prime_exp,kx, fs[j]);
                Int start = (-fsize-kx)/prime_exp_p;
                Int end = (fsize-1-kx)/prime_exp_p;
                if(-fsize>start*prime_exp_p+kx)start++;
                if(end*prime_exp_p+kx>=fsize)end--;
                if(start<=end) {
                    flag = true;
                    for(Int k = start; k <= end; k++) {
                        kfactor[k*prime_exp_p+kx+fsize][i]++;
                        kfactor_p[k*prime_exp_p+kx+fsize] *= prime;
                    }
                }
            }
            prime_exp*=prime;
            prime_exp_p*=prime;
            fs=fsn;
        }
    }
    //for(Int x = -fsize; x < fsize; x++) {
    //    if(f_a*x*x+f_b*2*x+f_c==kfactor_p[x+fsize]/f_a)
    //    printf("%lld: %lld, %lld\n", x, f_a*x*x+f_b*2*x+f_c,kfactor_p[x+fsize]/f_a);
    //}
    vector<Int> xs;
    vector<vector<Int> > mat;
    for(Int x = -fsize; x < fsize; x++) {
        Int y = f_a*x+f_b;
        if(y*y-n==kfactor_p[x+fsize]) {
            printf("%lld: %lld, %lld\n", x, y*y-n,kfactor_p[x+fsize]);
            xs.push_back(x);
            mat.push_back(vector<Int>(primes.size()));
            for(int i = 0; i < (int)primes.size(); i++) {
                mat.back()[i] = kfactor[x+fsize][i]%2;
            }
        }
    }
    int mat_resize = primes.size()+mat.size();
    for(int i = 0; i < (int)mat.size(); i++) {
        mat[i].resize(mat_resize, 0);
        mat[i][primes.size()+i] = 1;
    }
    for(int i = 0; i < (int)primes.size(); i++) {
        int k = -1;
        for(int j = 0; j < (int)mat.size(); j++) {
            if(mat[j][i]==1) {
                k = j;
            }
        }
        for(int j = 0; j < (int)mat.size(); j++) {
            if(mat[j][i]==1) {
                for(int l = 0; l < (int)mat[j].size(); l++) {
                    mat[j][l] ^= mat[k][l];
                }
            }
        }
    }
    printf("processed matrix.\n");fflush(stdout);
    for(int i = 0; i < (int)mat.size(); i++) {
        bool flag=false;
        for(int j = 0; j < (int)mat.size(); j++) {
            if(mat[i][j+primes.size()]) {
                flag=true;
            }
        }
        if(!flag)continue;
        Int xprod = 1;
        vector<int> pfactor(primes.size());
        for(int j = 0; j < (int)mat.size(); j++) {
            if(mat[i][j+primes.size()]) {
                Int x = xs[j];
                xprod *= f_a*x+f_b;
                for(int k = 0; k < (int)kfactor[x+fsize].size(); k++) {
                    pfactor[k] += kfactor[x+fsize][k];
                }
            }
        }
        Int prod = 1;
        for(int j = 0; j < (int)pfactor.size(); j++) {
            for(int k = 0; k < pfactor[j]/2; k++) {
                prod *= primes[j];
            }
        }
        Int mgcd=gcd(xprod-prod,n);
        Int pgcd=gcd(xprod+prod,n);
        if(mgcd<0)mgcd=-mgcd;
        if(pgcd<0)pgcd=-pgcd;
        if(1<mgcd && mgcd<n) {
            printf("%lld\n",mgcd);
        }
        if(1<pgcd && pgcd<n) {
            printf("%lld\n",pgcd);
        }
        fflush(stdout);
    }
    return 0;
}