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