Files
TempleOS-Rijn/Demo/MagicPairs.CPP
T
2018-09-08 21:47:19 +02:00

734 lines
14 KiB
C++

/*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 lo,hi,N,sqrt_primes,cbrt_primes;
Prime *primes;
CSrvCmd *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,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);
}
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;
k=2;
for (i=0;i<mp_cnt;i++) {
rj[i].doc=DocPut;
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],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;
return res;
}
MagicPairs(1000000);
/*N:10000000000 SqrtPrimes:9592 CbrtPrimes:325 PowersOfPrimes:19676
48:75
140:195
1050:1925
1575:1648
2024:2295
5775:6128
8892:16587
9504:20735
62744:75495
186615:206504
196664:219975
199760:309135
266000:507759
312620:549219
526575:544784
573560:817479
587460:1057595
1000824:1902215
1081184:1331967
1139144:1159095
1140020:1763019
1173704:1341495
1208504:1348935
1233056:1524831
1236536:1459143
1279950:2576945
1921185:2226014
2036420:2681019
2102750:2142945
2140215:2421704
2171240:3220119
2198504:3123735
2312024:3010215
2580864:5644415
2958500:3676491
4012184:4282215
4311024:7890575
5088650:6446325
5416280:7509159
6081680:9345903
6618080:12251679
7460004:10925915
7875450:16381925
2510904284:2819077155
8713880:13693959
8829792:18845855
1261742180:1532917659
3762043760:5406226575
9247095:10106504
12146750:16247745
12500865:12900734
13922100:31213899
14371104:28206815
22013334:37291625
22559060:26502315
23379224:26525415
23939685:31356314
1283921184:3105475295
26409320:41950359
27735704:27862695
6286924952:8651635047
28219664:32014575
2540165144:2920453095
6292254584:9645435015
33299000:58354119
34093304:43321095
37324584:80870615
40818855:42125144
41137620:84854315
5053883744:7240490655
6309307935:6333431264
1312252784:1516137615
49217084:52389315
8816125724:8954833635
5066210864:6535546575
52026920:85141719
52601360:97389039
5070571352:8580974247
3824655464:5223152535
61423340:88567059
62252000:93423519
3831852284:3991840515
64045904:70112175
66086504:69090615
66275384:87689415
68337324:141649235
72917000:115780599
1343056365:1450207250
2595146084:2641212315
5095542788:6361788411
2596144160:3585589215
76011992:87802407
77723360:145810719
1351707560:2139198039
3856140890:4239746469
6362067464:6445782135
89446860:197845235
93993830:99735705
94713300:240536075
94970204:96751395
97797104:114332175
7626796695:7940755304
5125032584:6113931255
100256480:174080799
101366342:105993657
101589320:168669879
1378574624:1996209375
102019644:162156995
106004024:129081735
106524264:231472535
5138808584:7179757815
114859884:164799635
116329136:183651663
116481464:128825415
116540304:204064175
117106064:133592175
117869940:261412235
2652982304:3552265695
1406830964:1862919435
128707425:205556894
130016744:166335255
130304888:184191111
2670130904:3556889895
134064584:140687415
139829624:165507975
142287992:193773447
147041504:150157215
152340944:233348655
153582624:374770655
154043505:210391694
154885920:421077215
157470404:199563195
5200486880:7236136479
157906430:161733825
158948000:237316959
160797560:335622279
3955132335:4184346704
2705286500:4823747739
161888544:322435295
5206687304:6978097335
164340224:189208575
165416048:214559631
168102044:171553635
170262560:281885919
2716915824:6217036175
1467581780:2313337515
172622505:175742294
6470954595:7115272604
173431335:197570264
1473700515:1506753884
178415048:226351671
2731700565:2784075434
3985532264:5323549335
1488607820:2343404979
191800455:198727544
195318500:269436699
3995575155:4415176844
2746338504:5962781495
198978584:248611815
201006414:311172785
2758505540:3947516859
207351704:237435495
208548624:350145775
210533904:516860015
1517220584:1654350615
215676320:308322399
2772801024:6327806975
1526672420:3045615579
226964240:313001199
1531953080:2411826759
227323250:316305549
227437280:396897759
1532741210:2367009189
230092730:271842885
237110120:376157079
1547906100:3705257675
242854884:449036315
2806455200:4618069599
1559119688:1783523511
1560890925:1704525074
252734468:268315515
260310375:290057624
261335025:338626574
262754384:265664175
7842803300:9309598875
2837602520:3735363879
4092444884:4677638955
274313564:350646435
4095184040:5537724759
4106225144:5370403335
285695190:559939625
4109899464:8797016375
6614603415:6669707624
293057960:419396439
5367444368:8311680111
296092784:348865935
300215895:308040104
4126705730:4364776125
4128822368:5786510751
5380195604:6754511595
4133786964:6235871915
5387747840:7749536319
313472880:953056175
6645302000:9441168399
5396517945:7144614854
4165376480:5797856799
2921793224:3770681655
341660220:794589635
1679061285:1727172314
344163140:400950459
1682066420:1903187979
353587304:420556695
354818060:390655395
2957731824:5398876175
367762395:436621604
379137590:443390409
380747960:508826439
1741044500:2573840619
4253975684:4333216635
401836344:933562055
407482460:659384739
3014994744:6565037255
4264348088:6767203911
1764267764:1858726155
417152295:442147544
5528653844:6755559915
1778869064:2144382135
3031422944:4121667615
428870144:478571775
3039432200:5445738999
1795408875:1816402964
447538784:728192415
5563884375:6135366824
460824848:638390511
5585101664:6246237855
3085130664:8057013335
1840746500:2258708859
8106645015:8257480424
479527244:530004915
5606173364:6149687115
486142620:987827555
5616810584:7529122215
499295384:592455015
4372365500:5770916739
506589020:816390819
4382760704:6802007295
1883876000:2776954719
4388870115:4587544604
3155451915:3215269364
3157070504:3958280535
1906051160:2849128359
1910285792:2105199135
3163223580:6981159395
3181556664:8096753735
548544744:1265855255
4438526780:5102023875
1939995596:2541717555
579753384:1224278615
4479663936:8868820415
4482524150:4516899849
597781820:676362435
1992599648:2732097951
599291000:1199626119
2006332664:2576707335
2012227784:2350626615
3269502495:3373298144
2020337564:2296563555
627398240:1155794079
7036928780:9857022195
3285728264:4537964535
3298072995:3495040604
645632204:654929715
2050561604:2604140475
5802929360:7904693295
653318420:794572779
7075066760:9357349239
662138295:721644104
665749784:851088615
7084112024:8644791975
2082593240:3971454759
2085372512:2327248287
7098019335:7477315064
713342504:844622295
2132355440:2898322575
724235784:1151128055
5895863792:9396728847
4643378844:8037099875
726558860:882424179
7160733944:9394335495
5905497584:7353653775
734590850:932552445
2160543104:2234106495
740532716:980389395
4666112264:6046105335
3419884580:4799782875
2169354564:4523120315
3424407504:6712800815
754727336:888514263
3442577624:3667278375
2193824940:4035324755
5950399364:7720741755
4696878944:5341833375
3450876624:6968768975
783839264:1106039775
4719190112:7016349087
787047488:951044031
4722238664:6332295735
2239676220:4668784835
817931828:856444875
820888964:921943035
3515685984:6970094495
2267171640:5701673159
4768579424:6631658655
2293334424:4236950375
6057750104:7455901095
6066178965:7843529834
4813637840:8865486639
2325991220:3373643979
4845714104:6985625415
2355309824:3119339775
2359000280:4197999399
902555024:1366942575
2361854565:2689435034
4864759760:8526383151
4878780704:6199240095
7397436500:9069574059
926979200:1782222399
2413813688:2931650631
3667561695:3712611104
2418204404:3199177995
949977644:1194885075
4921491824:4977874575
2423454264:6202721735
4923510500:5839241499
4924118024:6111546615
959131635:1010000396
4935148904:5646393495
2436393924:4013193275
2441496464:2838165615
3698827335:4129427384
6209309204:6443651115
8732225504:9818200095
2458574180:2887857819
6217560734:6263852385
987558704:1101593295
3715037920:6185489695
991087064:1131760935
3731353430:4332612585
6244641080:8560862919
2490395120:3301952895
2498946345:2589127190
1023736328:1051937271
1033002620:1526208579
1033763744:1307247711
1035813500:1370080899
1059627800:2052672999
1098115136:1298344383
1123418384:1252992015
1141074704:1648746735
1145890340:1699088859
1154706804:1982390795
1174809285:1212941114
1217353016:1489939143
1226825000:1742984919
Found:364 Time:24243.4199
*/