Tag Archives: brute force

Euler 91 – Right Triangles

Blue Abstraction - Triangle

I haven’t done a Euler problem for a while. This one caught my eye, so I gave it a try in Haskell. Yeah, yeah, yeah, I know the picture is an equilateral…

My solution started out as a brute force, but then I added a few smarts along the way. First I ran a list comprehension of all possible points, then created a routine with guards to eliminate obvious impossible situations. The angle of each was calculated using the inverse cosine of two vectors:

550px-Find-the-Angle-Between-Two-Vectors-Step-6

… then its was simply a matter of checking for a 90 degree angle && the other two angles equalling 90. Oh yeah and then dividing the total by two because the list comprehension finds doubled matches because of transposed points.

I compiled this routine using GHC with flag -O2 and it runs in about 3 seconds for [0..50]. I’m still waiting for GHCi to finish …


result :: Int -> [Int]
result n = [1 | px <- r, py <- r, qx <- r, qy <- r, isRight px py qx qy ]  
	where
		r = [0..n]

isRight :: Int -> Int -> Int -> Int -> Bool
isRight px py qx qy
	| (px,py) == (0,0) = False
	| (qx,qy) == (0,0) = False
	| (px,py) == (qx,qy) = False
	| ao == 90 && (round (ap + aq)) == 90 = True
	| ap == 90 && (round (ao + aq)) == 90 = True
	| aq == 90 && (round (ao + ap)) == 90 = True
	| otherwise = False
	where
		ao = calcAngle px py qx qy
		ap = 180 + (-1 * (calcAngle px py (qx-px) (qy-py)))
		aq = 180 + (-1 * (calcAngle qx qy (px-qx) (py-qy)))

calcAngle :: Int -> Int -> Int -> Int -> Double
calcAngle px py qx qy = acos (num / denom) * 180 / pi
	where
		num = fInt (px * qx + py * qy)
		denom = (calcMag px py) * (calcMag qx qy)

calcMag :: Int -> Int ->  Double
calcMag x y = sqrt $ fInt (x*x) + fInt (y*y)

fInt = fromIntegral

main = do
	putStrLn "Enter limit:"
	limit <- getLine
	let answer = round $ 0.5 * fInt (length $ result $ read limit)
	putStrLn ("Answer = " ++ show answer)

Euler 221 – Alexandrian Numbers


When I came upon this Euler problem, and saw it was based on my namesake, I couldn’t help going for it. As with many of the Project Euler problems, the trick is to find a technique that reduces the brute force solution into something more efficient.

The functions from the original problem can be re-written to:
A = p(p + d)((p^2 + 1)/d + p)

The Alexandrian numbers are calculated from the new function, where d runs over divisors of p^2+1 and p runs over all positive integers. Used a fast Divisor routine that I Googled, then started building a Set (which is speed optimized to avoid duplicates) until length is 15K. Convert to a List, then sort and display 150,000th number. Numbers 1 thru 6 also listed to verify computation working. Runs in about 10 mins on my quad core iMac.

# e221.py

from math import sqrt
from functools import reduce


def appendEs2Sequences(sequences, es):
    result = []
    if not sequences:
        for e in es:
            result.append([e])
    else:
        for e in es:
            result += [seq + [e] for seq in sequences]
    return result


