Tag Archives: prime numbers

Basket Making 101

Kalathos_Louvre_MNE1176
Praxis 2: Write a function that calculates the list of prime numbers using the Sieve of Eratosthenes.

Ummm, yeah. Done it a couple dozen times before over at Project Euler. Never hurts to review. I did manage to use more function composition this time around. Oh, and don’t forget to install module data.ordlist using Cabal.


import Data.List.Ordered

primes :: Integer -> [Integer]
primes m = 2 : sieve [3,5..m]
    where
        sieve [] = []
        sieve (p:xs) = p : sieve (minus xs [p*p, p*p+2*p..m])

main = do
    putStrLn "Enter n: "  
    input <- getLine
    putStrLn $ "# of primes: " ++ (show . length . primes $ read input)

Euler Problem 243 – Always the Totient

OK, so I had to look up the answer to this one. In my defence I did lots of research and determined that the resilience is related (once again) to Euler’s Totient function, specifically the relationship R(d) = phi(d)/(d-1). I figured that by using this relationship (and looking for R(d) lower than the ratio given) was the trick instead of having to brute force it. Nope. That is the brute force method. Ran for days and got nowhere.

So I broke down and looked up a solution. There aren’t very many for this one. Most were in interative languages like Ruby, which don’t help me. I did find one amazing solution is Haskell here, but for the life of me I don’t know why it works. I’ve adapted it and replaced as much as I can with clearer functions to help explain it. Perhaps not as efficient, but it still runs in less than a second.

Step One: Generate lists of primes and totients.
Step Two: Generate a list primorials consisting of tuples based on products of primes and totients.
Step Three: Use a recursive routine starting with the primorials list and a multiplier (m) of 1. If the current list head * m == the next primorial, then drop the head and loop. If the current list head * m is than the target then return it, otherwise increment m and loop.
Step Four: Shrug shoulders, open arms with palms upward, crinkle forehead, take a deep sign, and say, “WTF!”

import Data.List.Ordered

primes :: [Integer]
primes = 2 : sieve [3,5..] where
	sieve [] = []
	sieve (p:xs) = p : sieve (xs `minus` [p*p, p*p+2*p..])

totient :: Integer -> Integer
totient n = fromIntegral $ length [x | x <- [1..n], gcd x n == 1]

primorials :: [(Integer, Integer)]
primorials = [(product $ take x primes, product $ map totient $ take x primes) | x <- [1..]]

findRes :: [(Integer, Integer)] -> Integer -> Integer
findRes ((cpri,cphi):(npri,nphi):pps) m
	| cpri * m == npri = findRes ((npri,nphi):pps) 1
	| fromIntegral (cphi * m) / fromIntegral ((cpri * m) - 1) < ratio = cpri * m
	| otherwise = findRes ((cpri,cphi):(npri,nphi):pps) (m + 1)
	where ratio = 15499/94744

main = print $ findRes primorials 1

Euler Problem 123 – Remainder of Primes

Here is another prime number problem. Again, I re-used my prime number sieve, which I called whenever a remainder was to be calculated. I suppose I could have sped things up a bit by pre-calculating the primes, but this routine still runs in only a few minutes and does the trick.

Haskell’s list manipulation capabilities are simply amazing. The one thing I wasn’t sure about was if asking for the list head in main was sufficient to end the routine. Note that the list comprehension is based on [1,3..] which is open ended. Obviously, it did end, otherwise you wouldn’t be reading this and I’d still be staring at the prompt.

import Data.List.Ordered

pri :: [Integer]
pri = 2 : sieve [3,5..] where
    sieve [] = []
    sieve (p:xs) = p : sieve (xs `minus` [p*p, p*p+2*p..])

calcRem :: Integer -> Integer
calcRem n = (((pri!!(x-1))-1)^x + ((pri!!(x-1))+1)^x) `rem` (pri!!(x-1))^2
	where x = fromIntegral n
	
main = print $ head [x | x <- [1,3..], calcRem x > 10^10]

Euler Problem 124 – Prime Factors

When I was looking for a picture to represent Prime Number Factorization, I had no idea what to use. So I Googled images on the topic and discovered there was a Star Trek: Voyager episode called Prime Factors. The above picture is a scene from it — exciting I know! 😉 Too bad its not a picture of Seven of Nine.

So my solution consist of two parts. A prime number generator. I re-used the code from my solution of Euler Problem 131. The other part is a factorization routine that matches against the list of prime numbers. I found some concise Haskell code here.

This routine runs in less than 4 seconds on a single core 3 Ghz system.

import Data.List
import Data.List.Ordered

primes :: Integer -> [Integer]
primes m = 2 : sieve [3,5..m] where
    sieve [] = []
    sieve (p:xs) = p : sieve (xs `minus` [p*p, p*p+2*p..m])

primeFactors :: Integer -> [Integer]
primeFactors x = 
	unfoldr findFactor x where
		first (a,b,c) = a
		findFactor 1 = Nothing
		findFactor b = (\(_,d,p)-> Just (p, d))
			$ head $ filter ((==0).first) 
			$ map (\p -> (b `mod` p, b `div` p, p)) $ primes 100000

