Solution to ps1

[jetamors]

The solution I wrote is behind the cut. Please let me know if you have any questions about the math or the code itself, or feel free to share your own solution in the comments.

# ps1 for MIT OpenCourse 6.00 Introduction to Computing
# Written by Jetamors, January 2010

# Problem 1: compute and print the 1000th prime number

# Remember that a prime number is an integer (whole number) that is not evenly divisible by any integer
# other than itself and 1.  Therefore, my strategy is to divide every number bigger than 2 by all the integers smaller
# than it, except 1. For example, the number 5 would be divided by 4, 3, and 2.
# If every division gives me a remainder, then the number is prime; if even one remainder is zero, then the
# number cannot be prime.
# 
# There are two things I did to simplify this problem:
# 1. Every even number bigger than 2 cannot be a prime number; therefore I stripped out all the prime numbers.
# 2. A number can't be evenly divisible by any integer more than half its size.  For example, the biggest number you can
#    evenly divide 10 by is 5. (I'm not really sure how to explain the math behind this, but if you test it it'll be
#    true.)  Therefore, I don't bother dividing a potential prime number by any number more than half its size.
# You don't have to do either of these things, but your code will take much longer to run if you don't do them. This code
# could be made to run even faster by immediately going to the next candidate as soon as you get a zero remainder,
# but... I didn't do that, oh well.

print 'Problem 1:'
primeiter = 1 # This counts the number of primes. We initialize this to 1 to account for the prime number 2.
candidate = 2 # This is the variable we use for each candidate for prime number. This is initialized to 2 so the first number
              # the while loop evaluates will be 3.

while (primeiter<1000):
    candidate = candidate + 1                               # increment at the beginning so the last number won't be incremented
    # print 'Current iteration is',candidate
    remaindercheck = 1                                      # this is the variable we use when evaluating remainders
    if (candidate/2)*2 != candidate:                        # this strips out all the even numbers
        for i in range(3,candidate/2):                      # this strips out all the numbers bigger than 1/2.
            remaindercheck = remaindercheck * candidate%i   # This number is only non-zero if every remainder is non-zero.
            # print remaindercheck
        if remaindercheck != 0:
            # print 'Prime number is',candidate
            primeiter = primeiter + 1                       # increment this every time we see a prime number

# printing the output
print 'The 1000th prime number is', candidate
print ' '





# Problem 2: compute the sum of the logarithms of all the primes from 2 to some number n.
#            Print out the sum of the logs of the primes, the number n, and the ratio of these two quantities.

# Here, I modified the code from problem 1. The major difference is in the while loop.  Before, the
# while loop continued until we reached the 1000th prime number.  Now, we want the loop to continue until
# n is reached.  Notice also that we add the log of a number only if it's a prime number.

print 'Problem 2:'
from math import *      # so we can use the log function
candidate = 2 
logprimesum = log(2)    # this initalizes the sum of the logs with the log of the prime number 2.
n = int(raw_input('Please enter a number! '))

while (candidate<=n):   # note that the while statement has changed; now it's based on n, not the number of primes.
    candidate = candidate + 1
    remaindercheck = 1
    if (candidate/2)*2 != candidate:
        for i in range(3,candidate/2):
            remaindercheck = remaindercheck * candidate%i
        if remaindercheck != 0:
            logprimesum = logprimesum + log(candidate)
            
# printing the output
print ' '
print 'The sum of the logs of the prime numbers from 2 to',n,'is',logprimesum
print 'The ratio of that sum to',n,'is',logprimesum/n

Comments

[winterthunder]

Thank you for posting this. :) I was trying to do something with the sieve of Eratosthenes, and it wasn't working in a pretty spectacular fashion.

Reading through this, I'm following you up until this line: remaindercheck = remaindercheck * candidate%i

If I've understood correctly, this is the line that divides the number by all the numbers less than it, right? How exactly is it doing that?

[jetamors]

Sorry, I see I didn't give a good explanation for that part! C&Ping the relevant section:

for i in range(3,candidate/2):                      # this strips out all the numbers bigger than 1/2.
            remaindercheck = remaindercheck * candidate%i   # This number is only non-zero if every remainder is non-zero.

In case a little extra explanation is needed: a for loop is a specific kind of while statement. It increments the variable i from 1 to candidate/2, and then it finishes. (Generally speaking, I've been taught that it's better to use for loops where you can; a for loop is guaranteed to end at some point, but an incorrectly written while loop can loop infinitely.)

