/*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,result=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); result++; } } } Free(pow_cnt); Free(residue); Free(sumfacts); return result; } #define CORE_SUB_RANGE 0x1000 I64 TestCoreRange(RangeJob *r) { I64 i,result=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; result+=TestCoreSubRange(&rj); } return result; } I64 MagicPairs(I64 N) { F64 t0=tS; I64 result=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++) result+=JobResultGet(rj[i].cmd); Free(powprimes); Free(primes); "Found:%u Time:%9.4f\n",result,tS-t0; return result; } 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 */