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.
| 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | ||
| 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | ||
| << |
71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | >> |
| 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 | ||
| 91 | 92 | 93 | 94 | 95 | 96 | 97 | 98 | 99 | 100 |
Problem 51. This initially took over an hour and still did not finish. I had used too many loops and did not attempt to eliminate patterns. It appears the bulk of the time is taken by generating the 2556 acceptable patterns for the 6-digit primes. I modified it so that it tests each pattern family as they are generated which reduced the time to 28 seconds since the solution pattern family is positioned before 1492. Finally we removed more loops to reduce the time to 18 seconds.
prob51 <- function(primeraw)
{
# primeraw: list of primes in ascending order
# break into digits
primestr <- as.character(primeraw)
primen <- nchar(primestr)
maxprim <- nchar(primestr[length(primestr)])
for (n in 1:maxprim) {
primes <- primestr[which(primen == n)]
primex <- primeraw[which(primen == n)]
bset <- c()
for (i in 2:(2^n -1)) {
b <- expand(i, 1, F, base=2)
# for maxlen == 8, only 2 numbers can fail to be prime
# reject if remaining digit sum is divisible by 3
# because 3,6,9 will cause the number to be 0 mod 3
# reject if slot falls on the least significant position
# reject if 1 slot in pattern because remaining digit sum is 0,1,2 mod 3
# but 1,4,6 convert 2 mod 3 to 0 mod 3
# 2,5,8 convert 1 mod 3 to 1 mod 3
# reject if 2 slots in patterm because remaining digit sum is 0,1,2 mod 3
# but 11,44,77 convert 1 mod 3 to 0 mod 3
# 22,55,88 convert 2 mod 3 to 0 mod 3
# reject if 4 slots in pattern because remaining digit sum is 0,1,2 mod 3
# but 1111,4444,7777 convert 2 mod 3 to 0 mod 3,
# 2222,5555,8888 convert 1 mod 3 to 0 mod 3
# reject if 5 slots in pattern because remaining digit sum is 0,1,2 mod 3
# but 11111,44444,77777 convert 1 mod 3 to 0 mod 3,
# 22222,55555,88888 convert 2 mod 3 to 0 mod 3
# No pattern should be more than 5 as we expect to find the answer with
# 6 digits or less. And also I don't have all 7 digit primes :)
if (b[length(b)] == "0") {
if (length(which(b=="1")) == 3) {
bset <- c(bset, paste(b,collapse=""))
}
}
}
cat(n, ":", length(bset), " potential sets found.\n", sep="")
flush.console()
patterns <- c()
pslist <- sapply(primes, function(x) strsplit(x, ""))
for (bs in bset) {
b <- strsplit(bs, "")[[1]]
idx <- (n - length(b)) + which(b == "1")
zidx <- setdiff(1:n, idx)
is3 <- sapply(pslist, function(x) sum(as.numeric(x[zidx])) %% 3)
pslistFinal <- pslist[which(is3 > 0)]
patterns <- sapply(pslistFinal,
function(x) { x[idx] <- "x";paste(x,collapse="")})
patterns <- unique(patterns)
cat(" ", n, "-digit patterns generated: ", length(patterns), "\n", sep="")
flush.console()
minprime <- primex[1] -1
for (pat in patterns) {
patties <- sapply(0:9, function(i) as.numeric(gsub("x", i, pat)))
patties <- patties[which(patties > minprime)]
# binary search with findInterval()
found <- intersect(primex[findInterval(patties, primex)], patties)
if (length(found) == 8) {
cat(found,"\n")
flush.console()
return(found);
}
}
}
}
}
Problem 53. We calculate how many C(n,k) there are and subtract all C(n,k) <= 1,000,000. Also the Pascal triangle is almost symmetric. We should test k=0 to floor(n/2) and count all k twice except when k=floor(n/2) and n is even.
prob53 <- function(z, lim, y=1)
{
cnt <- 0
for (n in y:z) {
if (n > lim) {
break
}
subcnt <- 2
k <- 1
s <- n
tops <- floor(n/2)
while (k <= tops) {
if (s <= lim) {
if (k == tops && (n %% 2 == 0)) {
subcnt <- subcnt + 1
}
else {
subcnt <- subcnt + 2
}
}
else {
break
}
s <- s*(n-k)/(k+1)
k <- k + 1
}
cnt <- cnt + subcnt
}
(z+1)*(z+2)/2 - y*(y+1)/2 - cnt
}
Problem 57. This took 8.5 seconds.
prob57 <- function(k) {
d1 <- c(2)
d2 <- c(1)
cnt <- 0
for (i in 1:k) {
d <- addIt(multiplyIt(c(2), d1), d2)
n <- addIt(d, d1)
if (length(n) > length(d)) {
cnt <- cnt +1
}
d2 <- d1
d1 <- d
}
cnt
}
Problem 64. Brute force. prob64(10000) = 1322 took 11 seconds.
- \(\sqrt{x} = b_0 + \frac{\sqrt{x} - b_o}{1} =: b_0 + \frac{\sqrt{x} + n_0}{d_0}\)
- Define \(d_1 := x - n^2_0\). Then \[ \frac{d_0}{\sqrt{x} + n_0} = \frac{d_0(\sqrt{x} - n_0)}{x-n^2_0} = b_1 + \frac{d_0 \sqrt{x} - d_0n_0 - b_1d_1}{d_1} = b_1 + \frac{\sqrt{x} + (-n_1 - b_1d_1/d_0)}{d_1/d_0} \]
- Let \(n_1 := -n_0 -b_1d_1/d_0\) and \(d_1 := d_1/d_0\) so we obtain again \(b_1 + (\sqrt{x}+n_1)/d_1\).
- Check if we have encountered the pair \((n_1, d_1)\)
previously. If we have and it it is entry idx. Then the next value of b
is found in entry idx+1.
We rely on the fact (?) that \(d_1\) is always divisible by \(d_0\).
expandsqrt <- function(x, str=T)
{
# If *str* = T, will return result as a string N;P where P is the periodic portion
# else it returns a list with two attributes 'x' and 'p' where 'p' contains
# the periodic portion
kr <- c() # continued fraction
nr <- c() # numerator
dr <- c() # denominator
sqrx <- sqrt(x)
b0 <- trunc(sqrx)
d0 <- 1
n0 <- -b0
repeat {
idx <- which((nr == n0) && (dr == d0))
kr <- c(kr, b0)
nr <- c(nr, n0)
dr <- c(dr, d0)
if (length(idx) > 0) {
idx <- min(length(kr), idx + 1)
break
}
if (x == n0^2) {
# ----- continued fraction terminates
idx <- 0
break
}
d1 <- x - n0^2
b1 <- trunc(d0*(sqrx - n0)/d1)
n0 <- -b1*d1/d0 - n0
d0 <- d1/d0
b0 <- b1
}
if (str) {
if (idx == 0) {
z <- paste(kr, collapse=",")
}
else {
z <- paste(paste(kr[0:(idx-1)], collapse=","),
";",
paste(kr[idx:length(kr)], collapse=","), collapse="")
}
}
else {
if (idx == 0) {
z <- list(x=kr, p=c())
}
else {
z <- list(x=kr[0:(idx-1)], p=kr[idx:length(kr)])
}
}
z
}
prob64 <- function(n)
{
cnt <- 0
for (i in 1:n) {
p <- expandsqrt(i, F)
if (length(p$p) %% 2 == 1) {
cnt <- cnt + 1
}
}
cnt
}
Problem 65. Build continued fraction of e. Then run f <- convgfrac(e,100) which gives you a list of two strings: n (the numerator) and d (the denominator). We calculate the sum as sum(as.numeric(strsplit(f$n,"")[[1]])) = 272. This took 0.12 seconds.
e <- c(2, rep(1,120))
idx <- which(1:120 %% 3 == 0)
e[idx] <- idx*2/3
The trick is to start from the bottom and work up instead of from the top working down. So first we reverse the continued fraction resulting in R. Then we follow this recipe:
convgfrac <- function(alist, t) {
# alist: continued fraction as a vector (eg, \sqrt(3) = c(1,1,2,1,2,1,2,1,2,1,2,1,2))
# t: desired convergant.
#
# Returns f where f$n holds the numerator and f$d holds the denominator as strings.
f <- list(n=c(), d=c())
ra <- lapply(alist[t:1], function(z) expand0(z))
i <- 1
n <- c(1)
d <- ra[i][[1]]
if (t > 1) {
for (i in 2:t) {
n0 <- n
n <- d
d <- addIt(multiplyIt(d, ra[i][[1]]), n0)
}
}
f$d <- paste(n,collapse="")
f$n <- paste(d,collapse="")
f
}
Problem 84. Here we quickly ran into R 2.15's poor support for passing arguments by reference. I could pass around environments but it is a bit of a pain. Also I wanted to try OO-programming in R. I took a look at the S4 method of creating objects and knocked that out of consideration. It is so dang ugly. So I decided to use the R package R.oo which wraps around the S3 method of creating objects. This is lengthy indeed. The only remaining thing is to determine how many runs I would need to do to get a stable order of frequency. One way to do it is to run the game with 1000 turns and do this 25 times. You are then able to calculate the standard deviation of the probabilities of the top 4 squares (eg, \(\sigma_1, \sigma_2, \sigma_3, \sigma_4 \)). From that you can estimate the population standard deviation (eg, \(5\sigma_1, 5\sigma_2, 5\sigma_3, 5\sigma_4\)). Then you look at the probabilities and how much accuracy you want (\(\epsilon\)) and from that you can determine the number of turns to run (\(\ge (5\sigma/\epsilon)^2\)). I just ran it for 100,000 turns.
library(R.oo)
setConstructorS3("Prob84",
function(diceList=1:4, ccnum=16, chnum=16)
{
## board information
# board types: 1=CC, 2=CH, 3=R, 4=U, 5=G2J
#
# @note: R.oo requires all Objects to allow a default constructor
#
boardsq <- data.frame(typ=rep(0,40), name=c(""), cnt=rep(0,40))
boardsq[,2] = c('GO', 'A1','CC1','A2', 'T1','R1','B1', 'CH1','B2','B3',
'JAIL','C1','U1', 'C2', 'C3','R2','D1', 'CC2','D2','D3',
'FP', 'E1','CH2','E2', 'E3','R3','F1', 'F2' ,'U2','F3',
'G2J', 'G1','G2', 'CC3','G3','R4','CH3','H1', 'T2','H2')
boardsq[,1] = c(0, 0, 1, 0, 0, 3 , 0, 2, 0, 0,
0, 0, 4, 0, 0, 3 , 0, 1, 0, 0,
0, 0, 2, 0, 0, 3 , 0, 0, 4, 0,
5, 0, 0, 1, 0, 3 , 2, 0, 0, 0)
## community cards
#
ccList <- rep(0, ccnum)
ccList[sample(1:ccnum, 2, replace=F)] <- c(1,11)
ccListPtr <- 1
## chance cards
# 0 = do nothing
# >0 go to that square
# -3 = nearest railroad
# -4 = nearest utility
# -6 = go back 3 steps
chList <- rep(0, chnum)
chList[sample(1:chnum, 10, replace=F)] <- c(1,11,12,25,40,
6,-3,-3,-4,-6)
chListPtr <- 1
extend(Object(), "Prob84",
.boardsq = boardsq,
.chList = chList,
.ccList = ccList,
.chListPtr = chListPtr,
.ccListPtr = ccListPtr,
.doubleRolls = 0,
diceList = diceList)
})
setMethodS3("reshuffle", "Prob84",
function(this) {
## reshuffle cards
chnum <- length(this$.chList)
this$.chList <- rep(0, chnum)
this$.chList[sample(1:chnum, 10, replace=F)] <- c(1,11,12,25,40,
6,-3,-3,-4,-6)
this$.chListPtr <- 1
ccnum <- length(this$.ccList)
this$.ccList <- rep(0, ccnum)
this$.ccList[sample(1:ccnum, 2, replace=F)] <- c(1,11)
this$.ccListPtr <- 1
0
})
setMethodS3("run", "Prob84",
function(this, turns=10000) {
pos <- 1
rolls <- sample(this$diceList, 2*turns, replace=T)
for (i in 1:turns) {
pos <- nextSq(this, pos, c(rolls[2*i-1],rolls[2*i]))
this$.boardsq[pos,3] <- this$.boardsq[pos,3]+1
}
# -- print out sorted list
this$.boardsq[order(this$.boardsq[,3]),]
})
setMethodS3("nextSq", "Prob84",
function(this, n, roll) {
## calculate next position
nsq <- ((n + sum(roll) -1) %% 40) + 1
if (roll[1] == roll[2]) {
this$.doubleRolls <- this$.doubleRolls + 1
if (this$.doubleRolls == 3) {
this$.doubleRolls <- 0
return(11)
}
}
else {
this$.doubleRolls <- 0
}
ctype <- this$.boardsq[nsq,1]
if (ctype == 2) {
chCard <- this$.chList[this$.chListPtr]
this$.chListPtr <- (this$.chListPtr%%(length(this$.chList))) + 1
if (chCard == 0) {
return(nsq)
}
else if (chCard > 0) {
return(chCard)
}
else if (chCard == -3) {
return(Prob84$nearestR(nsq))
}
else if (chCard == -4) {
return(Prob84$nearestU(nsq))
}
else if (chCard == -6) {
return(advance(this, nsq, -3))
}
stop(simpleError(paste("unknown Community Card", chCard)))
}
else if (ctype == 1) {
ccCard <- this$.ccList[this$.ccListPtr]
this$.ccListPtr <- (this$.ccListPtr%%(length(this$.ccList))) + 1
if (ccCard > 0) {
return(ccCard)
}
else {
return(nsq)
}
}
else if (ctype == 5) {
return(11)
}
return(nsq)
}
)
setMethodS3("advance", "Prob84",
function(this, n, steps) {
## go forward 'steps'
nsq <- ((n + steps -1)%%(length(this$diceList))) +1
return(nsq)
})
setMethodS3("nearestR", "Prob84",
function(this, n) {
## find nearing railroad given current location
# @warn: Hardcoded
if (n == 8) {
return(16)
}
else if (n == 23) {
return(26)
}
else if (n == 37) {
return(6)
}
stop(simpleError(paste("Error nearest railroad:", n)))
},
static=TRUE,
protected=TRUE)
setMethodS3("nearestU", "Prob84",
function(this, n) {
## find nearest utility givven current location
# @warn Hardcoded
if (n == 8) {
return(13)
}
else if (n == 23) {
return(29)
}
else if (n == 37) {
return(13)
}
stop(simpleError(paste("Error nearest utility:", n)))
},
static=TRUE,
protected=TRUE)
Problem 91.
prob91 <-function(n) {
cnt <- 0
aseq <- 1:(n-1)
if (n == 1) {
aseq <- c()
}
for (a in aseq) {
for (b in 1:n) {
p <- a^2 + b^2
for (c in (a+1):n) {
for (d in 0:max(b-1,1)) {
if (p == a*c + b*d) {
cnt <- cnt + 1
}
}
}
}
}
3*n^2 + 2*cnt
}
Problem 92. Fairly straightforward. Just have to remember previous results so we don't redo too much work. But I don't want to remember everything! After calculating the first link in the chain, the maximum value of the sum is 9^2*(number of digits). So we will make space for 9^2*(number of digits). This version ran in 376 seconds and loops 19,984,731 times. When byte compiled with
compiler.cmpfun, it ran in 278 seconds.
prob92 <- function(n)
{
## Find number of chains that end at 89 from all integers less than n
#
# @return 2-element vector with the first element indicating how many loops were executed,
# and the second the number of integers that end at 89
loops <- 0
count <- 0
# For large enough n, the sum can't be larger than 9^2*(number of digits in n)
chainSize <- max(89, 81*round(1+log(n)/log(10)))
chains <- rep(NA, chainSize)
chains[1] <- 1
chains[89] <- 89
for (i in 1:(n-1)) {
v <- expand0(i)
s <- sum(v*v)
loops <- loops + 1
while (s != 1 && s != 89) {
olds <- s
v <- expand0(s)
s <- sum(v*v)
finalS <- chains[s]
if (!is.na(finalS)) {
chains[olds] <- finalS
s <- finalS
}
loops <- loops + 1
}
if (s == 89) {
count <- count + 1
}
}
c(loops, count, length(chains))
}
Surprisingly though, if we decided to remember the values for all integers, it actually runs faster: prob92b runs in 249 seconds and performs only 10,000,106 loops because we don't have to shrink the number beforehand.
prob92b <- function(n)
{
## Find number of chains that end at 89 from all integers less than n
# @return 2-element vector with the first element indicating how many loops were executed,
# and the second the number of integers that end at 89
loops <- 0;
chains <- rep(NA, n-1)
chains[1] <- 1
chains[89] <- 89
for (i in 1:(n-1)) {
s <- i
while (s != 1 && s != 89) {
v <- expand0(s)
s <- sum(v*v)
finalS <- chains[s]
if (!is.na(finalS)) {
chains[i] <- finalS
s <- finalS
}
loops <- loops + 1
}
}
c(loops, length(which(as.matrix(chains)[,1] == 89)))
}
And just for fun, we have this function for finding a chain:
sqchain <- function(n)
{
## Find chain of positive integer n
# @return chain as vector
result <- c(n)
s <- n
if (n == 89) {
v <- expand0(s)
s <- sum(v*v)
result <- c(result,s)
}
while (s != 1 && s != 89) {
v <- expand0(s)
s <- sum(v*v)
result <- c(result, s)
}
result
}
Problem 94. Heron's formula states that the area of a triangle with sides a, b, c is \[ \sqrt{s(s-a)(s-b)(s-c)} \] where the semi-perimeter \(s := (a+b+c)/2\).
We have two cases: case A: \((a, a, a+1)\) and case B: \((a, a, a-1)\). Then the area in case A is \[ \frac{a+1}{4}\sqrt{(3a +1)(a-1)} \] and the area in case B is \[ \frac{a-1}{4}\sqrt{(3a -1)(a+1)} \quad. \] . For case A, we need \((3a+1)(a-1) = x^2\) for some positive integer \(x\). Solving for a we have that \[ a = \frac{-1 + \sqrt{4 + 3x^2}}{3} \quad.\] Therefore we must solve the Diophantine equation \( 4+3x^2 = y^2 \) for some positive integer y. Being ignorant of ways to solve such an equation, we defer to someone with that knowledge: Dario Alpern's Generic Two integer variable equation solver. This gives us the recursive relationship to calculate (x,y) for values up to 1,000,000,000. The final answer sum(prob94()) took 0.02 seconds.
prob94 <- function(LIMIT=1e9) {
# Use recurrence relationship to solve Diophantine equation
x <- 0
y <- 2
ylist <- c(y)
xlist <- c(x)
while (y <= LIMIT) {
x1 <- 2*x + 1*y
y1 <- 3*x + 2*y
xlist <- c(xlist, x1)
ylist <- c(ylist, y1)
x <- x1
y <- y1
}
# Prepare
alist <- c() # values of a
plist <- c() # perimeter values
for (i in 1:length(ylist)) {
y <- ylist[i]
x <- xlist[i]
m <- y %% 3
if (m == 1) {
# this is case B. The modulo test is to ensure a is an integer
a <- (y -1)/3
if (a > 1) {
if ( ((a-1)*(x %% 4)) %% 4 == 0) {
if (3*a -1 <= LIMIT) {
alist <- c(alist, a)
plist <- c(plist, 3*a - 1)
}
}
}
}
else if (m == 2) {
# this is case A.
a <- (y +1)/3
if (a > 1) {
if ( ((a+1)*(x %% 4)) %% 4 == 0) {
if (3*a +1 <= LIMIT) {
alist <- c(alist, a)
plist <- c(plist, 3*a +1)
}
}
}
}
}
plist
}
Problem 97. Using our previously implemented functions we easily find 8739992577 in 0.04 sec.
lastDigits(addIt(c(1), multiplyIt(c(2,8,4,3,3), powerIt(2, 7830457, 10))), 10)
Problem 100. Let B be the number of blue discs and N be the total number of discs. Then
\[\begin{align} \frac{B}{N}\frac{B-1}{N-1} &= \frac{1}{2} \\ 2B(B-1) &= N(N-1) \\ 2B^2 - N^2 - 2B +N = 0 \end{align} \] B=1 and N=1 is one solution. We use Dario Alpern's Generic Two integer variable equation solver to find that the k-th solution is given by
\[\begin{align} B_{k+1} = 3B_k + 2N_k -2 \\ N_{k+1} = 4B_k + 3N_k -3 \end{align} \] prob100 <- function() {
b0 <-1
n <-1
trillion <- c(1, rep(0,12))
repeat {
b <- addIt(addIt(multiplyIt(c(3), b0), multiplyIt(c(2), n)), c(-2))
n <- addIt(addIt(multiplyIt(c(4), b0), multiplyIt(c(3), n)), c(-3))
cat(paste(b,collapse=""),paste(n,collapse=""),"\n")
if (cmpIt(n, trillion) > 0) {
break
}
b0 <- b
}
b
}
No comments:
Post a Comment