Skip to content

A little exercise in R!

8 messages · Philippe GROSJEAN, Petr Savicky, Bert Gunter +2 more

Original
(Ted Harding)
# 
Greetings all!
A recent news item got me thinking that a problem stated
therein could provide a teasing little exercise in R
programming.

http://www.bbc.co.uk/news/uk-england-cambridgeshire-17680326

  Cambridge University hosts first European 'maths Olympiad'
  for girls

  The first European girls-only "mathematical Olympiad"
  competition is being hosted by Cambridge University.
  [...]
  Olympiad co-director, Dr Ceri Fiddes, said competition questions
  encouraged "clever thinking rather than regurgitating a taught
  syllabus".
  [...]
  "A lot of Olympiad questions in the competition are about
  proving things," Dr Fiddes said.

  "If you have a puzzle, it's not good enough to give one answer.
  You have to prove that it's the only possible answer."
  [...]
  "In the Olympiad it's about starting with a problem that anybody
  could understand, then coming up with that clever idea that
  enables you to solve it," she said.

  "For example, take the numbers one up to 17.

  "Can you write them out in a line so that every pair of numbers
  that are next to each other, adds up to give a square number?"

Well, that's the challenge: Write (from scratch) an R program
that solves this problem. And make it neat.

NOTE: If there should happen to be some R package that can solve
this kind of problem already, without you having to think much,
then its use is illegitimate! (I.e. will be deemed "regurgitation").

Over to you.

With best wishes,
Ted.

-------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at wlandres.net>
Date: 13-Apr-2012  Time: 22:33:43
This message was sent by XFMail
Philippe GROSJEAN
# 
Hi all,

I got another solution, and it would apply probably for the ugliest one :-(
I made it general enough so that it works for any series from 1 to n (n 
not too large, please... tested up to 30).

Hint for a better algorithm: inspect the object 'friends' in my code: 
there is a nice pattern appearing there!!!

Best,

Philippe

..............................................<?}))><........
  ) ) ) ) )
( ( ( ( (    Prof. Philippe Grosjean
  ) ) ) ) )
( ( ( ( (    Numerical Ecology of Aquatic Systems
  ) ) ) ) )   Mons University, Belgium