i is the number I'm going to be dividing candidate by. The statement candidate%i returns the remainder of candidate/i; so for example, if candidate = 5 and i = 3, then candidate%i = 2, because 5/3 = 1 and 2/3.

For the statement remaindercheck = remaindercheck * candidate%i , I'm taking advantage of the fact that any number multiplied by zero is equal to zero. What this means is that if any remainder is equal to zero, then remaindercheck will be set to zero until the end of the for loop. If none of the remainders are zero, then remaindercheck will be some positive number.

For example:
if candidate = 9:
when i = 3: remainder of 9/3 = 0, therefore remaindercheck = 1 * 0 = 0
when i = 4: remainder of 9/4 = 1, therefore remaindercheck = 0 * 1 = 0

When the for loop ends, remaindercheck = 0 and we know this is not a prime number.

if candidate = 11:
when i = 3: remainder of 11/3 = 2, therefore remaindercheck = 1 * 2 = 2
when i = 4: remainder of 11/4 = 3, therefore remaindercheck = 2 * 3 = 6
when i = 5: remainder of 11/5 = 1, therefore remaindercheck = 6 * 1 = 6

When the for loop ends, remaindercheck != 0 and we know this is a prime number.

Does that make sense?

[winterthunder]

I think so... Or it's getting closer to making sense, at least. :)

Are you defining i in the line: for i in range(3,candidate/2)?

Thanks so much for posting and being willing to explain!

[jetamors]

Yes, that's correct. If it helps, here's a while loop that would give exactly the same result as my for loop:

i = 3
while (i < candidate/2):                      
    remaindercheck = remaindercheck * candidate%i
    i = i + 1

Let me know if you have any other questions or want me to explain something in a different way. It's so hard to explain math when we can't talk face to face >.<

[winterthunder]

Oooh, that does help. It also tells me that I need to go back to the while loop/for loop sections to get a better handle on how the two are constructed and the differences between them.

Yeah, I'm pretty sure at least half of my issues with this is the math; I got through AP Calc BC by the skin of my teeth in high school, six years ago, and I haven't looked at anything other than statistics since then. Thank you again!

[jetamors]

Oh and now having looked up the sieve of Eratosthenes... that would work fine for problem 2, but I think it would be pretty difficult for problem 1. The trouble is that the sieve assumes that you have a specific largest integer in mind: for example, all the primes between 1 and 300. But for problem 1, we're trying to find our largest integer: the 1000th prime number.

You could probably kludge it by getting, say, every prime number between 1 and 10,000 and then hoping that the 1000th prime is in there. But that would be pretty inefficient.

[winterthunder]

I think my thought process with the sieve was that you could use the concept of eliminating all multiples of a prime to select the numbers to test. I never actually got to workable code, but maybe if you created an infinite list and used the multiples to eliminate numbers from the list... then with each number that wasn't a multiple, run the remainder test and eliminate any number that fails it, and then when your list reaches 1000 items you just ask for the last one. Except, now that I type this, I'm pretty sure that lists have to be finite, and even if they weren't, your list would never actually reach 1000 because you'd still have all the untested numbers above the thousandth integer in it. :P

[jetamors]

That makes sense in terms of the mathematical theory underlying it, I think. Unfortunately, as you realized, in programming we can't just write an elegant equation and call it a day; ultimately, any computer can only deal with actual numbers, so you can never do an infinite anything.

If we used your approach with a finite list, though, it would still work: for example, you could take the numbers from 1 to 10,000 and sieve out all the non-primes. You'd be left with more than 1000 prime numbers, so you could just print the 1000th. It's not how I would do it personally, but the cool thing about programming is that even for simple problems like this, different people will write very different programs to solve them. There's more personal style in it than you might think :)

[aranthe]

Since primes require so much computing power, I tried to close down the number of candidates and tests as much as I could. To do that, I did three things:

• Since 2 is the first prime, generated only odd number candidates.
• Stored the primes and used a slice based on the prime nearest 1/3 of the candidate value as test divisors. (Since all the even numbers are eliminated, 1/3 becomes the upper limit.)
• Used the break keyword to exit the loop as soon as a candidate number was weeded out.

It runs pretty fast; here are the stats:

• Total candidates tested: 3959
• Total tests: 203,514
• Average tests per candidate: 51

# Initialize state variables.
primes = [2]
count = 1
test_slice = primes[0:]

