Tag Archives: Math

Prime testing function in R

I was hoping to begin tinkering a bit with the multicore package in R beyond some extremely trivial examples.  Thanks to a combination of R’s dumb quirkiness (for example, being worthless on loops), my poor planning, and general bad programming, my Saturday afternoon tinkering project is ultimately worthless in fulfilling that purpose.

I was really hoping to take the prime searcher I had been using to solve Project Euler problems and use it to make a prime testing function which would utilize the 4 cores on my adorable little core i5.  I succeeded in creating it, but it wasn’t until I had finished (which thankfully was only about an hour of time wasted) that I realized that my plan was stupid and I was going to have to loop over entries (or use sapply(), which is the same thing; don’t kid yourself).  This, as we all know, is where R is fantastically bad.

A very hasty test showed that relying on a single core, the original way I was doing things is a little over 2 times faster than this new implementation.  This doesn’t bode well for the possible speedup when scaling up to 4 cores, and since I have more pressing projects at the moment, I’m going to have to dump this one barring some illumination.

However, it is a workable prime tester.  And for testing against single primes, it’s honestly not that slow.  But, as pointed out above, if you’re using it to test primality of a (reasonably large) range of (reasonably large) values, then you’re going to have to loop and R is going to catch fire in the process.

The syntax is simple enough; simply conjure your integer of choice n and ask R

IsPrime(n)

which is like the old Maple syntax back when I used it extensively (which was ~ Maple 8 or 9, I believe).  Functions and sample output below:

IsPrime <- function(n){ # n=Integer you want to know if is/not prime
  if ((n-floor(n)) > 0){
    cat(sprintf("Error: function only accepts natural number inputs\n"))
  } else if (n < 1){
      cat(sprintf("Error: function only accepts natural number inputs\n"))
    } else
  # Prime list exists
  if (try(is.vector(primes), silent=TRUE) == TRUE){

    # Prime list is already big enough
    if (n %in% primes){
      TRUE
    } else
    if (n < tail(primes,1)){
      FALSE
    } else
    if (n <= (tail(primes,1))^2){
      flag <- 0
      for (prime in primes){
        if (n%%prime == 0){
          flag <- 1
          break
        }
      }
    
      if (flag == 0){
        TRUE
      }
      else {
        FALSE
      }
    }

    # Prime list is too small; get more primes
    else {
      last.known <- tail(primes,1)
      while ((last.known)^2 < n){
        assign("primes", c(primes,GetNextPrime(primes)), envir=.GlobalEnv)
        last.known <- tail(primes,1)
      }
      IsPrime(n)
    }
  } else {
    # Prime list does not exist
    assign("primes", PrimesBelow(n,below.sqrt=TRUE), envir=.GlobalEnv)
    IsPrime(n)
  }
}

# Get next prime
GetNextPrime <- function(primes){ # primes=Known prime list
  i <- tail(primes,1)
  while (TRUE){
    flag <- 0
    i <- i+2
    if (i%%6 == 3){
      flag <- 1
    }
    if (flag == 0){
      s <- sqrt(i)+1
      possible.primes <- primes[primes<s]
      for (prime in possible.primes){
        if ((i%%prime == 0)){
          flag <- 1
          break
        }
      }
      if (flag == 0){
        break
      }
    }
  }
  i
}

# Primes below specified integer n; optionally only those below sqrt(n)
PrimesBelow <- function(n, below.sqrt=FALSE){
  if (below.sqrt == TRUE){
    m <- ceiling(sqrt(n))
  } else {
    m <- n
  }
  
  primes <- c(2,3)
  i <- 3
  while (i < m-1){
    flag <- 0
    i <- i+2
    if (i%%6 == 3){
      flag <- 1
    }
    if (flag == 0){
      s <- sqrt(i)+1
      possible.primes <- primes[primes<s]
      for (prime in possible.primes){
        if ((i%%prime == 0)){
          flag <- 1
          break
        }
      }
      if (flag == 0){
        primes <- c(primes, i)
      }
    }
  }
  primes
}

