|
PROGRAM
|
Program from David Wasserman (dwasserm(AT)earthlink.net), Sep 19 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;
{
a(n) =
local(r, num2, num5, d, M, pLeft, mainP, searchP, fixed, x, rNeeded, y, nextP);
r = n;
d = digitcount(n);
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 - num2 - num5;
mainP = if (num2 == d, 2, 3);
searchP = min(4, pLeft);
fixed = 2^num2*5^num5;
x = mainP^(pLeft - searchP);
rNeeded = lift(Mod(r, M)/Mod(x, M));
while (bigomega(rNeeded) != searchP,
rNeeded += M
);
y = fixed*x*rNeeded;
if (mainP == 2,
nextP = 3,
nextP = if (num5 == d, 5, 7)
);
while (searchP < pLeft && fixed*x*nextP^(1 + searchP)/mainP < y,
searchP++;
x /= mainP
);
rNeeded = lift(Mod(r, M)/Mod(x, M));
while (bigomega(rNeeded) != searchP,
rNeeded += M
);
return(fixed*x*rNeeded);
} (End)
|
