#!/usr/bin/env python

# Above statement tells your terminal that the script should be run in 
# a python enviroment.

# --------------------------------------------------------- #
#  Script calculating all prime numbers up to 10000
#
#  Author:  Troels C. Petersen (Niels Bohr Institute)
#  Date:    06-10-2015 (originally)
# --------------------------------------------------------- #

from __future__ import print_function, division
import numpy as np               # Import the numpy library

Nmax = 10000              # Maximum prime number
N    = 2                  # Starting value


# Calculating prime numbers comes down to taking a number (here N),
# and testing if it is prime by dividing it by all numbers up the
# the square root of it. If no number divides N (tested using the 
# "%" operator, which returns the remainder: 6%3 = 0, 5%2 = 1, etc.),
# then N is prime!

# See if you can follow this in python program form below. There 
# are only 9 lines in total! The central part consists of a While-loop,
# which loops as long as the statement is true.
# It is true to begin with (N=2, Nmax = 10000), but at some point
# N will have increased enough to equal Nmax, and the looping stops.

# Create empty list for primes. We dont create a numpy list since we
# don't know the size of it before creation:
primes = []

while (N < Nmax):
    # print "Potential Prime number I want to test: ", N
    hit = False                # Define a variable "hit", with starting value "False".

    for i in range (2,int(np.sqrt(N))+1) :   # Loop over the possible divisors.
        # print "  I'm testing the first number with the second and get the third: ", N, i, N%i

        if (N%i == 0):         # If the remainder of the integer division is 0
            hit = True         # hit = "True" (i.e. the number being tested had a divisor, and is not a prime)
            break              # Skip the rest of the possible divisors.

    if (not hit) :             # If no number gave a perfect integer division (i.e. with remainder 0)...
        primes.append(N)
        print("I got a prime", N)   # ...the it is a prime number. Fill it into the histogram.

    N += 1                     # Increase N by one.

# All primes up to Nmax as a numpy array
primes = np.array(primes)
print()
print("First 10 primes, ", primes[:10]) # prints the first 10 primes

# Check that there are no even primes except 2:
print("  The even prime numbers are listed here: ", primes[(primes % 2) == 0])

# This algorithm gives the prime numbers up to Nmax = 10000!
# Think it through, and see if you understand it...