/*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_cnt;
  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_cnt++;
    }
    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_cnt];

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_cnt =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_cnt,0,range*sizeof(U16));
    m=1;
    for (i=0; i<p->pow_cnt; i++)
    {
      m*=p->prime;
      n=m-r->lo%m;
      while (n<range)
      {
        pow_cnt[n]++;
        n+=m;
      }
    }
    for (n=0; n<range; n++)
      if (i=pow_cnt[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_cnt);
  Free(residue);
  Free(sumfacts);
  return res;
}

#define CORE_SUB_RANGE  0x1000

I64 TestCoreRange(RangeJob *r)
{
  I64 i,n,res=0;
  RangeJob rj;
  MemCpy(&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_cnt+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_cnt; 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=JobQue(&TestCoreRange,&rj[i],mp_cnt-1-i,0);
  }
  for (i=0; i<mp_cnt; 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);