Some sample output:

 > primes
Error: object 'primes' not found
> IsPrime(100)
[1] FALSE
> IsPrime(101)
[1] TRUE
> primes
[1]  2  3  5  7 11
>

After screwing around with the prime tester for a few minutes, I was able to find this adorable little gem

1> IsPrime(1234567891)
[1] TRUE

So 1234567891 is prime, but none of 1, 12, 123, 1234, 12345, 123456, 1234567, 12345678, 123456789, 12345678910, 123456789101, or 1234567891011 are (of course, nearly half of those are even, but the point stands).


A Math Degree in the Kitchen

So today I was making some rice, and I decided based on calorie load that I wanted more than 1/3 cup of rice, but less than 1/2.  I decided that 5/12 cup was acceptable, but my only tools were a 1/4 cup scoop and a 1/3 cup scoop.

I didn’t have to think very hard about this; I instantly knew it was possible to get 5/12 cup of rice with these tools, because 4 and 3 are relatively prime–i.e., they share no common factors.  More on this in a minute, but first, how to do it.

It’s a very uncomplicated idea, and it should look very similar to the Diehard 3 water puzzle.  In the film, John McClain has a 3 liter jug and a 5 liter jug, and wants to get exactly 4 liters of water.  My rice puzzle is basically the same game.    To get my 5/12 cup of rice, I fill the 1/3 cup scoop and pour that rice into my rice cooker.  Then I again fill the 1/3 scoop and pour that (without spilling) into the 1/4 scoop.  How much remains in the 1/3 scoop?  Well, I had 4/12 cup of rice in it, and I took out 3/12 cup of rice.  So only 1 remains in the 1/3 cup scoop.  Add the 1/12 to the rice cooker, and since 1/3+1/12=5/12, we’re done.  And although that’s how I did it, I actually didn’t even need to get the rice cooker involved.  I can form 5/12 cups inside the scoops without having to use another vessel (though the thing holding all my rice is allowed).  I could have filled the 1/3 cup scoop, poured that into the 1/4 cup scoop, emptied the 1/4 cup scoop, poured the 1/12 cup in the 1/3 cup scoop into the 1/4 cup scoop, then filled the 1/3 cup scoop.  Tada!

I warned you that it was very uncomplicated; the above is really not very interesting on any level.  What is interesting is why this works, and how I knew so quickly that it must be possible to do it.  If you don’t know why, then all that stuff I just wrote might seem like a cute curiosity, when really there is a very basic underlying principle involved.  Instead of needing to ask “I wonder if it’s possible to get 1/6 cups of rice with those tools”, I can instantly say that the answer is “yes”.  I didn’t have to think about it at all.  Here’s how it works.

Instead of thinking of the scoops as 1/4 and 1/3 cup, I’ll do you a favor and eliminate the fractions for you.  Let’s think of them as 3 and 4 serving scoops (out of 12), respectively.  Believe it or not, we’re almost done.  As it turns out, if the greatest common divisor (written gcd) of these two values (in this case 3 and 4 from 1/4 cup and 1/3 cup scoops, respectively) is 1, then it’s possible to get any of 1/12, 2/12, 3/12, 4/12, 5/12, 6/12, or 7/12 cups of rice (maxing out at 3+4=7, because that’s the maximum amount of stuff our two scoops can hold).  Said another way, I can get 1, 2, 3, 4, 5, 6, or 7 servings (out of 12) in this fashion.  Why?  Because if gcd(3,4)=1 (and thank goodness that’s the case), then there exist integers x and y satisfying

3x+4y=1

But that isn’t what I wanted!  I wanted something like 3a+4b=5.  But this is easy.  Multiply the above equation on both sides by 5 and we get

5(3x+4y)=3(5x)+4(5y)=5

So just call 5x=a and 5y=b.  So it must be possible.

And people say that pure math has no practical application.