hschumann2/TempleOS-Source-Code

1 2/*The magic pairs problem:3 4Let SumFact(n) be the sum of factors5of n.6 7Find all n1,n2 in a range such that8 9SumFact(n1)-n1-1==n2  and10SumFact(n2)-n2-1==n111 12-----------------------------------------------------13To find SumFact(k), start with prime factorization:14 15k=(p1^n1)(p2^n2) ... (pN^nN)16 17THEN,18 19SumFact(k)=(1+p1+p1^2...p1^n1)*(1+p2+p2^2...p2^n2)*20(1+pN+pN^2...pN^nN)21 22PROOF:23 24Do a couple examples -- it's obvious:25 2648=2^4*327 28SumFact(48)=(1+2+4+8+16)*(1+3)=1+2+4+8+16+3+6+12+24+4829 3075=3*5^231 32SumFact(75)=(1+3)*(1+5+25)    =1+5+25+3+15+7533 34Corollary:35 36SumFact(k)=SumFact(p1^n1)*SumFact(p2^n2)*...*SumFact(pN^nN)37 38*/39 40//Primes are needed to sqrt(N).  Therefore, we can use U32.41class PowPrime42{43  I64 n;44  I64 sumfact; //Sumfacts for powers of primes are needed beyond sqrt(N)45};46 47class Prime48{49  U32 prime,pow_cnt;50  PowPrime *pp;51};52 53I64 *PrimesNew(I64 N,I64 *_sqrt_primes,I64 *_cbrt_primes)54{55  I64 i,j,sqrt=Ceil(Sqrt(N)),cbrt=Ceil(N`(1/3.0)),sqrt_sqrt=Ceil(Sqrt(sqrt)),56        sqrt_primes=0,cbrt_primes=0;57  U8 *s=CAlloc((sqrt+1+7)/8);58  Prime *primes,*p;59 60  for (i=2;i<=sqrt_sqrt;i++) {61    if (!Bt(s,i)) {62      j=i*2;63      while (j<=sqrt) {64        Bts(s,j);65        j+=i;66      }67    }68  }69  for (i=2;i<=sqrt;i++)70    if (!Bt(s,i)) {71      sqrt_primes++; //Count primes72      if (i<=cbrt)73        cbrt_primes++;74    }75 76  p=primes=CAlloc(sqrt_primes*sizeof(Prime));77  for (i=2;i<=sqrt;i++)78    if (!Bt(s,i)) {79      p->prime=i;80      p++;81    }82  Free(s);83 84  *_sqrt_primes=sqrt_primes;85  *_cbrt_primes=cbrt_primes;86  return primes;87}88 89PowPrime *PowPrimesNew(I64 N,I64 sqrt_primes,Prime *primes,I64 *_num_powprimes)90{91  I64 i,j,k,sf,num_powprimes=0;92  Prime *p;93  PowPrime *powprimes,*pp;94 95  p=primes;96  for (i=0;i<sqrt_primes;i++) {97    num_powprimes+=Floor(Ln(N)/Ln(p->prime));98    p++;99  }100 101  p=primes;102  pp=powprimes=MAlloc(num_powprimes*sizeof(PowPrime));103  for (i=0;i<sqrt_primes;i++) {104    p->pp=pp;105    j=p->prime;106    k=1;107    sf=1;108    while (j<N) {109      sf+=j;110      pp->n=j;111      pp->sumfact=sf;112      j*=p->prime;113      pp++;114      p->pow_cnt++;115    }116    p++;117  }118  *_num_powprimes=num_powprimes;119  return powprimes;120}121 122I64 SumFact(I64 n,I64 sqrt_primes,Prime *p)123{124  I64 i,k,sf=1;125  PowPrime *pp;126  if (n<2)127    return 1;128  for (i=0;i<sqrt_primes;i++) {129    k=0;130    while (!(n%p->prime)) {131      n/=p->prime;132      k++;133    }134    if (k) {135      pp=p->pp+(k-1);136      sf*=pp->sumfact;137      if (n==1)138        return sf;139    }140    p++;141  }142  return sf*(1+n); //Prime143}144 145Bool TestSumFact(I64 n,I64 target_sf,I64 sqrt_primes,I64 cbrt_primes,Prime *p)146{147  I64 i=0,k,b,x1,x2;148  PowPrime *pp;149  F64 disc;150  if (n<2)151    return FALSE;152  while (i++<cbrt_primes) {153    k=0;154    while (!(n%p->prime)) {155      n/=p->prime;156      k++;157    }158    if (k) {159      pp=p->pp+(k-1);160      if (ModU64(&target_sf,pp->sumfact))161        return FALSE;162      if (n==1) {163        if (target_sf==1)164          return TRUE;165        else166          return FALSE;167      }168    }169    p++;170  }171/*  At this point we have three possible cases to test1721)n==p1         ->sf==(1+p1)        ?1732)n==p1*p1      ->sf==(1+p1+p1^2)   ?1743)n==p1*p2      ->sf==(p1+1)*(p2+1) ?175 176*/177  if (1+n==target_sf) {178    while (i++<sqrt_primes) {179      k=0;180      while (!(n%p->prime)) {181        n/=p->prime;182        k++;183      }184      if (k) {185        pp=p->pp+(k-1);186        if (ModU64(&target_sf,pp->sumfact))187          return FALSE;188        if (n==1) {189          if (target_sf==1)190            return TRUE;191          else192            return FALSE;193        }194      }195      p++;196    }197    if (1+n==target_sf)198      return TRUE;199    else200      return FALSE;201  }202 203  k=Sqrt(n);204  if (k*k==n) {205    if (1+k+n==target_sf)206      return TRUE;207    else208      return FALSE;209  } else {210// n==p1*p2 -> sf==(p1+1)*(p2+1) ?  where p1!=1 && p2!=1211    // if p1==1 || p2==1, it is FALSE because we checked a single prime above.212 213    // sf==(p1+1)*(n/p1+1)214    // sf==n+p1+n/p1+1215    // sf*p1==n*p1+p1^2+n+p1216    // p1^2+(n+1-sf)*p1+n=0217    // x=(-b+/-sqrt(b^2-4ac))/2a218    // a=1219    // x=(-b+/-sqrt(b^2-4c))/2220    // b=n+1-sf;c=n221    b=n+1-target_sf;222// x=(-b+/-sqrt(b^2-4n))/2223    disc=b*b-4*n;224    if (disc<0)225      return FALSE;226    x1=(-b-Sqrt(disc))/2;227    if (x1<=1)228      return FALSE;229    x2=n/x1;230    if (x2>1 && x1*x2==n)231      return TRUE;232    else233      return FALSE;234  }235}236 237U0 PutFactors(I64 n) //For debugging238{239  I64 i,k,sqrt=Ceil(Sqrt(n));240  for (i=2;i<=sqrt;i++) {241    k=0;242    while (!(n%i)) {243      k++;244      n/=i;245    }246    if (k) {247      "%d",i;248      if (k>1)249        "^%d",k;250      '' CH_SPACE;251    }252  }253  if (n!=1)254    "%d ",n;255}256 257class RangeJob258{259  CDoc *doc;260  I64 num,lo,hi,N,sqrt_primes,cbrt_primes;261  Prime *primes;262  CJob *cmd;263} rj[mp_cnt];264 265I64 TestCoreSubRange(RangeJob *r)266{267  I64 i,j,m,n,n2,sf,res=0,range=r->hi-r->lo,268        *sumfacts=MAlloc(range*sizeof(I64)),269        *residue =MAlloc(range*sizeof(I64));270  U16 *pow_cnt =MAlloc(range*sizeof(U16));271  Prime *p=r->primes;272  PowPrime *pp;273  MemSetI64(sumfacts,1,range);274  for (n=r->lo;n<r->hi;n++)275    residue[n-r->lo]=n;276  for (j=0;j<r->sqrt_primes;j++) {277    MemSet(pow_cnt,0,range*sizeof(U16));278    m=1;279    for (i=0;i<p->pow_cnt;i++) {280      m*=p->prime;281      n=m-r->lo%m;282      while (n<range) {283        pow_cnt[n]++;284        n+=m;285      }286    }287    for (n=0;n<range;n++)288      if (i=pow_cnt[n]) {289        pp=&p->pp[i-1];290        sumfacts[n]*=pp->sumfact;291        residue [n]/=pp->n;292      }293    p++;294  }295 296  for (n=0;n<range;n++)297    if (residue[n]!=1)298      sumfacts[n]*=1+residue[n];299 300  for (n=r->lo;n<r->hi;n++) {301    sf=sumfacts[n-r->lo];302    n2=sf-n-1;303    if (n<n2<r->N) {304      if (r->lo<=n2<r->hi && sumfacts[n2-r->lo]-n2-1==n ||305            TestSumFact(n2,sf,r->sqrt_primes,r->cbrt_primes,r->primes)) {306        DocPrint(r->doc,"%u:%u\n",n,sf-n-1);307        res++;308      }309    }310  }311  Free(pow_cnt);312  Free(residue);313  Free(sumfacts);314  return res;315}316 317#define CORE_SUB_RANGE  0x1000318 319I64 TestCoreRange(RangeJob *r)320{321  I64 i,n,res=0;322  RangeJob rj;323  MemCpy(&rj,r,sizeof(RangeJob));324  for (i=r->lo;i<r->hi;i+=CORE_SUB_RANGE) {325    rj.lo=i;326    rj.hi=i+CORE_SUB_RANGE;327    if (rj.hi>r->hi)328      rj.hi=r->hi;329    res+=TestCoreSubRange(&rj);330 331    n=rj.hi-rj.lo;332    lock {progress1+=n;}333 334    Yield;335  }336  return res;337}338 339I64 MagicPairs(I64 N)340{341  F64 t0=tS;342  I64 res=0;343  I64 sqrt_primes,cbrt_primes,num_powprimes,344        i,k,n=(N-1)/mp_cnt+1;345  Prime *primes=PrimesNew(N,&sqrt_primes,&cbrt_primes);346  PowPrime *powprimes=PowPrimesNew(N,sqrt_primes,primes,&num_powprimes);347 348  "N:%u SqrtPrimes:%u CbrtPrimes:%u PowersOfPrimes:%u\n",349        N,sqrt_primes,cbrt_primes,num_powprimes;350  progress1=0;351  *progress1_desc=0;352  progress1_max=N;353  k=2;354  for (i=0;i<mp_cnt;i++) {355    rj[i].doc=DocPut;356    rj[i].num=i;357    rj[i].lo=k;358    k+=n;359    if (k>N) k=N;360    rj[i].hi=k;361    rj[i].N=N;362    rj[i].sqrt_primes=sqrt_primes;363    rj[i].cbrt_primes=cbrt_primes;364    rj[i].primes=primes;365    rj[i].cmd=JobQue(&TestCoreRange,&rj[i],mp_cnt-1-i,0);366  }367  for (i=0;i<mp_cnt;i++)368    res+=JobResGet(rj[i].cmd);369  Free(powprimes);370  Free(primes);371  "Found:%u Time:%9.4f\n",res,tS-t0;372  progress1=progress1_max=0;373  return res;374}375 376MagicPairs(1000000);377