def primefactors(n):
    i = 2
    while i <= sqrt(n):
        if n % i == 0:
            l = primefactors(n // i)
            l.append(i)
            return l
        i += 1
    return [n]


def factorGenerator(n):
    p = primefactors(n)
    factors = {}
    for p1 in p:
        try:
            factors[p1] += 1
        except KeyError:
            factors[p1] = 1
    return factors


def divisors(n):
    factors = factorGenerator(n)
    divisors = []
    listexponents = [list(map(lambda x:k ** x, list(range(0, factors[k] + 1)))) for k in factors.keys()]
    listfactors = reduce(appendEs2Sequences, listexponents, [])
    for f in listfactors:
        divisors.append(reduce(lambda x, y: x * y, f, 1))
    divisors.sort()
    return divisors


if __name__ == "__main__":
    A = set()
    p = 1
    while len(A) <= 1500001:
        for d in divisors((p * p) + 1):
            A.add(int(p * (p + d) * (p + ((p * p) + 1) / d)))
        p += 1
    L = list(A)
    L.sort()
    for i in [1,2,3,4,5,6,150000]:
        print (str(i) + ": " + str(L[i-1]))

Euler Problem 105 – Disjoint Sets

The interesting thing about this one was that I had to figure out how to load a file in Haskell by using IO inside a Monad. Maybe something pretty basic, but still a novel concept. And a little weird.

This solution could be considered a brute force solution, but it still runs in about 45 seconds on an old PC. The routine first reads the file, splits the lines, then converts the strings into Integers.

For each set, it calculates the perfect series of sub-sets (tricky understanding disjoint sets), then creates a test series against the criteria. Then it simply sums up each set if the sub-set counts are the same.

Now to solve for related problems 103 and 106.

import Data.List
import Data.List.Split

testSubs :: [Integer] -> [Integer] -> Bool
testSubs b c
	| sb /= sc && lb > lc && sb > sc = True
	| sb /= sc && lb <= lc = True 
	| otherwise = False
	where
		sb = sum b
		sc = sum c
		lb = length b
		lc = length c

isSpecial :: [Integer] -> Bool
isSpecial set
	| perfectSet == testSet = True
	| otherwise = False
	where
		list = init $ tail $ subsequences set
		perfectSet = [(b,c)| b <- list, c <- list, b /= c]
		testSet = [(b,c)| b <- list, c <- list, b /= c && testSubs b c]
		
main = do
	raw <- readFile "sets.txt"
	print $ sum [sum $ clean x | x <- [splitOn "," y | y <- lines raw], isSpecial $ clean x]
	where
		clean = map read

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 145 – Haskell Revisit

Oh my head. Is that toast I smell? I’m actually getting headaches from trying to understand Haskell. Functional programming is totally different fron Imperative. But its also very powerful, especially for mathematically inclined Euler problems. That being said, what the hell is a Monad?

Admittedly I found much of the following code online, but I changed it a bit for my purposes. I compiled it using GHC on a dual-core 6-year old PC and it ran in 13 minutes, 11 seconds! For those of you playing at home, that’s a 819% speed improvement over Python! (see previous post). Where’s the Advil?

import Data.Time

oddDigits :: Integer -> Bool
oddDigits n
	| n < 10 = mod n 2 == 1
	| otherwise = (mod (mod n 10) 2) == 1 && oddDigits (div n 10)

rev :: Integer -> Integer -> Integer
rev n acc
	| n < 10 = n + (10 * acc)
	| otherwise = rev (div n 10) ((mod n 10) + (10 * acc))

isReverse :: Integer -> Bool
isReverse n
	| mod n 10 == 0 = False
	| otherwise = oddDigits (n + (rev n 0))

main = do
	start <- getCurrentTime
	print (length (filter isReverse [1..10^9]))
	stop <- getCurrentTime
	print $ diffUTCTime stop start

Euler 145 – Brute Force Solution

Testing for reversible properties of a billion numbers sounds like a simple premise. And it is. But it takes time. There are a number of shortcuts you can find on the Internet, like if a number is reversible then its reverse is also reversible. But after messing around with whether or not a number and its reverse are part of a set with certain starting or ending numbers etc, so as to avoid even numbers, I came to the conclusion that good old-fashioned brute force gets the job done best.

The routine below runs in about 1.8 hours on a 6-year old dual-core machine. Not blazingly fast. Python is SLOOOOW. I have a need for speed, but hate C/C++. Hmmm, being a sucker for punishment, I’m going to translate the same code into Haskell and see what happens. Haskell compiler speeds have a reputation of being epic. We shall see…

import time
if __name__=="__main__":
    t0 = time.time()
    c = 0
    i = 1
    while i < 10**9:
        if i % 10**6 == 0:
            print "i:"+str(i)+" -> c:"+str(c)
        if i % 10 == 0:
            i += 1
            continue
        rev = str(int(str(i))+int(str(i)[::-1]))
        if '0' in rev or '2' in rev or '4' in rev or '6' in rev or '8' in rev:
            i += 1
            continue
        c += 1
        i += 1
    print "Answer = "+str(c)
    print time.time() - t0, "seconds"

Euler Problem 92

This one takes a while to run, but that’s no surprise since I used brute force. Nothing too complicated, just iterate 10 million time and keep a counter. But if you use psyco, it does speed things up a bit. Roughly 10 mins.

#Euler92
import psyco
psyco.full()

class Problem:
    def digits(self, n):
        sum = 0                
        for c in str(n):
            sum += int(c)*int(c)
        return sum
    def Solution(self):
        count = 0
        for i in range (1, 10000001):
            if i % 100000 == 0:
                print i, count
            num = i
            f1 = 0
            f89 = 0             
            while f89 < 2 and f1 < 2:
                num = self.digits(num)
                if num == 89:
                    f89 += 1
                    if f89 == 2:             
                        count += 1
                        break
                elif num == 1:
                    f1 += 1
                    if f1 == 2:
                        break
        return "Answer = " + str(count)
if __name__ == '__main__':
    P = Problem()
    print P.Solution()