# Do this until we have 1000 primes
while len(primes) < 1000:

    # Increment counter by 2. (Eliminates all divisible by 2.)
    count += 2
    
    # Set prime switch true.
    is_prime = True

    # Loop through the testing slice.
    for prime in test_slice:

        # Check to see if count is prime.
        if count % prime == 0:
            
            # Not a prime, set switch to false.
            is_prime = False

            # If divisible by 3, figure new test slice upper limit.
            if prime == 3:

                # Get the divisor limit.
                divisor_limit = count / 3

                # Check to see if it's in the prime array.
                # (If divisor_limit is a power of 3, it won't be.)
                if divisor_limit in primes:

                    # If it's a prime, find its index and add 1.
                    i = 1 + primes.index(divisor_limit)
                    # set the new test slice range from 1 to i.
                    test_slice = primes[1:i]

            # Don't waste any more time on this number.
            break

    # Check the switch; if it's still true...
    if is_prime:
        #...count is prime. Add it to the primes array.
        primes.append(count)

        # This is only for the first loop.
        if count == 3:
            test_slice = [count]

# print
ordinal = str(len(primes)) + 'th'
print 'The', ordinal, 'prime is', count

I'd love to see if I could improve on that, but I'm really out of time this week.

[jetamors]

Haha, that runs a lot faster than mine--good catch on the 1/3rd and using the break command.

[owl]

Here's my version. I decided to store the primes I'd already calculated in a list and use them as the divisors for the next candidate, rather than determining the divisors each time through the loop. Python's high-level, I imagine its list implementation is efficient enough for me not to worry about it :)
A number can't have a factor that's greater than its square root, so if it has reached that point, it's a prime and the control can break out of the loop.

#Problem 1
primes = [2] #list of primes. We know that 2 is the first prime
integer = 3 #so start checking at 3

while len(primes) < 1000:
   for p in primes: #we want to divide the candidate integer by all the primes we already have
       is_prime = 1 #we'll assume it's prime until proven otherwise
           if (integer % p == 0):
	       is_prime = 0 #if it can be divided equally by a prime, it's not prime
	       break        #so break out of the loop here
	   if p*p >= integer:
		break #don't keep checking with divisors greater than the square root
    if is_prime: #stopped dividing; is it prime?
        primes.append(integer) #add it to the array of primes
    integer = integer + 1 #and then move on to the next integer
#we're out of the while loop, so print the last prime, 
#and let's check how many we have while we're at it
print "The %dth prime is: %d" %(len(primes), primes[len(primes)-1]) #Don't forget the list is 
#zero-indexed, or we'll get a lovely index out of bounds error here

And the slightly adapted version for the second problem:

from math import *
n = int(raw_input("Enter a number: ")) #ask the user what is n
primes = [2] #list of primes. We know that 2 is the first prime
integer = 3 #so start checking at 3
sum_logs = log(2) #the log of the primes we have so far

while len(primes) < n:
    for p in primes: #we want to divide the candidate integer by all the primes we already have
        is_prime = 1 #we'll assume it's prime until proven otherwise
            if (integer % p == 0):
        is_prime = 0 #if it can be divided equally by a prime, it's not prime
        break        #so break out of the loop here
        if p*p >= integer:
            break #don't keep checking with divisors greater than the square root
     if is_prime: #stopped dividing; is it prime?
        primes.append(integer) #add it to the array of primes
        sum_logs = sum_logs + log(integer) #find the log and add it to the total
    integer = integer + 1 #and then move on to the next integer
#print out the sum of the logs, and the ratio
print "The sum of the logs of the first",len(primes), "primes is",sum_logs
print "The ratio of the sum to",n,"equals",(sum_logs / n)

I don't think I've put so many comments in code since I was in first year of uni ;)
You might also notice that I like to put double quotes around strings, and have a space on each side of the operators. I just think it's easier to read that way.

[jetamors]

This one also runs pretty fast! I do see an error in your part 2, though; the question asks you to find all the primes between 1 and n, not n number of primes. If you do it that way, the ratio you get should always be between 0 and 1.

[owl]

Ooops! Requirements reading fail on my part, you'd think I'd know better by now. If I change while len(primes) < n: to while integer <= n: that should do the trick.

[aranthe]

Very elegant! Since I used the same basic logic outline, I merged our approaches: Used your square root limit with my increment of 2 and knocked about another 4000 tests off. Cool!