|
PROGRAM
|
Contribution from David Wasserman (dwasserm(AT)earthlink.net), Sep 30 2008: (Start)
(PARI) digitcount(n, base = 10) = local(d); if (n == 0, return(1)); d = 1 + floor(log(n)/log(base)); while (n >= base^d, d++); while (n < base^(d - 1), d--); d;
leastCombo(numPrimes, firstPrime, limit, M, r) = lcHelper(numPrimes, firstPrime, limit, prod(i = firstPrime, firstPrime + numPrimes - 1, prime(i)), 1, firstPrime + numPrimes, M, r);
lcHelper(left, i, limit, minimum, found, nextP, M, r) = local(p); if (left == 0, return(if (found%M == r, found, limit))); p = prime(i); limit = lcHelper(left - 1, i + 1, limit, minimum, found*p, nextP, M, r); if (minimum*prime(nextP)/p <= limit, limit = lcHelper(left, i + 1, limit, minimum*prime(nextP)/p, found, nextP + 1, M, r)); limit;
{
\This function returns 0 if its assumptions fail; otherwise it returns the correct answer. With the defaults 3 and 40, no failures occur in the first 1000 terms.
a(n, searchP = 3, searchSize = 40) =
local(r, num2, num5, d, M, pLeft, assumedP, fixed, x, rNeeded, y, best, z, zz, j, z3, pLimit);
r = n;
d = digitcount(n);
if (n < 7,
while (omega(r) != n, r += 10);
return(r)
);
while (num2 < d && !(r%2),
num2++;
r = r/2
);
while (num5 < d && !(r%5),
num5++;
r = r/5
);
M = 10^d/2^num2/5^num5;
pLeft = n - if (num2, 1, 0) - if (num5, 1, 0);
searchP = 3;
searchSize = 40;
fixed = 2^num2*5^num5;
assumedP = pLeft - searchP;
x = 3*prod(i = 4, assumedP + 2, prime(i));
rNeeded = lift(Mod(r, M)/Mod(x, M));
best = prime(n + 40)^searchP;
forvec (v = vector(searchP, i, [assumedP + 3, assumedP + searchSize + 2]), y = prod(i = 1, searchP, prime(v[i])); if (y%M == rNeeded, best = min(best, y)), 2);
if (best == searchP ^prime(n + 40), return(0));
y = fixed*x*best;
z = fixed*x*prod(i = assumedP + 3, pLeft + 2, prime(i));
if (z >= 2*y, return(0));
rNeeded = lift(Mod(r, M)/3);
return(fixed*3*leastCombo(pLeft - 1, 4, y/fixed/3, M, rNeeded));
} (End)
|