rad :: Integer -> Integer
rad n = product $ Data.List.nub $ primeFactors n

main = print $ snd $ sort [(rad x, x) | x <- [1..100000]] !! (9999)

Ah, what the heck — Why not?

Euler Problem 131 – Prime Numbers

Now that I’m on a Haskell kick, I opted to solve this one because of the need for calculating prime numbers. I found tons of resources online that detail different ways to calculate it. I settled on a modified Sieve of Eratosthenes that also takes into account squares, making it faster. Calculating the 78K prime numbers under a million, only takes a few seconds.

The other part of the solution comes from the fact that by playing with the values of the relationship from the problem, we see that n^2 and (n+p) have to be cubes. This makes the prime number a difference of cubes such that prime = (a+1)^3–a^3. This new value stays under a million when a is less than 600. So we just test over this smaller series and count the number of resulting primes. This solution runs in less than 10 secs.

Oh by the way, just for kicks the above image represents the distribution of prime numbers. The source code can be found here.

import Data.List
import Data.List.Ordered

primes :: Integer -> [Integer]
primes m = 2 : sieve [3,5..m] where
    sieve [] = []
    sieve (p:xs) = p : sieve (xs `minus` [p*p, p*p+2*p..m])

testEquation :: Integer -> Bool
testEquation a = (((a+1)^3)-(a^3)) `elem` (primes 1000000)

main = print $ length [(((a+1)^3)-(a^3))| a <- [1..600], testEquation a ]

Euler Problem 87 – More Primes

I started this one by iterating 50 million times. The code worked but it was damn slow. I figured it would take 27 hours. Not good. Then … lightbulb! Flip the question over by iterating on only the possible primes. Simple, elegant, fast — works in 4 seconds.

#Euler87
import math
class Problem: 
    def primeGen (self, lowerLimit, upperLimit):
            isPrime = [False] + [True]*(upperLimit-1)
            p = 2
            while p * p < upperLimit:
                if isPrime[p] == True:
                    m = p * p
                    while m < upperLimit:
                        if m % p == 0:
                            isPrime[m] = False
                        m += p
                p += 1
            P = []
            count = 0
            for p in isPrime:
                if p and count >= lowerLimit :
                    P.append(count)
                count += 1
            return P    
    def Solution(self):
        P2 = self.primeGen(2,int(math.exp(.5*math.log(50e6))))
        P3 = self.primeGen(2,int(math.exp(math.log(50e6)/3)))
        P4 = self.primeGen(2,int(math.exp(.25*math.log(50e6))))
        print "Number of Primes for ", "2^:", len(P2), "3^:", len(P3), "4^:", len(P4)
        N = []
        for p2 in P2:
            for p3 in P3:
                for p4 in P4:
                    if p2**2 + p3**3 + p4**4 <= 50e6:
                        N.append(p2**2 + p3**3 + p4**4)
        print "\nAnswer =", len(set(N))
if __name__ == '__main__':
    P = Problem()
    P.Solution()

Euler Problem 77

So this one is simply a variant of the previous problem, but with the additional of a prime number filter. Its necessary to work backwards until you have a count that closest to 5000, without going under.

#Euler77
class Algorithms:
    def partitionsGen(self, n):
        ms = {n: 1}
        keys = [n]
        yield ms
        while keys != [1]:
            if keys[-1] == 1:
                del keys[-1]
                reuse = ms.pop(1)
            else:
                reuse = 0
            i = keys[-1]
            newcount = ms[i] = ms[i] - 1
            reuse += i
            if newcount == 0:
                del keys[-1], ms[i]
            i -= 1
            q, r = divmod(reuse, i)
            ms[i] = q
            keys.append(i)
            if r:
                ms[r] = 1
                keys.append(r)
            yield ms
    def primeGen (self, lowerLimit, upperLimit):
        isPrime = [False, True]
        i = 2
        while i < upperLimit:
            isPrime.append(True)
            i += 1
        p = 2
        while p * p < upperLimit:
            if isPrime[p] == True:
                m = p * p
                while m < upperLimit:
                    if m % p == 0:
                        isPrime[m] = False
                    m += p
            p += 1
        P = []
        count = 0
        for p in isPrime:
            if p and count >= lowerLimit :
                P.append(count)
            count += 1
        return P
class Problem77:
    def __init__(self, limit):
        self.limit = limit
        self.Euler = Algorithms()
    def Solution(self):
        answer = 0
        Primes = self.Euler.primeGen(2,self.limit+1)
        print Primes
        while answer == 0:
            count = 0
            print "Prime Partitioning "+str(self.limit)
            for part in self.Euler.partitionsGen(self.limit):
                flag = True
                for n in part:
                    if n not in Primes:
                        flag = False
                if flag == True:
                    count += 1
            print "Count = "+str(count)
            if count < 5000:
                answer = self.limit + 1
            self.limit -= 1
        return "Answer = " + str(answer)
if __name__ == '__main__':
    P = Problem77(73)
    print P.Solution()