Euler 49

There is a sequence 1487, 4817 and 8147 that has three properties:
1. All numbers are primes
2. They are all permutations of each other
3. Each term increases by a certain number (in this case 3330)

Project Euler claims that there are two such sequences with four-digit numbers. Our task is to find the other.
The answer is the three primes with the required properties, pasted together in ascending order.

As usual I’ll load the numbers library. And the first task must be to find all four-digit primes.

library(numbers)
cand <- Primes(1000,9999)

But we will not need these prime. We already know that answer.

cand <- cand[-which(cand %in% c(1487, 4817, 8147))]

If one of these four-digit primes is part of such a sequence, there must exist at least two other primes, which are permutations of it.
Lets begin by generating that list
The library gtools provides the function permutations()

library(gtools)

listing <- function(x){
  b <- unlist(strsplit(as.character(x),""))
  c <- permutations(4,4,v=b, set=FALSE)
  d <- apply(c, 1, paste, collapse="")
  d <- as.numeric(d)
  sort(unique(d[d %in% cand])) 
}

We take a four-digit number as input, converts it to characters via as.character(), and uses strsplit() to get a character vector with for numbers.
Then permutations. The first “4” indicates that the source vector (b) has 4 elements. The second “4” that we want results that contains 4 elements. That we take “b” as the source vector – i.e. the elements from which the permuations should be made. And set=FALSE is a lgical flag. Default is TRUE, and will remove duplicates. That is not what we want, so it is set to FALSE.
Then all the permutations are collapsed to four-digit strings, converted to numeric, and duplicates are removed with unique(), returning a list of all permuations of the digits in the input, that are in the list of candidate primes.

Thats nice. But we know that there should be three numbers in the sequence. This function simply tests if there are at least three numbers in the listing:

poss <- function(x){
  length(listing(x)) >= 3
}

That was the simple part. We can now take all the four-digit primes, excluding the three we already know about, make permutations of them, and see if there are at least three permutations among them, that are prime.

Now it gets a bit more complicated. Given a set of four-digit numbers, we need to locate a subset of three, A, B and C. The differences C-B and B-A should be equal.
There are probably more elegant ways to do that, but this is the way I found.
Given a number x, get all possible permutations of that:
listing(x)
Get all the combinations of three of these numbers:
comb(listing(x),3)
Now we have a matrix, with all the possible combinations in the columns.
Find the differences between the first and the second row, and the second and the third row, for all columns:
apply(apply(combn(listing(x), 3),2,diff)
If the two differences are equal, they fullfil the third property these numbers should have. Therefor we need to find the differende between the two differences. It that is 0, we have our result:
apply(apply(combn(listing(x), 3),2,diff),2,diff)==0])
That gives a logical vector, that is TRUE whenever the difference between the three numbers is the same. We get the numbers by:
combn(listing(x), 3)[,apply(apply(combn(listing(x), 3),2,diff),2,diff)==0]
I’ll define that as a function:

corrComb <- function(x){
  combn(listing(x), 3)[,apply(apply(combn(listing(x), 3),2,diff),2,diff)==0]  
}

It would be nice to have that as a logical test as well. I’m running out of ideas for reasonable function names:

test <- function(x){
  as.logical(sum(corrComb(x)))
}

Very simple – if there is a combination fulfilling the requirements, return TRUE.

Now, lets put it all together. I’ll use the library purrr, which allows me to use these nifty pipes.

Taking the candidate primes, we keep the ones that have three or more permutations that are prime, and in our candidate list. Of those, we keep the ones that have fulfill the requirement that the pairwise differences between three of them are equal.
Those are piped to corrComb, that returns the three numbers. They are piped to sort, so we get them in ascending order, and then I paste them together.

library(purrr)
answer <- cand %>%
  keep(poss) %>%
  keep(test) %>%
  corrComb%>%
  sort %>%
  paste(.,collapse="")

Lessons learned.

  1. Those pipes are pretty neat.
  2. The apply family is also pretty cool
  3. Read the documentation. I spent some time being confused about the results from permuations() until I discovered the set-flag

Euler 80

Problem 80 from Project Euler.

The problem tells us that if the square root of a natural number is not an integer, it is irrational.
Project Euler also claims that it is well known. I did not know it.

We are then told that the square root of 2 is 1.4142135623730950… And that the digital sum of the first 100 digits is 475.

The task is now to take the first 100 natural numbers. And find the total of the digital sums for the first 100 digits for all the irrational square roots.

Lets begin by figuring out how to handle that many digits. R does not support more than around 15 places after the decimal point.

The library Rmpfr can handle arbitrary precision:

library(Rmpfr)
## Loading required package: gmp
## 
## Attaching package: 'gmp'
## The following objects are masked from 'package:base':
## 
##     %*%, apply, crossprod, matrix, tcrossprod
## C code of R package 'Rmpfr': GMP using 64 bits per limb
## 
## Attaching package: 'Rmpfr'
## The following objects are masked from 'package:stats':
## 
##     dbinom, dnorm, dpois, pnorm
## The following objects are masked from 'package:base':
## 
##     cbind, pmax, pmin, rbind
a <- sqrt(mpfr(2,500))

The variable a now contains the square root of 2 with a precision of 500 bytes. I’m not quite sure how many decimal places that actually translates to. But testing the following code allows me to conclude with confidence that it is at least 100.

A thing to note here is, that

a <- mpfr(sqrt(2),500)

and

a <- sqrt(mpfr(2,500))

are not equal. In the first exampel sqrt(2) is evaluated before saving the value with the high precision. Start by converting the number 2 to a high precision representation, before doing math on it.

Next is writing a function that will return the digital sum of the first 100 digits of a number.

digitsum <- function(x){
  s <- 0
  for(i in 1:100){
    s <- s + floor(x)
    x <- (x - floor(x))*10
  }
  s
}

First s is initialized to 0. Then floor(x) gives us the first digit in x. We add that to s, and subtract it from x, and multiply by 10. Repeat that 100 times, and you get the sum of the first 100 digits in x.

Let us test that it works. Project Euler told us what the result for sqrt(2) is:

digitsum(a)
## 1 'mpfr' number of precision  500   bits 
## [1] 475

Nice, the correct result (not that that guarantees that I’ve done everything correctly).

Now, lets find all the irrational square roots we need to look at:

library(purrr)
t <- 1:100
s <- t %>%
  keep(~as.logical(sqrt(.x)%%1))

I need to practice this way of coding a bit more. t contains the first 100 natural numbers. I pass that to the keep()-function, and the predicate function takes the square root of each number, take the modulus 1, and convert it to a logical value. If the modulus of the square root is 0, the square root is an integer, and 0 i false. So we’re keeping all the non-integer squareroots.

Now I’ll convert all the natural numbers to the mpfr-class. The next line takes the square root. The third line calculate the digitalsum. And the final line gives us the result.

s <- mpfr(s,500)
r <- sqrt(s)
r <- digitsum(r)
sum(r)
## 1 'mpfr' number of precision  500   bits 
## [1] Censored

Lessons learned:
Rmpfr allows us to work with (more or less) arbitrary precision.
But we need to convert numbers to the relevant class before doing math on it.