Computation R program

问题:

I want to compute the smallest possible number that is divisible evenly by all the natural numbers from 1-20; I have written the following program in R and am not getting the desired output (rather it seems my loop is almost never ending).

My program is as follows:

a = 21
c = 0
while ( c < 20){
    c = 0
    n = 1
    while ( n < 21 ){
        if (a%%n == 0) c = c + 1
        n = n+1
    }
    a = a + 1
}
print (a)

Where did I go wrong?

回答1:

Here's a more R-like solution using that fact that the answer will be the product of primes, p <= 20, each raised to an index i such that p^i <=20

sMult <- function(x)
# calculates smallest number that 1:x divides
{
    v <- 1:x
    require(gmp) # for isprime
    primes <- v[as.logical(isprime(v))]
    index <- floor(log(x)/log(primes))
    prod(rep(primes,index))
}

Which yields:

> sMult(20)
[1] 232792560
> sMult(20)%%1:20
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

While this solution is general, it should be noted that for large x, isprime is probabalistic. Of course, when this is likely to give false results, you are also likely to have a number so large that it will not be able to be stored accurately. Fortunately the gmp package implements a large integer class, bigz. To use this change to final line of the function to:

prod(as.bigz(rep(primes,index)))

Compare the following results:

> sMult(1000)
[1] Inf
> sMult2(1000)
[1] "7128865274665093053166384155714272920668358861885893040452001991154324087581111499476444151913871586911717817019575256512980264067621009251465871004305131072686268143200196609974862745937188343705015434452523739745298963145674982128236956232823794011068809262317708861979540791247754558049326475737829923352751796735248042463638051137034331214781746850878453485678021888075373249921995672056932029099390891687487672697950931603520000"

回答2:

Its not clear what you are doing with n and c

Here's a refactored routine (I don't know R syntax, so you'll need to convert this, but the logic is still relavent)

For a number to be evenly divisable by all 1..20 it follows that it is also evenly divisable by all the primes <= 20 (obviously), so it must be divisable by the product of the primes <= 20 (which is 9699690)

So it is only necassary to test multiples of 9699690.

Start testing at 20 so the loop breaks sooner, fewer net iterations

You should add a overflow check on ans in case the answer is > max value of the data type of ans

ans = 9699690
Do
    Found = True
    For i = 20 To 2 Step -1
        If ans Mod i <> 0 Then
            Found = False
            ans = ans + 9699690
            Exit For
        End If
    Next
    If i < best Then best = i
    If Found Then Exit Do
Loop
Debug.Print ans

回答3:

Using Chris' logic to only test multitudes of 9699690 I can find the answer with:

found <- FALSE
test <- 9699690 

while(!found)
{
    test <- test + 9699690 
    found <- all(test%%(1:20)==0)
}

cat("The number is: ",test,"\n")
The number is:  232792560 

Edit:

As for why the OPs code isn't working, which is probably what you are interested in rather than solving this little puzzle, there is only one small problem with the code and it is probably not what you thing. If we enter one value before the true answer:

a = 232792559
c = 0 
while ( c < 20){     
    c = 0     
    n = 1     
    while ( n < 21 ){         
        if (a%%n == 0) c = c + 1         
        n = n+1     
    }     
    a = a + 1 
} 
print (a) 
[1] 232792561

We get one too much, which is because you add 1 to a even if the answer is correct. Move it to the front of the outer loop and it works. It just seems like an infinite loop because it is soooo slow in R to compute things like this. In R, we need to do things differently than languages like C because R code is not compiled (means loops take a long time) and is optimized for vectorized input.

Looking at your code tells me that you are new to R but experienced in other languages maybe, but this experience will not help you that much in R. I would advice you to read up on how to vectorize things.