| | |
Triangles and Squares Revisited...
Recall a number of dots, $T_n$, of the form $\displaystyle{T_n = 1 + 2 + 3 + \cdots + n = \frac{n(n+1)}{2}}$ where $n$ is a positive integer, can be arranged in a triangular form with $n$ dots on each side.
Recall also, that some numbers of dots, like $36$, can form both a triangle and a square.
Suppose you wish to find all such numbers from 1 to $T_{1500000}$ using R. (Notably -- this is a range of values significantly longer than any row or column allowed in Excel!)
To this end, write the following three functions:
T(n), which returns $T_n$ for the given value of $n$.
is.perfect.square(x) which returns $TRUE$ if $x$ is a perfect square and $FALSE$ otherwise.
Hint: you may find it useful to use the R function trunc(x), which removes the decimal part of any real value, as squaring the truncated square root of $x$ will only equal $x$ if it is a perfect square.
nums.that.can.form.a.triangle.and.square(n), which will return all values of $T_k$ where $1 \le k \le n$ and $T_k$ is a perfect square.
What is the sum of the (only eight) $T_k$ values that are also perfect squares, where $k$ ranges from 1 to $1500000$?
is.perfect.square = function(x) {
return(trunc(sqrt(x))^2 == x)
}
T = function(n) {
return(n*(n+1)/2)
}
nums.that.can.form.a.triangle.and.square = function(n) {
triangular.nums = T(1:n)
return(triangular.nums[sapply(triangular.nums,is.perfect.square)])
}
sum(nums.that.can.form.a.triangle.and.square(1500000))
Chebyshev's Theorem puts limits on the proportion of a data set that can be beyond some number of standard deviations from its mean.
Recall, the mean and standard deviations associated with a sample of size $n$ are given by the following, respectively. $$\overline{x} = \frac{\sum x}{n} \quad \quad \textrm{and} \quad \quad s = \sqrt{\frac{\sum(x-\overline{x})^2}{n-1}}$$ where for each sum above, $x$ ranges over all data values in the sample.Also note that the built-in R functions mean(v) and sd(v) calculate these two sample statistics, $\overline{x}$ and $s$.
Use these two functions, to create a function proportion.within.k.std.dev.of.mean(sample.data,k), that finds the proportion of values in the vector sample.data that lie strictly between $(\overline{x} - k \cdot s)$ and $(\overline{x} + k \cdot s)$.
Store the following data in a vector named data, and then use this function to verify the proportion of the following data that falls within 1.1 standard deviations of the mean is 0.64.
16, 22, 31, 31, 35, 55, 72, 45, 11, 4, 70, 41, 71, 88, 21, 5, 86, 23, 91, 19, 63, 70, 17, 12, 49, 82, 40, 56, 25, 40, 37, 5, 83, 26, 16, 24, 3, 28, 10, 20, 61, 41, 10, 48, 61, 8, 69, 48, 26, 13, 82, 77, 11, 92, 88, 65, 40, 75, 88, 15, 34, 81, 15, 43, 69, 4, 58, 70, 50, 40, 77, 63, 42, 95, 52, 60, 79, 79, 73, 93, 18, 36, 83, 22, 53, 39, 10, 87, 78, 37, 34, 65, 100, 90, 100, 80, 28, 5, 59, 23Note that the difference between the actual proportion found ($0.64$) and the bound established by Chebyshev's Theorem ($1 - 1/1.1^2 \doteq 0.1736$) is surprisingly large.
One might wonder how much lower one can make the actual proportion (and thus get closer to Chebyshev's bound), if you limit yourself to adding only a single value to the data set.
To this end, you run the following code:
proportion.with.extra.value = function(extra.val) {
proportion.within.k.std.dev.of.mean(c(data,extra.val),1.1)
}
min(sapply(1:100,proportion.with.extra.value))
Rounded to the nearest thousandth, what minimum proportion results?
proportion.within.k.std.dev.of.mean = function(sample.data,k) {
s = sd(sample.data)
x.bar = mean(sample.data)
return(sum((sample.data > x.bar - k*s) &
(sample.data < x.bar + k*s))/length(sample.data))
}
data = ...
proportion.within.k.std.dev.of.mean(data,1.1)
proportion.with.extra.value = function(extra.val) {
proportion.within.k.std.dev.of.mean(c(data,extra.val),1.1)
}
min(sapply(1:100,proportion.with.extra.value))
Sometimes one can find it useful to use a piecewise-linear function to approximate some other function.
Above we see the function shown in red, and given by
$$f(x) = 0.9 - \frac{18}{5} (x- 0.5)^2$$
can be fairly well approximated by the piecewise-linear function shown in black, and given by
$$f_{approx}(x) = \left\{
\begin{array}{cl}
3x & \textrm{ if } \quad x \lt 0.2\\
0.4 + x & \textrm{ if } \quad 0.2 \le x \lt 0.5\\
1.4 - x & \textrm{ if } \quad 0.5 \le x \lt 0.8\\
3 - 3x & \textrm{ otherwise }
\end{array}
\right.$$
Write two functions f(x) and f.approx(x) that evaluate $f(x)$ and $f_{approx}(x)$ as defined above, respectively.
Upon doing this correctly, the following R code should produce the same plot as shown at the beginning of this problem:
xs = seq(from=0,to=1,by=0.01) ys = sapply(xs,f.approx) plot(xs,ys,type="l",xlab="",ylab="") y2s = sapply(xs,f) lines(xs,y2s,col="red")One should be careful, however, about combining multiple applications of the approximating function -- the little errors that result from each application can add up.
To see this, use R to find the sum, to the nearest thousandth, of the differences $(f(x)-f_{approx}(x))$ as $x$ ranges over the values $0, 0.01, 0.02, 0.03, \ldots, 0.99, 1$.
As a check of your code, the same sum when $x$ ranges over the values $0, 0.1, 0.2, \ldots, 0.9, 1$ should be $0.24$.
f.approx = function(x) {
if (x < 0.2) {
return(3*x)
} else if ((0.2 <= x) && (x < 0.5)) {
return(0.4+x)
} else if ((0.5 <= x) && (x < 0.8)) {
return(1.4 - x)
} else {
return(3-3*x)
}
}
f = function(x) {
.9-(18/5)*(x-.5)^2
}
xs = seq(from=0,to=1,by=0.01)
sum(sapply(xs,f)-sapply(xs,f.approx))
Consider the data set of integers given here. In an effort to visualize where the outliers fall for a given data set, create a function outliers.of(data,lump.num) that produces a histogram of the elements in the vector named data in the following way:
The breaks used in the histogram should be constructed so that lump.num consecutive integers are "lumped together" in each rectangle of the histogram.
For example, if lump.num is 1 and the minimum value in the data set is 7, then the breaks between the rectangles should occur at $\{6.5, 7.5, 8.5, \ldots\}$, ending with max(data)+0.5. Thus, 7 is the center of the first rectangle, 8 the center of the second, 9 for the third, and so on.
If instead the lump.num is 2 and and the minimum is again 7, the breaks should occur at $\{6.5, 8.5, 10.5, \ldots\}$, again ending with max(data)+1.5. In this case, both 7 and 8 are in the first rectangle, while 9 and 10 are in the second, and so on.
Beyond getting lump.num consecutive integers associated with each rectangle, note that the "lumped together" $x$ values being counted will be centered in the rectangle to which they correspond.
Note that, as a bit of an exception, the last value should always be chosen to be lump.num-0.5 past the maximum value in the data set. This is so R doesn't complain that there are data values outside the range breaks specifies.
One can pass a col argument that consists of a vector of color names (as strings of text) to many functions that plot things (e.g., boxplot(), hist(), plot(), etc.). An example is shown below using the boxplot() function:
binom.region = function(a,b,n,p) {
xs = 0:n
probs = dbinom(xs,size=n,prob=p)
cols = ifelse((xs >= a) & (xs <= b),"green","gray")
barplot(names.arg=xs,height=probs,space=0,col=cols)
}
binom.region(4,7,10,0.5)
which results in the following plot:
col argument works similarly for the hist() function, coloring consecutive bars with the colors the col vector specifies.
Use this argument to color the bars "red" if any of the numbers to which they correspond would be outlier values by the IQR test, and "green" otherwise.
The function should return a vector of all outliers values (by the IQR test) present in the data set.
If this function is used on this example data set with lump.num=2, the result should look like this:
outliers.of(example.data,lump.num=2) [1] 92 94 23 24 25 25 25 22 20 25 20 26 22 98 [15] 22 25 25 97 22 93 20 20 95 25 20 20 22 96 [29] 25 99 100
Note that above, there are 10 bars that are colored red and 31 outlier values found. When one applies the function to the data set provided at the beginning of this problem, also using a lump.num=2, what is the sum of the number of bars that are red and the number of outlier values found?
outliers.of = function(data,lump.num) {
lower.bound = quantile(data,0.25)-1.5*IQR(data)
upper.bound = quantile(data,0.75)+1.5*IQR(data)
iqr.outliers = data[data < lower.bound | data > upper.bound]
bar.bounds = seq(from=min(data)-0.5,to=max(data)+lump.num-0.5,by=lump.num)
right.bar.bounds = bar.bounds[-1]
left.bar.bounds = bar.bounds[-length(bar.bounds)]
cols = ifelse((left.bar.bounds > upper.bound |
right.bar.bounds < lower.bound),"red","green")
h = hist(data,breaks=bar.bounds,col=cols)
return(iqr.outliers)
}
data = ...
outliers.of(data,lump.num=2)
length(outliers.of(data,lump.num=2))
In computer science, a hash function is a function that turns its input into a number $h$ from $0$ to $(m-1)$, so that the input can be stored at position $h$ in some array of $m$ positions. For a hash function to work efficiently, the outputs for randomly chosen inputs should be as uniformly distributed between $0$ and $(m-1)$ as possible.
The following R function is equivalent to what the Java language uses as the default hash function when the input is a string of text.
hash = function(s,m) {
a.vals = utf8ToInt(s)
h = 0
for (i in 1:length(a.vals)) {
h = (31*h + a.vals[i]) %% m
}
return(h)
}
Test the claim that this function produces values that are uniformly distributed between $0$ and $9$ when $m=10$ by applying it to each word found here, and then conducting an appropriate hypothesis test.
To the nearest thousandth, what is the $p$-value for the test?
hash = function(s,m) {
a.vals = utf8ToInt(s)
h = 0
for (i in 1:length(a.vals)) {
h = (31*h + a.vals[i]) %% m
}
return(h)
}
url = "http://math.oxford.emory.edu/site/math117/labFactorsTablesListsAndDataFrames/words.txt"
words = read.table(url)
hash.vals = as.vector(table(factor(sapply(as.vector(words[[1]]),hash,m=10))))
chisq.test(hash.vals,p=rep(0.1,times=10))
# Ans: 0.267 (rounded to nearest thousandth)
It can be difficult in R to view and make sense of a table with a very large number of rows and columns. Often, one is interested in identifying the cells of the table with the largest values. This is especially true when these represent frequency counts of certain combinations of categorical variable values.
Write a function dominant.cells(t,k) that reports, in decreasing order, the top $k$ cell values in a two-dimensional table $t$, along with their respective factor combinations, in a manner similar to that shown below:
> # First, let's create a new (random) table: > categoryA = factor(sample(1:4,size=1000,replace=TRUE)) > categoryB = factor(sample(5:8,size=1000,replace=TRUE)) > t=table(catA=categoryA,catB=categoryB) > t catB catA 5 6 7 8 1 55 67 61 54 2 66 50 64 81 3 48 65 71 64 4 58 58 65 73 > # Now, we apply our dominant.cells() function: > dc = dominant.cells(t,5) > dc catA catB Freq 14 2 8 81 16 4 8 73 11 3 7 71 5 1 6 67 2 2 5 66
Towards this end, you may find the order(v) function useful, which creates a permutation of vector $v$ that would (through subsetting) put its elements in order. An example follows:
> v = sample(1:10) > v [1] 4 3 8 2 9 1 7 10 6 5 > order(v) [1] 6 4 2 1 10 9 7 3 5 8 > v[order(v)] [1] 1 2 3 4 5 6 7 8 9 10
Import the data set dominant.txt, and convert it to a table with the table() function. Then use your function to find the cells of this table with the top 5 counts.
Based on the result, what is the sum of both the top 5 frequencies seen in the table and the row and column pairs to which they belong?
(Note: the corresponding sum for the top 5 frequencies in the $4 \times 4$ table generated and shown above is 404.)
url =
"http://math.oxford.emory.edu/site/math117/labFactorsTablesListsAndDataFrames/dominant.txt"
dominant.df = read.table(url,sep="")
t = table(dominant.df)
dominant.cells = function(tbl,k) {
tbl.df = data.frame(tbl)
freq.order = order(tbl.df$Freq,decreasing = TRUE)
reordered.tbl = tbl.df[freq.order,]
return(reordered.tbl[1:k,])
}
dominant.cells(t,7)
# Ans: 567
In doing so, make the assumption that during the peak weeks of the shower, the times meteors will be visible to your telescope are uniformly distributed.
Note, as the time differences between successive sightings of meteors you have seen in the past are not distributed in a manner that you recognize (i.e., they are not normally distributed, nor do they follow a $\chi^2$-distribution, or $t$-distribution, etc.), you decide to simulate a bunch of time differences by generating $n$ randomly chosen, but uniformly distributed times that on average are $k$ minutes apart, and then calculating the differences desired.
Armed with these simulated time differences, one can estimate the probability that one must wait $x$ minutes (or longer) to see a meteor. Write an R function p.value.estimate(x,k,n) that simulates the appropriate differences and returns this probability.
Rounded to two decimal places, what is the $p$-value estimate produced by this function for the probability of waiting 45 minutes or longer between meteor sightings when the true difference between meteor sightings is 12 minutes, based on 100000 simulated differences?
Hint: There are multiple approaches you can take, but a particularly expeditious one makes use of the cut() and table() functions in R!
p.value.estimate = function(x,rate,n) {
times = sort(rate*n*runif(n))
times.between = times[-1] - times[-length(times)]
counts = table(cut(times.between,breaks=c(0,x,max(times.between)),labels=FALSE))
return(counts[2]/n)
}
p.value.estimate(45,12,100000)
1-pexp(45,1/12) # <-- (for confirmation) this is the true p.value
# Ans: 0.02 (rounded to two decimal places)
A concordance is an alphabetical list of words present in a text, that shows where in the text these words occur.
Complete the body of the function concordance(url) shown below, so that it returns a list whose component names (i.e., tags) are the words encountered in the source text (specified by a url - for example, "http://math.oxford.emory.edu/myfile.txt"), and whose values are vectors whose elements are the positions of the words, as given by the numbers of words into the source text where the word in question is encountered.
As an example, if the source text started with "to be or not to be...", the concordance list would have "to" as one of its component tags, with the associated vector starting with positions 1, 5, ...
The tags in the list should be appear in order by their frequency in the source material, so that the beginning of the list shows the most frequently occurring words.
concordance = function(url) {
txt = scan(url,"")
# As a result of the above command, txt[i] now gives the ith word in the source text.
# add your code here...
}
As an example, below shows the first 3 components of the list output by this function when applied to the text of the Gettysburg Address:
> concordance("gettysburg.txt")
Read 271 items
$that
[1] 25 42 63 73 87 88 98 206 216 228 233 242 255
$the
[1] 23 122 143 168 176 200 222 258 261 264 270
$we
[1] 32 55 65 99 108 112 116 152 211 229
...
Once written, use this function to construct the corresponding concordance for Edgar Allen Poe's The Raven.
What is the sum of the elements in the vectors that correspond to the top 21 most frequent words in this poem?
Note: For the Gettysburg Address concordance, the sum of the elements in the vectors that correspond to the top 17 most frequent words would be 15337.
concordance = function(url) {
txt = scan(url,"")
word.list = list()
for (i in 1:length(txt)) {
word = txt[i]
word.list[[word]] = c(word.list[[word]],i)
}
freqs = sapply(word.list,length)
return(word.list[order(freqs,decreasing=TRUE)])
}
url = paste0(c("http://math.oxford.emory.edu/site/math117/",
"labFactorsTablesListsAndDataFrames/raven.txt"))
c = concordance(url)
sum(sapply(c[1:21],sum))
# Ans: 186013
Benford's Law can often be applied to data that spans several orders of magnitude -- things as diverse as populations of countries, loan amounts approved by a bank, river drainage rates, etc. Shockingly, it tells us that in these situations the leading digit of any particular element has a much larger chance of being a 1 than any other digit -- somewhere around 30% of the time! Indeed, Benford's law suggests that 1 will be the most common leading digit, followed by 2 (around 17.6% of the time), and 3 (12.5%), and so on -- with leading digits of 9 only happening around 4.8% of the time.
Check out this video to learn more:
You decide you want to test the claims made in the video above, so you write a function first.digit(n) that returns the leading digit of $n$, and then -- knowing that you need a data set that spans several orders of magnitude -- you apply it to all powers of $2$ from $2$ to $2^{500}$
As a hint for how to find the first digit of a number $n$, note that $\log_{10} n$ rounded down gives one less than the number of digits of $n$. (By the way, the R function trunc() will round down positive integers.).
Use your function to find the probability that a randomly chosen power of $2$ from $2$ to $2^{500}$ starts with a $1$. Then, find the probability that such a number will start with a $9$. What is the sum of these two probabilities (to the nearest hundredth)?
first.digit = function(x) {
num.digits = trunc(log(x,base=10))+1
return(trunc(x/(10^(num.digits-1))))
}
powers = 2^(1:500)
first.digits = sapply(powers,first.digit)
(sum(first.digits == 1)/500 + sum(first.digits == 9)/500)
Under certain circumstances in R (and other programming languages) one can use a function being defined in the definition of that function. When a function is defined in this way, we say the definition is recursive. Allowing such definitions might initially sound absurd , but it is actually perfectly reasonable. Consider how one might define the factorial function $n!$ with the following:
$$n! = \left\{ \begin{array}{ll} 1 & \textrm{ if } n = 1\\ n \cdot (n-1)! & \textrm{ otherwise} \end{array} \right.$$
Note how $(n-1)!$ appears in the definition of $n!$. Using the above, we can find $4!$ by applying the above successively to $4!$, $3!$, $2!$, and finally $1!$, as shown below: $$\begin{array}{rcll} 4! &=& 4 \cdot 3!\\ &=& 4 \cdot (3 \cdot 2!)\\ &=& 4 \cdot (3 \cdot (2 \cdot 1!))\\ &=& 4 \cdot (3 \cdot (2 \cdot 1)))\\ &=& 4 \cdot 3 \cdot 2 \cdot 1\\ \end{array}$$ We could define the factorial function recursively in R in a similar way:
fact = function(n) {
ifelse(n==1,1,n*fact(n-1))
}
Applying this function to find $5!$, we have
fact(5) [1] 120Here's another recursive definition. This time, the function reverses the order of the elements of a vector:
reverse = function(v) {
if (length(v) == 1) {
return(v)
} else {
return(c(v[length(v)],reverse(v[1:(length(v)-1)])))
}
}
reverse(1:10)
[1] 10 9 8 7 6 5 4 3 2 1
To see how the above function works, consider the following sequences formed during its invocation:
reverse(c(1,2,3,4,5)) = c(5, reverse(c(1,2,3,4)))
= c(5, c(4, reverse(c(1,2,3))))
= c(5, c(4, c(3, reverse(c(1,2)))))
= c(5, c(4, c(3, c(2, reverse(c(1))))))
= c(5, c(4, c(3, c(2, 1))))
= c(5, c(4, c(3,2,1)))
= c(5, c(4,3,2,1))
= c(5,4,3,2,1)
As shown above, the key to writing a recursive function has two parts:
First, one has to deal with the smallest case that could happen. In the case of the factorial function, the smallest number we decided to worry about was $1$†. In the case of reverse(), the smallest case was a list of only one element.
In both of these "smallest cases", the answer was immediate, $1! = 1$ and the reverse of a list of only one element is itself.
Second, after dealing with the smallest case, one defines the function for some input in terms of the same function applied to some "smaller" case.
In the case of the factorial function, we defined $f(n)$ in terms of $f(n-1)$, and clearly $(n-1)$ is smaller in magnitude than $n$ for the inputs our function considers.
In the case of reversing a list, we defined reverse(n) in terms of the same function applied to a vector with one less element (i.e., v[1:(length(v)-1)]).
With all of this in mind, suppose one is interested in automating the process of finding out how often certain sums occur when $n$ dice are rolled.
For example, if only $n=1$ die is rolled, then there are six equally likely possibilities: $$1,2,3,4,5,\textrm{ and } 6$$ If two dice are rolled, one often organizes all of the equally likely possibilities in a table: $$\begin{array}{c|c|c|c|c|c|c|} & 1 & 2 & 3 & 4 & 5 & 6\\\hline 1 & 2 & 3 & 4 & 5 & 6 & 7\\\hline 2 & 3 & 4 & 5 & 6 & 7 & 8\\\hline 3 & 4 & 5 & 6 & 7 & 8 & 9\\\hline 4 & 5 & 6 & 7 & 8 & 9 & 10\\\hline 5 & 6 & 7 & 8 & 9 & 10 & 11\\\hline 6 & 7 & 8 & 9 & 10 & 11 & 12\\\hline \end{array}$$ If more than two dice are rolled, however, organizing things into a table form becomes more difficult.
Imagine a function that would produce the same values, but in a sorted vector of equally likely rolls -- as the following suggests:
sums.of.n.dice(1) [1] 1 2 3 4 5 6 sums.of.n.dice(2) [1] 2 3 3 4 4 4 5 5 5 5 6 6 6 6 6 7 7 7 7 7 7 8 [23] 8 8 8 8 9 9 9 9 10 10 10 11 11 12 sums.of.n.dice(3) [1] 3 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 7 7 [23] 7 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 [45] 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 [67] 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 10 10 10 10 10 10 10 [89] 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 11 11 [111] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 [133] 11 11 11 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 [155] 12 12 12 12 12 12 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 [177] 13 13 13 13 13 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 15 15 [199] 15 15 15 15 15 15 15 15 16 16 16 16 16 16 17 17 17 18Use a recursive definition in R to define such a function and then use it to find out how many ways one can roll a 19 with 5 standard dice.
†: One could make the argument that we should have started with $n=0$ as our smallest case, since $0!=1$. This is perfectly reasonable, and modifying the function to accomplish this is easily done. As a side note related to this, know that there is a way to define the factorial function for non-integers as well! Doing so involves something called the Gamma function, $\Gamma(x)$. Compare gamma(1:10) to factorial(1:10) in R to see the relationship between these two functions, and then evaluate gamma(sqrt(2)) to see the gamma function applied to a non-integer. If you are curious and wish to learn more, do an internet search for "gamma function", and see what you find!
sums.of.n.dice = function(n) {
if (n == 1) {
return(1:6)
} else {
return(sort(rep(sums.of.n.dice(n-1),each=6)+(1:6)))
}
}
rolls = sums.of.n.dice(5)
sum(rolls == 19)
A palindrome is a word or phrase that reads the same forwards or backwards, such as "madam" or "nurses run". Similarly, A palindromic number is a number that is the same when written forwards or backwards. The first few palindromic numbers are $$0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 22, 33, 44, 55, 66, 77, 88, 99, 101, 111, 121, \ldots$$ Suppose you are curious about the proportion of palindromic numbers from 1 to some given $n$.
You suspect the following function, that reverse the order of elements in a vector might be helpful:
reverse = function(v) {
if (length(v) == 1) {
return(v)
} else {
return(c(v[length(v)],reverse(v[1:(length(v)-1)])))
}
}
reverse(1:10)
[1] 10 9 8 7 6 5 4 3 2 1
With this in mind, write the following two functions in R to aid you in deciding if a given value is a palindromic number:
digits.of(n), which returns a vector with the digits of $n$ as its elements. The following provides an example of its application:
digits.of(56789) [1] 5 6 7 8 9
is.palindromic.num(n), which returns either $TRUE$ or $FALSE$ depending on whether $n$ is a palindromic number or not, respectively, consistent with the below examples:
is.palindromic.num(122787221) [1] TRUE is.palindromic.num(122687221) [1] FALSE
reverse = function(v) {
if (length(v) == 1) {
return(v)
} else {
return(c(v[length(v)],reverse(v[1:(length(v)-1)])))
}
}
digits.of = function(n) {
num.digits = ifelse(n==0,0,trunc(log(n,base=10))+1)
ks = num.digits:1
digits = (n %% (10^ks)) %/% (10^(ks-1))
return(digits)
}
is.palindromic.num = function(n) {
return(all(digits.of(n) == reverse(digits.of(n))))
}
sum(sapply(1:54321,is.palindromic.num))
At a magic show, a magician hands two sealed envelopes to a man and a woman in the audience, respectively. He then asks a child in the audience to pick some random value between 1 and 999, and shout it out to the audience. The child randomly chooses the number 547. (To make things more exciting, try the trick with your own chosen number.) After the audience hears the child's number, the magician instructs the man to open up his envelope and read aloud and follow the instructions inside. He does so, saying the following:
Write the digits of the chosen number so that they decrease from left to right.
Find the number whose digits are the reverse of the number formed in the previous step.
Given the two numbers just found, subtract the smaller one from the larger one.
Find the number whose digits are the reverse of the difference just found.
Add the last two numbers together (i.e., the difference and the "reversed" difference) and then shout this number to the audience.
Following the instructions, the man yells out "1089". The magician then directs the woman to open her envelope and reveal what is written inside. Of course, to the amazement of the crowd, the number "1089" is written on the sole piece of paper inside her envelope. (What happened with your chosen number?)
Wondering if every possible number the child could have chosen results in the same $1089$, you decide to replicate thngs in R...
Write the following functions:
digits.of(n), which returns a vector with the digits of the integer argument $n$ as its elements. Importantly, if $n=0$, the function should return a one-element vector equivalent to c(0). The following provides an example of its application:
digits.of(56789) [1] 5 6 7 8 9 digits.of(0) [1] 0
num.from.digits(digits), which does the reverse of digits.of(). Given a vector whose elements are single digits, this function returns the value they would form if placed side by side in the same order they appear in the vector. An example of this function's application is shown below
num.from.digits(c(5,6,7,8,9)) [1] 56789
trick(n), which applies the directions the man read to the crowd so that it can return the last number he calculated and shouted to the audience. In the story, the initial value picked by the audience member was 547. As such, it should be the case that
trick(547) [1] 1089Just in case you have any doubts, note that the final value is indeed $1089$ as predicted, as the following calculations show:
Number with digits decreasing: 754 754 - 457 = 297 297 + 792 = 1089
You apply the trick to all of the values from 1 to 999 with
results = sapply(1:999,trick)As you guessed, many of the answers are 1089 -- but surprisingly not all of them!
What is the second most common result when the trick is applied to these numbers?
Also -- and just for fun -- look at all of the values that didn't produce 1089 by running
which(results != 1089)Is there some easy-to-state restriction the magician could have insisted upon when asking the child to pick his initial number, so that this trick ALWAYS works?
digits.of = function(n) {
num.digits = ifelse(n==0,0,trunc(log(n,base=10))+1)
ks = num.digits:1
digits = (n %% (10^ks)) %/% (10^(ks-1))
return(digits)
}
num.from.digits = function(digits) {
num.digits = length(digits)
summands = digits*10^((num.digits-1):0)
return(sum(summands))
}
trick = function(n) {
step1Num = num.from.digits(sort(digits.of(n),decreasing = TRUE))
step2Num = num.from.digits(reverse(digits.of(step1Num)))
step3Num = step1Num - step2Num
step4Num = num.from.digits(reverse(digits.of(step3Num)))
step5Num = step3Num + step4Num
return(step5Num)
}
results = sapply(1:999,trick)
summary(factor(results))
which(results != 1089)