/*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);