Chinese Reminder for non co primes

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;

#define SZ 1005001
int phi[SZ];
int factor[SZ];
char isnp[SZ];

void initPhi(){
    for ( int i = 0 ; i < SZ ; i++ ) phi[i] = factor[i] = i;
    for ( int i = 4 ; i < SZ ; i += 2 ) isnp[i] = 1, factor[i] = 2;

    for ( int i = 9 ; i < SZ ; i += 3 ) isnp[i] = 1, factor[i] = 3;

    int nxt = 4;
    for ( int i = 5 ; i * i < SZ ; ){
        if ( isnp[i] == 0 ) {
            for ( int j = i * i ; j < SZ ; j += i )
                isnp[j] = 1, factor[j] = i;
        }

        nxt ^= 6;
        i += nxt;
    }

    cout << "done" << endl;
    for ( int i = 2 ; i < SZ ; i++ ){
        if ( isnp[i] == 0 ){
            for ( int j = i ; j < SZ ; j += i ){
                phi[j] -= phi[j]/i;
            }
        }
    }

    cout << "done" << endl;
}

// ax + by = g
// a'x' + b'y' = g , a'= b, b' = a - (a/b) * b
// -> bx' + ay' - (a/b)*by' = g
// -> ay' + b(x' - (a/b) * y') = g
// -> x = y', y = x' - (a/b) * y'
void extendedEuclid(int a, int b, int &x, int &y, int &g){
    int x1, y1;

    if ( b == 0 ){
        x = 1, y = 0, g = a;
        return;
    }

    int quo = a / b;
    extendedEuclid(b, a - quo * b, x1, y1, g);
    x = y1;
    y = x1 - quo * y1;
}

// ax = 1 mod b -> ax + by = 1
// if x < 0, ab + ax = 1 mod b -> a(b + x) = 1 mod b
int modInv(int a, int b){
    int x, y, g;

    a %= b;
    extendedEuclid(a, b, x, y, g);
    if ( x < 0 ) x = b + x;

    //cout << "ModInv(" << a << "," << b <<") : " << x << endl;
    return x;
}

ll multMod(ll a, ll b, ll m){
    ll res = 0;

    while ( b ){
        if ( b & 1 ) res = ( res + a );
        if ( res >= m ) res -= m;
        a <<= 1;
        if ( a >= m ) a -= m;
        b >>= 1;
    }

    return res;
}

ll ChineseRemainder(vector<int> &residue, vector<int> &mods){
    ll M = 1;
    ll res = 0;

    for ( int i = 0 ; i < mods.size() ; i++ ){
        M *= mods[i];
    }

    for ( int i = 0 ; i < residue.size() ; i++ ){
        ll temp = residue[i] * modInv(M / mods[i], mods[i]);
        temp = multMod(temp, M / mods[i], M);
        res = res + temp;
        if ( res >= M ) res -= M;
    }

    return res;
}

map<int, int> factors[5001];
// n is assumed to be between 1000000 and 1005000
map<int, int> getFactors(int n){
    int offset = 1000000;
    if ( factors[n - offset].size() == 0 ){
        int index = n - offset;
        int m = n;
        while ( m > 1 ) {
            int v = 0;
            int f = factor[m];
            while ( m % f == 0 ) v++, m /= f;
            factors[index].insert(pair<int, int>(f, v));
        }
    }

    return factors[n - offset];
}

//int powers[10005000][20];
int findPower(int b, int p) {
    int res = 1, tmp = b;

    while (p){
        if (p & 1) res = res * tmp;
        tmp *= tmp;
        p >>= 1;
    }

    return res;
}

long long chineseNonCoPrime(int n, int modn, int m, int modm) {
    map<int, int>::iterator ptrMap, ptrMap1;
    vector<int> res;
    vector<int> mods;
    map<int, int> fn = getFactors(n);
    map<int, int> fm = getFactors(m);

    for ( ptrMap = fn.begin(), ptrMap1 = fm.begin() ; ptrMap != fn.end() && ptrMap1 != fm.end() ; ) {
        if ( ptrMap->first == ptrMap1->first ){
            int p1 = findPower(ptrMap->first, ptrMap->second);
            int p2 = findPower(ptrMap1->first, ptrMap1->second);
            int m1 = modn % p1;
            int m2 = modm % p2;

            //printf("case 1 -> (%d mod %d), (%d mod %d)\n", m1, p1, m2, p2);
            int diff = ptrMap->second - ptrMap1->second;

            if ( diff == 0 && m1 != m2 ) return 0;

            if ( diff < 0 ){
                diff = -diff;
                if ( m2 != m1 ) m2 ^= m1 ^= m2 ^= m1;
                if ( p2 != p1 ) p2 ^= p1 ^= p2 ^= p1;
            }

            //printf("Now, case 1 -> (%d mod %d), (%d mod %d).... diff %d\n", m1, p1, m2, p2, p1/p2);
            if (m1 >= m2 && ((m1 - m2) % (p1/p2) == 0)) {
                res.push_back(m1);
                mods.push_back(p1);
            } else return 0;

            ptrMap++, ptrMap1++;
        } else if ( ptrMap->first < ptrMap1->first ) {
            int p1 = findPower(ptrMap->first, ptrMap->second);

            res.push_back(modn % p1);
            mods.push_back(p1);

            ptrMap++;
        } else {
            int p2 = findPower(ptrMap1->first, ptrMap1->second);

            res.push_back(modm % p2);
            mods.push_back(p2);

            ptrMap1++;
        }
    }

    while ( ptrMap != fn.end()) {
        int p1 = findPower(ptrMap->first, ptrMap->second);

        res.push_back(modn % p1);
        mods.push_back(p1);

        ptrMap++;
    }

    while ( ptrMap1 != fm.end() ){
        int p2 = findPower(ptrMap1->first, ptrMap1->second);

        res.push_back(modm % p2);
        mods.push_back(p2);

        ptrMap1++;
    }

    //printf("Found all res and mods, now going to remainder\n");
    return ChineseRemainder(res, mods);
}

int main(){
    int n, m, a, b;
    map<int, int>::iterator ptrMap;

    initPhi();


    // n is assumed to be between 1000000 and 1005000
    /*while ( cin >> n ){
        map<int, int> f = getFactors(n);
        for ( ptrMap = f.begin() ; ptrMap != f.end() ; ptrMap++ )
            cout << ptrMap->first << "^" << ptrMap->second << " ";
        cout << endl;
    }*/

    /*while ( cin >> a >> n >> b >> m ){
    //while ( cin >> n >> m ){
        cout << chineseNonCoPrime(n,a,m,b) << endl;
    }*/
    long long res = 0;
    for ( int n = 1000000 ; n < 1005000 ; n++ ){
        for ( int m = n + 1 ; m < 1005000 ; m++ )
            res += chineseNonCoPrime(n, phi[n], m, phi[m]);
        cout << n << " ended" << endl;
    }
    cout << res << endl;
    getchar();

    return 0;
}