The equation images were all generated by the wonderful MathJax javascript LaTeX rendering library. Another convenient resource is Online LaTeX Equation Editor. See Common Routines for functions used in the solutions.
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | |
| 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | |
| 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | >> |
| 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | |
| 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 |
Problem 2.
prob2 <- function(z) {
f0 <- 1
f1 <- 1
s <- 0
repeat {
f <- f0 + f1;
if (f > z) {
break;
}
if (f %% 2 == 0) {
s <- s + f
}
f0 <- f1;
f1 <- f;
}
s
}
The answer is then prob2(4e6) = 4613732. Can't really see a faster way of doing this but probably doesn't matter as it was fast enough that system.time() didn't register.
Problem 9.
\[
\begin{align}
a+b+c = a+b + \sqrt{a^2 + b^2} &\stackrel{\triangledown}{=} K \\
1 + (b/a) + \sqrt{1 + (b/a)^2} &= K/a
\end{align}
\]
which we can simplify to find
\[
\begin{align}
b &= \frac{K}{2}\frac{2a -K}{a -K} \\
&= K - \frac{K^2}{2(K-a)}
\end{align}
\]
We see \(K^2/2\) must be divisible by \(K - a\). So in our case 500,000 must be divisible by 1000 - a. 500,000 = 2^5 * 5^6 and so there are only 6*7 = 42 possible values of a which we have to check. And out of those 42, we want the one where \(0 < a = K/2\) and \(a = b\). Runtime was quick enough that it didn't register with system.time().
prob9 <- function() {
k <- 1000
answ <- c(0, 0, 0);
for (i in 0:5) {
for (j in 0:6) {
f <- (2^i) * (5^j)
a <- k - f
if (a > 0 && a <= k/2) {
b <- k - k^2/2/f
if (b >= a) {
answ <- c(a, b, sqrt(a^2 + b^2));
return
}
}
}
}
answ
}
prob9() = (200, 375, 425) so abc = 31875000.
Problem 12. Brute force using a vector of primes as input. prob12(500) = 12375 which means the triangular number is 76576500 = (2^2)*(3^3)*(5^3)*7*11*13*17 with 576 divisors. This however took 108 seconds.
prob12 <- function(minDiv, primes) {
i <- 1
divs <- 0
repeat {
i <- i + 1
x <- i*(i+1)/2
divs <- numDivisors(x, primes)
if (divs > minDiv) {
break
}
}
i
}
A more direct way of calculating divisors without factoring is the following which resulted in a 19 second runtime.
numDivisors <- function(n) {
cnt <- 2
m <- trunc(sqrt(n))
if (m == 1) {
return (cnt)
}
if (n %% m == 0) {
cnt <- cnt + 1
}
d <- 2:(m-1)
cnt <- cnt + 2*length(which (n %% d == 0))
return (cnt)
}
Finally, consider x and y where gcd(x,y) = 1. Then the numDivisors(xy) = numDivisors(x)*numDivisors(y). In our case gcd(i, i+1) = 1 so gcd(i/2, i+1) = 1. numDivisors(N) runtime is \(O(\sqrt{N})\). Then \(O(\sqrt{X*Y}) > O(\sqrt{X}) + O(\sqrt{Y})\). The following version took 2.4 seconds.
prob12b <- function(minDiv) {
i <- 1
repeat {
i <- i + 1
divsA <- numDivisors(i/2)
divsB <- numDivisors(i+1)
if (divsA*divsB > minDiv) {
break
}
}
i
}
Problem 14. Another brute force method. However I remember some of the prior Collatz calculations to save time. Specifically I remember those from 1 to 2*N where N in our case is 1,000,000. I kept track of how many iterations I saved by using the table (savings) versus actually iterations (miss + savings). This function took 55 seconds to find the answer with a savings of 96.8% (that's 127,208,163 iterations).
prob14 <- function(x) {
savings <- 0
miss <- 0
maxc <- 0
maxn <- 0
maxt <- 2*x
collatzTable <- rep(0, maxt);
collatzTable[1] <- 1;
for (i in 1:x) {
len <- 1
n <- i
while (n > 1) {
if (n %% 2 == 0) {
n <- n/2
}
else {
n <- 3*n + 1
}
if (n <= maxt && collatzTable[n] > 0) {
collatzTable[i] <- collatzTable[n] + len
savings <- savings + collatzTable[n]
break
}
len = len + 1
miss <- miss + 1
}
if (collatzTable[i] > maxc) {
maxn <- i
maxc <- collatzTable[i]
}
}
cat("savings = ", savings, " (", 100*savings/(savings+miss), "%), misses = ", miss, "\n");
cat("n = ", maxn, " with chain length ", maxc, "\n");
c(maxn, maxc)
}
Problem 20. I couldn't really derive a formula for this. So I did it the brute force way by implementing infinite precision integer multiplication.
prob20 <- function(n) {
k <- 1
fact <- c(1)
if (n > 1) {
for (k in 2:n) {
kc <- sapply(strsplit(as.character(k),'')[[1]], as.numeric)
fact <- multiplyIt(fact, kc)
print(sum(fact))
}
}
sum(fact)
}
Then prob20(100) = 648. This took 0.3 seconds on my laptop.
I thought it would be interesting to see this function for up to 500! This took 14 seconds to calculate. It is slightly bigger than linear but at the same time it was definitely not N^p for p>1. My guess is it is something like O(N log (N)).
Problem 25. One could use the fact that
\[ F(n) = \frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n \right] \]
Since the second term is approximately (-0.11)^n, it rapidly approaches zero. Then for large n we only need to consider the first term. We can estimate the number of digits of F(n) as
\[ \text{digits} = 1 + \lfloor \log_{10}\frac{1}{\sqrt{5}} + n\log_{10} \frac{1+\sqrt{5}}{2} \rfloor \stackrel{\triangledown}{=} 1000 \] \[ \begin{align} \therefore n &= \left( 999 - \log_{10}\frac{1}{\sqrt{5}}\right)/\log_{10} \frac{1+\sqrt{5}}{2} \\ &\approx 4781.86 \end{align} \] This implies 4781 would not have 1000 digits; therefore n = 4782 which gives us the correct answer.
We can implement a brute force method by using infinite precision integer arithmetic. This program took 30 seconds to find that prob25(1000) = 4782.
prob25 <- function(n) {
if (n == 1) {
return (1)
}
f0 <- c(1)
f1 <- c(1)
i <- 3
repeat {
f <- addIt(f0, f1)
if (length(f) >= n) {
break
}
i <- i+1
f0 <- f1
f1 <- f
}
c(i, paste(f, collapse=""))
}
However we can write a different addIt where it breaks numbers up into blocks (of 8 digits) and uses arithmetic addition on each block before combining them. The following procedure took 4 seconds to find the first 1000 digit Fibonnaci number.
prob25b <- function(n) {
if (n == 1) {
return (1)
}
f0 <- 1
f1 <- 1
i <- 3
blocks <- ceiling(n/8)
repeat {
f <- addIt2(f0, f1, F, 8)
if (length(f) >= blocks) {
if (digitsIt2(f, 8) >= n) {
break
}
}
i <- i+1
f0 <- f1
f1 <- f
}
c(i, stringIt(f, 8))
}
Problem 26. Brute force. We implement fraction expansion. The answer is 1/983 which has an enormous 982 period. I'm guessing there is a algebraic explanation for this. But I don't know Algebra very well. The function took 4 seconds.
prob26 <- function(n) {
maxd <- 0
maxl <- 0
for (d in 1:n) {
f <- expand(1,d)
idx <- which(f == ";")
if (length(idx) > 0) {
newlen <- length(f) - idx
if (newlen > maxl) {
maxd <- d
maxl <- newlen
}
}
}
c(maxl, maxd)
}
Problem 30. We notice two things. First, no number can be larger than 6 digits because 10^6 > 6*9^5 = 354294 (ie, if all 6 digits are 9s). We found 6 numbers: 4150, 4151, 54748, 92727, 93084, and 194979. The sum totals to 443839. This took 8 seconds. We also noticed that x^5 mod 10 = x.
prob30 <- function() {
powers <- c((1:9)^5)
x <- c(0, 0, 0, 0, 0, 2)
v <- 2
maxv <- 6*9^5 + 1
while (v < maxv) {
if (sum(powers[x]) == v) {
cat(x, " = ", v, "\n")
}
carry <- 1
i <- length(x)
while (i > 0 && carry > 0) {
p <- x[i] + carry
x[i] <- p %% 10
carry <- trunc(p/10)
i <- i-1
}
v <- v + 1
}
}
Problem 35. Straightforward brute force. This took 46 seconds.
bsearch <- function(x, v) {
## Simple binary search of a integer vector
# @param x value to search
# @param v vector to search from
# @return index within v.
a <- 1
b <- length(v)
while (b >= a) {
# cat(a,"-",b,"\n")
idx <- a + (b-a+1)%/%2
if (v[idx] == x) {
return(idx)
}
else if (v[idx] > x) {
b <- idx-1
}
else {
a <- idx+1
}
}
0
}
isprimeLite <- function(x, primes) {
## This is called primeLite because it requires a predefined list
# of primes
maxprime <- primes[length(primes)]
if (x <= maxprime) {
return (bsearch(x, primes) > 0)
}
else if (sqrt(x) <= maxprime) {
return (length(which(x %% primes[1:round(sqrt(maxprime))] == 0)) == 0)
}
NA
}
rotL.int <- function(x, d) {
## rotate integer x which has d digits
f <- x%/%(10^(d-1))
(x - f*10^(d-1))*10 + f
}
prob35 <- function(n, primes) {
# d and dtrigger used to increment nbr of digits
d <- 1
dtrigger <- 10
count <- 0
for (x in 1:(n-1)) {
if (x == dtrigger) {
d <- d+1
dtrigger = 10*dtrigger
}
if (isprimeLite(x, primes)) {
cirprime <- T
if (d > 1) {
for (j in 1:(d-1)) {
x <- rotL.int(x, d)
if (!isprimeLite(x, primes)) {
cirprime <- F
break
}
}
}
if (cirprime) {
count <- count + 1
}
}
}
count
}
Problem 36. Another brute force method. It is clear we should only look at odd numbers since the binary representation should end in 1. This took 16 seconds.
prob36 <- function(k) {
s <- 0
x <- c(1)
xlen <- length(x)
for (xv in seq(1, k, by=2)) {
spt <- max(1, trunc(xlen/2))
# for example, compare 1:5 and 11:(11-5+1)
if (identical(x[1:spt],x[xlen:(xlen - spt +1)])) {
b <- expand(xv, 1, F, 2)
blen <- length(b)
bspt <- max(1, trunc(blen/2))
if (identical(b[1:bspt], b[blen:(blen - bspt +1)])) {
s <- s + xv
# cat(x, "=", b, "\n", sep="")
}
}
carry <- 2
i <- 1
while (carry > 0) {
if (i > xlen) {
x[i] <- carry
xlen <- xlen + 1
break
}
p <- x[i] + carry
x[i] <- p %% 10
carry <- trunc(p/10)
i <- i+1
}
}
s
}
Problem 40. We generate a table containing the start and end indices. We then calculate prob40(10^(0:6)) = c(1, 1, 5, 3, 7, 2, 1). The product is then 210.
generateTable <- function(n) {
mydata <- list(digits=c(), # number of digits
count=c(), # number of such numbers
idx=c(), # start index in irrational number
jdx=c()) # end index in irrational number
if (n < 1) {
return(mydata)
}
jdx <- 0
for (i in 1:n) {
mydata$digits <- c(mydata$digits, i)
cnt <- 10^i - 10^(i-1)
mydata$count <- c(mydata$count, cnt)
idx <- jdx + 1
mydata$idx <- c(mydata$idx, jdx + 1)
jdx <- idx + i*(cnt) -1
mydata$jdx <- c(mydata$jdx, jdx)
}
mydata
}
prob40 <- function(idx, n = 6)
{
mydata <- generateTable(n)
p <- sapply(idx, function(k) rev(which(mydata$idx <= k))[1])
n <- mydata$digits[p]
nbr <- 10^(n-1) + trunc((idx - mydata$idx[p])/n)
digit <- (idx - mydata$idx[p]) %% n
trunc(nbr/10^(n-digit-1)) %% 10
}
Problem 48. Brute force again. Here we implement powerIt(x, n, y) which computes x^n and retains last y digits. Prob48(1000) = 9110846700 took 16 seconds.
prob48 <- function(n, digits=10)
{
s <- c(0)
for (i in 1:n) {
s <- lastDigits(addIt(s, powerIt(i, i, digits)), digits)
}
s
}
Problem 49. Brute force. We find 2969, 6299, 9629. This took 0.5 seconds.
prob49 <- function(primes) {
idx <- which(primes > 1000)[1]
jdx <- which(primes > 10000)[1] -1
code <- data.frame(c=rep(0, jdx), # hash code
n=rep(0, jdx), # number of occurences
done=rep(0, jdx)) # done flag
# we generate a hash code for each number which is just the digits sorted
# (eg, 4909 --> 499)
#
for (i in idx:jdx) {
code$c[i] <- as.numeric(paste(sort(strsplit(as.character(primes[i]),"")[[1]]),collapse=""))
}
for (i in idx:jdx) {
code$n[i] <- length(which(code$c == code$c[i]))
}
for (i in idx:jdx) {
if (code$done[i] == 0 && code$n[i] >= 3) {
kperms <- which(code$c == code$c[i])
code$done[kperms] <- 1
kpermslen <- length(kperms)
for (bidx in 1:(kpermslen-2)) {
b <- primes[kperms[bidx]]
for (bjdx in (bidx+1):(kpermslen-1)) {
delta <- primes[kperms[bjdx]] - b
mymatch <- which (2*delta + b == primes[kperms])
if (length(mymatch) > 0) {
cat (c(b, primes[kperms[bjdx]], primes[kperms[mymatch]]), "\n")
}
}
}
}
}
}
No comments:
Post a Comment