/*The magic pairs problem:

Let SumFact(n) be the sum of factors
of n.

Find all n1,n2 in a range such that

SumFact(n1)-n1-1 == n2  and
SumFact(n2)-n2-1 == n1

-----------------------------------------------------
To find SumFact(k), start with prime factorization:

k=(p1^n1)(p2^n2) ... (pN^nN)

THEN,

SumFact(k)=(1+p1+p1^2...p1^n1)*(1+p2+p2^2...p2^n2)*
(1+pN+pN^2...pN^nN)

PROOF:

Do a couple examples -- it's obvious:

48=2^4*3

SumFact(48)=(1+2+4+8+16)*(1+3)=1+2+4+8+16+3+6+12+24+48

75=3*5^2

SumFact(75)=(1+3)*(1+5+25)      =1+5+25+3+15+75

Corollary:

SumFact(k)=SumFact(p1^n1)*SumFact(p2^n2)*...*SumFact(pN^nN)

*/

//Primes are needed to sqrt(N).  Therefore, we can use U32.
class PowPrime
{
    I64 n;
    I64 sumfact; //Sumfacts for powers of primes are needed beyond sqrt(N)
};

class Prime
{
    U32          prime, pow_count;
    PowPrime    *pp;
};

I64 *PrimesNew(I64 N, I64 *_sqrt_primes, I64 *_cbrt_primes)
{
    I64      i, j, sqrt = Ceil(Sqrt(N)), cbrt = Ceil(N ` (1 / 3.0)),
             sqrt_sqrt = Ceil(Sqrt(sqrt)), sqrt_primes = 0, cbrt_primes = 0;
    U8      *s = CAlloc((sqrt + 1 + 7) / 8);
    Prime   *primes, *p;

    for (i = 2; i <= sqrt_sqrt; i++)
    {
        if (!Bt(s, i))
        {
            j = i * 2;
            while (j <= sqrt)
            {
                Bts(s, j);
                j += i;
            }
        }
    }
    for (i = 2; i <= sqrt; i++)
        if (!Bt(s, i))
        {
            sqrt_primes++; //Count primes
            if (i <= cbrt)
                cbrt_primes++;
        }

    p = primes = CAlloc(sqrt_primes * sizeof(Prime));
    for (i = 2; i <= sqrt; i++)
        if (!Bt(s, i))
        {
            p->prime = i;
            p++;
        }
    Free(s);

    *_sqrt_primes = sqrt_primes;
    *_cbrt_primes = cbrt_primes;

    return primes;
}

PowPrime *PowPrimesNew(I64 N, I64 sqrt_primes, Prime *primes, I64 *_num_powprimes)
{
    I64          i, j, k, sf, num_powprimes = 0;
    Prime       *p;
    PowPrime    *powprimes, *pp;

    p = primes;
    for (i = 0; i < sqrt_primes; i++)
    {
        num_powprimes += Floor(Ln(N) / Ln(p->prime));
        p++;
    }

    p = primes;
    pp = powprimes = MAlloc(num_powprimes * sizeof(PowPrime));
    for (i = 0; i < sqrt_primes; i++)
    {
        p->pp = pp;
        j = p->prime;
        k = 1;
        sf = 1;
        while (j < N)
        {
            sf += j;
            pp->n = j;
            pp->sumfact = sf;
            j *= p->prime;
            pp++;
            p->pow_count++;
        }
        p++;
    }
    *_num_powprimes = num_powprimes;

    return powprimes;
}

I64 SumFact(I64 n, I64 sqrt_primes, Prime *p)
{
    I64          i, k, sf = 1;
    PowPrime    *pp;

    if (n < 2)
        return 1;
    for (i = 0; i < sqrt_primes; i++)
    {
        k = 0;
        while (!(n % p->prime))
        {
            n /= p->prime;
            k++;
        }
        if (k)
        {
            pp = p->pp + (k - 1);
            sf *= pp->sumfact;
            if (n == 1)
                return sf;
        }
        p++;
    }

    return sf * (1 + n); //Prime
}

Bool TestSumFact(I64 n, I64 target_sf, I64 sqrt_primes, I64 cbrt_primes, Prime *p)
{
    I64          i = 0, k, b, x1, x2;
    PowPrime    *pp;
    F64         disc;

    if (n < 2)
        return FALSE;
    while (i++ < cbrt_primes)
    {
        k = 0;
        while (!(n % p->prime))
        {
            n /= p->prime;
            k++;
        }
        if (k)
        {
            pp = p->pp + (k - 1);
            if (ModU64(&target_sf, pp->sumfact))
                return FALSE;
            if (n == 1)
            {
                if (target_sf == 1)
                    return TRUE;
                else
                    return FALSE;
            }
        }
        p++;
    }
/*  At this point we have three possible cases to test
1)n==p1                 ->sf==(1+p1)                ?
2)n==p1*p1          ->sf==(1+p1+p1^2)   ?
3)n==p1*p2          ->sf==(p1+1)*(p2+1) ?

*/
    if (1 + n == target_sf)
    {
        while (i++ < sqrt_primes)
        {
            k = 0;
            while (!(n % p->prime))
            {
                n /= p->prime;
                k++;
            }
            if (k)
            {
                pp = p->pp + (k - 1);
                if (ModU64(&target_sf, pp->sumfact))
                    return FALSE;
                if (n == 1)
                {
                    if (target_sf == 1)
                        return TRUE;
                    else
                        return FALSE;
                }
            }
            p++;
        }
        if (1 + n == target_sf)
            return TRUE;
        else
            return FALSE;
    }

    k  =Sqrt(n);
    if (k * k == n)
    {
        if (1 + k + n == target_sf)
            return TRUE;
        else
            return FALSE;
    }
    else
    {
// n==p1*p2 -> sf==(p1+1)*(p2+1) ?  where p1 != 1 && p2 != 1
        // if p1==1 || p2==1, it is FALSE because we checked a single prime above.

        // sf==(p1+1)*(n/p1+1)
        // sf==n+p1+n/p1+1
        // sf*p1==n*p1+p1^2+n+p1
        // p1^2+(n+1-sf)*p1+n=0
        // x=(-b+/-sqrt(b^2-4ac))/2a
        // a=1
        // x=(-b+/-sqrt(b^2-4c))/2
        // b=n+1-sf;c=n
        b = n + 1 - target_sf;
// x=(-b+/-sqrt(b^2-4n))/2
        disc = b * b - 4 * n;
        if (disc < 0)
            return FALSE;
        x1 = (-b - Sqrt(disc)) / 2;
        if (x1 <= 1)
            return FALSE;
        x2 = n  / x1;
        if (x2 > 1 && x1 * x2 == n)
            return TRUE;
        else
            return FALSE;
    }
}

U0 PutFactors(I64 n) //For debugging
{
    I64 i, k, sqrt = Ceil(Sqrt(n));
    for (i = 2; i <= sqrt; i++)
    {
        k = 0;
        while (!(n % i))
        {
            k++;
            n /= i;
        }
        if (k)
        {
            "%d", i;
            if (k > 1)
                "^%d", k;
            '' CH_SPACE;
        }
    }
    if (n != 1)
        "%d ", n;
}

class RangeJob
{
    CDoc    *doc;
    I64      num, lo, hi, N, sqrt_primes, cbrt_primes;
    Prime   *primes;
    CJob    *cmd;

} rj[mp_count];

I64 TestCoreSubRange(RangeJob *r)
{
    I64          i, j, m, n, n2, sf, res = 0, range = r->hi - r->lo, 
                *sumfacts   = MAlloc(range * sizeof(I64)), 
                *residue    = MAlloc(range * sizeof(I64));
    U16         *pow_count  = MAlloc(range * sizeof(U16));
    Prime       *p = r->primes;
    PowPrime    *pp;

    MemSetI64(sumfacts, 1, range);
    for (n = r->lo; n < r->hi; n++)
        residue[n - r->lo] = n;
    for (j = 0; j <r->sqrt_primes; j++)
    {
        MemSet(pow_count, 0, range * sizeof(U16));
        m = 1;
        for (i = 0; i < p->pow_count; i++)
        {
            m *= p->prime;
            n = m - r->lo % m;
            while (n < range)
            {
                pow_count[n]++;
                n += m;
            }
        }
        for (n = 0; n < range; n++)
            if (i = pow_count[n])
            {
                pp = &p->pp[i - 1];
                sumfacts[n] *= pp->sumfact;
                residue [n] /= pp->n;
            }
        p++;
    }

    for (n = 0; n < range; n++)
        if (residue[n] != 1)
            sumfacts[n] *= 1 + residue[n];

    for (n = r->lo; n < r->hi; n++)
    {
        sf = sumfacts[n - r->lo];
        n2 = sf - n - 1;
        if (n < n2 < r->N)
        {
            if (r->lo <= n2<r->hi && sumfacts[n2 - r->lo] - n2 - 1 == n ||
                TestSumFact(n2, sf, r->sqrt_primes, r->cbrt_primes, r->primes))
            {
                DocPrint(r->doc, "%u:%u\n", n, sf - n - 1);
                res++;
            }
        }
    }
    Free(pow_count);
    Free(residue);
    Free(sumfacts);

    return res;
}

#define CORE_SUB_RANGE  0x1000

I64 TestCoreRange(RangeJob *r)
{
    I64      i, n, res = 0;
    RangeJob rj;

    MemCopy(&rj, r, sizeof(RangeJob));
    for (i = r->lo; i < r->hi; i += CORE_SUB_RANGE)
    {
        rj.lo = i;
        rj.hi = i + CORE_SUB_RANGE;
        if (rj.hi > r->hi)
            rj.hi = r->hi;
        res += TestCoreSubRange(&rj);

        n = rj.hi - rj.lo;
        lock {progress1 += n;}

        Yield;
    }

    return res;
}

I64 MagicPairs(I64 N)
{
    F64          t0 = tS;
    I64          res = 0;
    I64          sqrt_primes, cbrt_primes, num_powprimes, 
                 i, k, n = (N - 1) / mp_count + 1;
    Prime       *primes = PrimesNew(N, &sqrt_primes, &cbrt_primes);
    PowPrime    *powprimes = PowPrimesNew(N, sqrt_primes, primes, &num_powprimes);

    "N:%u SqrtPrimes:%u CbrtPrimes:%u PowersOfPrimes:%u\n", N, sqrt_primes, cbrt_primes, num_powprimes;
    progress1 = 0;
    *progress1_desc = 0;
    progress1_max = N;
    k = 2;
    for (i = 0; i < mp_count; i++)
    {
        rj[i].doc           = DocPut;
        rj[i].num           = i;
        rj[i].lo            = k;
        k += n;
        if (k > N)
            k = N;
        rj[i].hi            = k;
        rj[i].N             = N;
        rj[i].sqrt_primes   = sqrt_primes;
        rj[i].cbrt_primes   = cbrt_primes;
        rj[i].primes        = primes;
        rj[i].cmd           = JobQueue(&TestCoreRange, &rj[i], mp_count - 1 - i, 0);
    }
    for (i = 0; i < mp_count; i++)
        res += JobResGet(rj[i].cmd);
    Free(powprimes);
    Free(primes);
    "Found:%u Time:%9.4f\n", res, tS - t0;
    progress1 = progress1_max = 0;

    return res;
}

MagicPairs(1000000);