( ( ( ( (
..............................................................

findSerie <- function (n, tmax = 500) {
   ## Check arguments
   n <- as.integer(n)
   if (length(n) != 1 || is.na(n) || n < 1)
     stop("'n' must be a single positive integer")
	
   tmax <- as.integer(tmax)
   if (length(tmax) != 1 || is.na(tmax) || tmax < 1)
     stop("'tmax' must be a single positive integer")
	
   ## Suite of our numbers to be sorted
   nbrs <- 1:n
	
   ## Trivial cases: only one or two numbers
   if (n == 1) return(1)
   if (n == 2) stop("The pair does not sum to a square number")
	
   ## Compute all possible pairs
   omat <- outer(rep(1, n), nbrs)	
   ## Which pairs sum to a square number?
   friends <- sqrt(omat + nbrs) %% 1 < .Machine$double.eps
   diag(friends) <- FALSE # Eliminate pairs of same numbers
	
   ## Get a list of possible neighbours
   neigb <- apply(friends, 1, function(x) nbrs[x])
	
   ## Nbr of neighbours for each number
   nf <- sapply(neigb, length)
	
   ## Are there numbers without neighbours?
   ## then, problem impossible to solve..
   if (any(!nf))
     stop("Impossible to solve:\n    ",
       paste(nbrs[!nf], collapse = ", "),
       " sum to square with nobody else!")
	
   ## Are there numbers that can have only one neighbour?
   ## Must be placed at one extreme
   toEnds <- nbrs[nf == 1]
   ## I must have two of them maximum!
   l <- length(toEnds)
   if (l > 2)
     stop("Impossible to solve:\n    ",
       "More than two numbers form only one pair:\n    ",
       paste(toEnds, collapse = ", "))
	
   ## The other numbers can appear in the middle of the suite
   inMiddle <- nbrs[!nbrs %in% toEnds]
	
   generateSerie <- function (neigb, toEnds, inMiddle) {
     ## Allow to generate serie by picking candidates randomly
     if (length(toEnds) > 1) toEnds <- sample(toEnds)
     if (length(inMiddle) > 1) inMiddle <- sample(inMiddle)
		
     ## Choose a number to start with
     res <- rep(NA, n)
		
     ## Three cases: 0, 1, or 2 numbers that must be at an extreme
     ## Following code works in all cases
     res[1] <- toEnds[1]
     res[n] <- toEnds[2]
		
     ## List of already taken numbers
     taken <- toEnds
		
     ## Is there one number in res[1]? Otherwise, fill it now...
     if (is.na(res[1])) {
         taken <- inMiddle[1]
         res[1] <- taken
     }
		
     ## For each number in the middle, choose one acceptable neighbour
     for (ii in 2:(n-1)) {
       prev <- res[ii - 1]
       allpossible <- neigb[[prev]]
       candidate <- allpossible[!(allpossible %in% taken)]
       if (!length(candidate)) break # We fail to construct the serie
       ## Take randomly one possible candidate
       if (length(candidate) > 1) take <- sample(candidate, 1) else
         take <- candidate
       res[ii] <- take
       taken <- c(taken, take)
     }
		
     ## If we manage to go to the end, check last pair...
     if (length(taken) == (n - 1)) {
       take <- nbrs[!(nbrs %in% taken)]
       res[n] <- take
       taken <- c(take, taken)
     }
     if (length(taken) == n && !(res[n] %in% neigb[[res[n - 1]]]))
     res[n] <- NA # Last one pair not allowed

     ## Return the series
     return(res)
   }
	
   for (trial in 1:tmax) {
     cat("Trial", trial, ":")
     serie <- generateSerie(neigb = neigb, toEnds = toEnds,
       inMiddle = inMiddle)
     cat(paste(serie, collapse = ", "), "\n")
     flush.console() # Print text now
     if (!any(is.na(serie))) break
   }
   if (any(is.na(serie))) {
     cat("\nSorry, I did not find a solution\n\n")
   } else cat("\n** I got it! **\n\n")
   return(serie)
}

findSerie(17)
On 13/04/12 23:34, (Ted Harding) wrote:
Petr Savicky
#
On Fri, Apr 13, 2012 at 10:34:49PM +0100, Ted Harding wrote:
Hi.

Is recursion acceptable? Using recursion, i obtained
two solutions.

  extend <- function(x)
  {
      y <- setdiff((1:17), x)
      if (length(y) == 0) {
          cat(x, "\n")
          return
      }
      y <- y[(y + x[length(x)]) %in% (1:5)^2]
      for (z in y) {
          extend(c(x, z))
      }
  }

  for (i in 1:17) extend(i)

  16 9 7 2 14 11 5 4 12 13 3 6 10 15 1 8 17 
  17 8 1 15 10 6 3 13 12 4 5 11 14 2 7 9 16 

Petr.
Bert Gunter
# 
Folks:

IMHO this is exactly the **wrong** way t go about this. These are
mathematical exercises that should employ mathematical thinking, not
brute force checking of cases.

Consider, for example, the 1 to 17 sequence given by Ted. Then 17
**must** be one end of the sequence and 16 the other. (Why?) Hence,
starting from the 17 end, the values ** must** be 17  8 1 ...
Proceeding in this way, it takes only a couple of minutes to solve.

The more interesting point which I think the question was really
about, is can this always be done? I haven't given this any thought,
but there may be an easy proof or counterexample. If the answer to
this latter is no, then perhaps even more interesting is to
characterize the set of numbers where it can/cannot be done.

But this is all way off topic, no?

Cheers,
Bert



On Fri, Apr 13, 2012 at 6:26 PM, Philippe Grosjean
<phgrosjean at sciviews.org> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Bert Gunter
    

    

      
      #
    

  


  

      
... and a moment's more consideration immediately shows it cannot be
done for n = 18, since 16,17, and 18 cannot all be at an end.

-- Bert

        
On Fri, Apr 13, 2012 at 9:59 PM, Bert Gunter <bgunter at gene.com> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      (Ted Harding)
    

    

      
      #
    

  


  

      
Well, this "Olympiad" challenge led to some interesting responses.
First, Bert Gunther noted that the arragement of 1:17 must have
17 at one end, allowing it to be solved on paper in a few minutes.
That would definitely be in the spirit of the Olympiad, where
'"In the Olympiad it's about starting with a problem that anybody
could understand, then coming up with that clever idea that enables
you to solve it," she [Ceri Fiddes] said.'

However, it was not in the spirit of my own challenge, which was
to write neat self-contained R code to solve it. So while Bert gets
special mention for clearing the pole-vault bar without a pole,
starting from a hand-stand position, there's no Gold there!

Justin Haynes posted it to stackoverflow, and gave the link to
Vincent Zoonekynd's answer. This recognised the underlying
graph-theoretic issue. However, his code wheeled on the igraph
package and used its graph.adjacency() function to solve it.
Neat, but flouted my rule that use of such an existing R package
was illegitimate. You are not allowed to win the 100 metres by
riding a motorcycle.

Petr Savicky's solution is the approach I preferred, and indeed
is along the lines of my own approach (and much neater than what
I achieved). Yes, recursion is indeed allowable, since it is built-in
to R's functionality. In this line of approach, starting from any
given integer and proceeding by selecting from allowable next
integers, you sooner or later encounter a multiple choice. So you
are tracing down the branches of a tree to see how far you can get. 

Some minor changes to Petr's code generalise it to being able to
answer Bert's later speculation about "can this always be done?",
i.e. is it just for 1:17, or what? Here is Petr's code as revised
(Petr's original commands, where changed, are preceded by "##",
and followed by my changes):

##extend <- function(x)
  extend <- function(x,N)
  {
##    y <- setdiff((1:17), x)
      y <- setdiff((1:N), x)    # for any 1:N, not just 1:17
      if (length(y) == 0) {
##        cat(x, "\n")
          cat(N) ; cat(": ") ;  cat(x, "\n")  # To show which N
          return
      }
##    y <- y[(y + x[length(x)]) %in% (1:5)^2]
      M <- ceiling(sqrt((N-1)+N))   # To include accessible squares
      y <- y[(y + x[length(x)]) %in% (1:M)^2]  # as here
      for (z in y) {
          extend(c(x, z),N)
      }
  }
##for (i in 1:17) extend(i)
##  16 9 7 2 14 11 5  4 12 13 3 6 10 15 1 8 17
##  17 8 1 15 10 6 3 13 12 4 5 11 14 2 7 9 16

Now to see what comes out for 1:N, N=1:17

  for( N in (1:17) ){for (i in (1:N)) extend(i,N)}
  # 1: 1
  # 15: 8 1 15 10 6 3 13 12 4 5 11 14 2 7 9
  # 15: 9 7 2 14 11 5 4 12 13 3 6 10 15 1 8
  # 16: 8 1 15 10 6 3 13 12 4 5 11 14 2 7 9 16
  # 16: 16 9 7 2 14 11 5 4 12 13 3 6 10 15 1 8
  # 17: 16 9 7 2 14 11 5 4 12 13 3 6 10 15 1 8 17
  # 17: 17 8 1 15 10 6 3 13 12 4 5 11 14 2 7 9 16

So, in addition to the original 1:17 case, we can also do it
for 1:15 and 1:16 (as well as the trivial 1:1). But how
far up can we take it? It turns out we can get suprisingly
prolific results! Take just going from 18:25 (which follows
on from the 1:17 above):

  for( N in (18:25) ){for (i in (1:N)) extend(i,N)}
  # 23: 2 23 13 12 4 21 15 10 6 19 17 8 1 3 22 14 11 5 20 16 9 7 18
  # 23: 9 16 20 5 11 14 22 3 1 8 17 19 6 10 15 21 4 12 13 23 2 7 18
  # 23: 18 7 2 23 13 12 4 21 15 10 6 19 17 8 1 3 22 14 11 5 20 16 9
  # 23: 18 7 9 16 20 5 11 14 2 23 13 12 4 21 15 10 6 19 17 8 1 3 22
  # 23: 18 7 9 16 20 5 11 14 22 3 1 8 17 19 6 10 15 21 4 12 13 23 2
  # 23: 22 3 1 8 17 19 6 10 15 21 4 12 13 23 2 14 11 5 20 16 9 7 18
  # 25: 2 23 13 12 24 25 11 14 22 3 1 8 17 19 6 10 15 21 4 5 20 16 9 7 18
  # 25: 3 22 14 2 23 13 12 4 21 15 10 6 19 17 8 1 24 25 11 5 20 16 9 7 18
  # 25: 4 21 15 10 6 19 17 8 1 3 22 14 2 23 13 12 24 25 11 5 20 16 9 7 18
  # 25: 8 17 19 6 10 15 21 4 12 13 23 2 14 22 3 1 24 25 11 5 20 16 9 7 18
  # 25: 9 16 20 5 4 21 15 10 6 19 17 8 1 3 22 14 11 25 24 12 13 23 2 7 18
  # 25: 10 15 21 4 12 13 23 2 14 22 3 6 19 17 8 1 24 25 11 5 20 16 9 7 18
  # 25: 11 25 24 12 13 23 2 14 22 3 1 8 17 19 6 10 15 21 4 5 20 16 9 7 18
  # 25: 13 23 2 14 22 3 1 8 17 19 6 10 15 21 4 12 24 25 11 5 20 16 9 7 18
  # 25: 18 7 2 23 13 12 24 25 11 14 22 3 1 8 17 19 6 10 15 21 4 5 20 16 9
  # 25: 18 7 9 16 20 5 4 21 15 10 6 19 17 8 1 3 22 14 2 23 13 12 24 25 11
  # 25: 18 7 9 16 20 5 4 21 15 10 6 19 17 8 1 3 22 14 11 25 24 12 13 23 2
  # 25: 18 7 9 16 20 5 11 25 24 1 3 22 14 2 23 13 12 4 21 15 10 6 19 17 8
  # 25: 18 7 9 16 20 5 11 25 24 1 8 17 19 6 3 22 14 2 23 13 12 4 21 15 10
  # 25: 18 7 9 16 20 5 11 25 24 1 8 17 19 6 10 15 21 4 12 13 3 22 14 2 23
  # 25: 18 7 9 16 20 5 11 25 24 1 8 17 19 6 10 15 21 4 12 13 23 2 14 22 3
  # 25: 18 7 9 16 20 5 11 25 24 12 4 21 15 10 6 19 17 8 1 3 13 23 2 14 22
  # 25: 18 7 9 16 20 5 11 25 24 12 4 21 15 10 6 19 17 8 1 3 22 14 2 23 13
  # 25: 18 7 9 16 20 5 11 25 24 12 13 23 2 14 22 3 1 8 17 19 6 10 15 21 4
  # 25: 22 14 2 23 13 3 1 8 17 19 6 10 15 21 4 12 24 25 11 5 20 16 9 7 18
  # 25: 23 2 14 22 3 13 12 4 21 15 10 6 19 17 8 1 24 25 11 5 20 16 9 7 18

So you also get solutions (non-unique) for 1:23 (6, basically 3,
of them) and for 1:25 (20, basically 10, of them). The case
  for( N in (26:100) ){for (i in (1:N)) extend(i,N)}
is left as an exercise for the reader.

Meanwhile, I await solutions from the teams who are competing in
the Circuit of the Universe, hurdling over all the (approx) 10^14.6
permutations of 1:17 (reports suggest that CERN may be in the lead).

Thanks for the responses, and best wishes to all,
Ted.

        
On 14-Apr-2012 04:59:48 Bert Gunter wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
-------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at wlandres.net>
Date: 14-Apr-2012  Time: 16:00:10
This message was sent by XFMail