#    Copyright (C) 2008 Toby Thain, toby@telegraphics.com.au
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by  
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License  
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

# Solutions to problems described at:
# http://projecteuler.net/index.php?section=problems

# return the sum of 1..n
procedure sum_series(n)
	return n*(n+1)/2  # Gauss' formula
end

# Add all the natural numbers below one thousand that are multiples of 3 or 5.
procedure problem1a(n)
	# dumb iterative method
	local s,sum
	s := set()
	every insert(s, 3*(1 to (n-1)/3) | 5*(1 to (n-1)/5))
	sum := 0
	every sum +:= !s
	write(sum)
end

procedure problem1(n)
	# use formula for sum of series
	n -:= 1
	write(3*sum_series(n/3) + 5*sum_series(n/5) - 15*sum_series(n/15))
end

# generate the Fibonacci numbers not greater than n
procedure fib(n)
	local a,b
	a := b := 1
	while a <= n do {
		suspend a
		a +:= b
		b :=: a
	}
end

# Find the sum of all the even-valued terms in the Fibonacci sequence which do not exceed four million.
procedure problem2(n)
	# note that every third Fibonacci number is even, but we don't use that fact
	# there is also a formula for sum of n Fibonacci terms, but we don't use that either :)
	local sum,q
	sum := 0
	every q := fib(n) & q%2 = 0 & sum +:= q
	write(sum)
end

# generate prime factors from largest to smallest
# I don't think this is very efficient
procedure prime_factor(n)
	local j
	j := 2
	while j <= n do
		if n%j = 0 then {
			suspend prime_factor(n/j)
			return j
		}
		else j +:= 1
end

# Find the largest prime factor of a composite number.
procedure problem3(n)
	write(prime_factor(n))
end

# Find the largest palindrome made from the product of two 3-digit numbers.
procedure problem4()
	# brute force, check every combination
	local max,n,a,b
	max := n := 0
	every a := 100 to 999 & b := 100 to a do {
		prod := string(a*b)
		if prod == reverse(prod) then {
			n +:= 1
			max := max < a*b
		}
	}
	write("found ",n," palindromes, largest is ",max)
end

# What is the smallest number divisible by each of the numbers 1 to 20?
procedure problem5(n)
	# a.k.a. Least Common Multiple, check @ http://www.research.att.com/~njas/sequences/A003418
	local t,i,ft,k,lcm

	t := table(0)
	every i := 2 to n do {
		#write(i)
		ft := table(0)
		# find exponent of each prime factor
		every ft[ prime_factor(i) ] +:= 1
		# if exponent is higher than previously known, keep it
		every k := key(ft) do {
			#write("  ",k," ^ ",ft[k])
			t[k] := t[k] < ft[k];
		}
	}
	# multiply all factors
	#write("SOLUTION:")
	lcm := 1
	every k := key(t) do {
		#write(k," ^ ",t[k])
		lcm *:= k ^ t[k]
	}
	write(lcm)

	every write("ERROR! ",i," ",0 ~= lcm % (i := 1 to n))  # test each divisor
end

# What is the difference between the sum of the squares and the square of the sums?
procedure problem6(n)
	# expand square of sums (a+b+c+d)(a+b+c+d)
	# = sums of squares + 2*( sum of products of each permutation of terms )
	# since we are subtracting the squares, we can ignore that part entirely
	# just compute the sum of products
	# this becomes:
	#   2 * (
	#     1 x (2 + ... + n) +
	#     2 x (3 + ... + n) + ...
	#     (n-1) x n
	#   )
	# no doubt this can be analysed further, but I just decided to stop
	# and iterate over each series, using the formula.

	local i,j,sum

	i := sum_series(n)
	sum := 0
	every j := 1 to n-1 do
		sum +:= j*(i - sum_series(j))
	write(2*sum)
end

# What is the 10001st prime number?
procedure problem7a(q)
	local p,primes
	# inefficient method
	p := 2
	primes := list()
	until *primes = q do {
		# if p is not prime (not divisible by any previous prime), then we must ignore it
		if not ( q := !primes & p%q = 0 ) then {
			# we have a new prime!
			put(primes,p)
		}
		p +:= 1
	}
	write(pull(primes))
end

# using Eratosthenes' sieve, return the set of primes not greater than n
procedure set_primes(n)
	local s,i
	s := set()
	every insert(s, 2 to n)
	every i := 2 to sqrt(n) & member(s, i) & every delete(s, i*(j := 2 to n/i))
	return s
end

procedure problem7(q)
	write(sort(set_primes(105000))[q])  # value greater than 10001'st prime (I cheated)
end

procedure problem7b(q)
	# direct translation of e.e.coli's Ruby solution
	# slower than sieve but requires no estimate of value
	local n,primes,max,p
	n := 3
	primes := [2]
	while *primes < q do {
		max := sqrt(n)
		every p := !primes do {
			if p > max then {
				put(primes, n)
				break
			}
			if n % p = 0 then break
		}
		n +:= 2  # slight improvement on e.e.coli's
	}
	write(pull(primes))
end

procedure gen_digits(n)
	while n > 0 do {
		suspend n % 10
		n /:= 10
	}
end

# Discover the largest product of five consecutive digits in the 1000-digit number.
procedure problem8(n)
	local max,prod,digits
	max := prod := 0
	digits := list()
	every put(digits, gen_digits(n)) do {
		if *digits > 5 then pop(digits)
		prod := 1
		every prod *:= !digits
		max := max < prod
	}
	write(max)
end

procedure mul_digits(s)
	local prod
	prod := 1
	every prod *:= integer(!s)
	return prod
end

procedure problem8b(n)
	local max,s
	# alternative method, using string
	max := 0
	s := string(n)
	every max := max < mul_digits(s[ (i := 1 to *s-4) +: 5 ])
	write(max)
end

# Find the only Pythagorean triplet, {a, b, c}, for which a + b + c = 1000.
procedure problem9(t)
	local a,b,c
	# a^2 + b^2 = c^2
	# a + b + c = 1000
	# a < b < c
	# brute force solution
	until b := 2 to t/2 & a := 1 to b-1 & c := t-a-b & a*a + b*b = c*c
	write(a*b*c)
end

# Calculate the sum of all the primes below two million.
procedure problem10(n)
	local sum
	sum := 0
	every sum +:= !set_primes(n-1)
	write(sum)
end

procedure gen_4(t)
	every i := 1 to *t[1]-3 & j := 1 to *t-3 do {
		suspend t[j][i +: 4]  # horizontal
		suspend [ t[j][i],   t[j+1][i],   t[j+2][i],   t[j+3][i]   ]  # vertical
		suspend [ t[j][i],   t[j+1][i+1], t[j+2][i+2], t[j+3][i+3] ]  # diagonal down-right
		suspend [ t[j][i+3], t[j+1][i+2], t[j+2][i+1], t[j+3][i]   ]  # diagonal down-left
	}
end

# What is the greatest product of four numbers on the same straight line in the 20 by 20 grid?
procedure problem11()
	local max,prod,t
	t := [[08,02,22,97,38,15,00,40,00,75,04,05,07,78,52,12,50,77,91,08],
		  [49,49,99,40,17,81,18,57,60,87,17,40,98,43,69,48,04,56,62,00],
		  [81,49,31,73,55,79,14,29,93,71,40,67,53,88,30,03,49,13,36,65],
		  [52,70,95,23,04,60,11,42,69,24,68,56,01,32,56,71,37,02,36,91],
		  [22,31,16,71,51,67,63,89,41,92,36,54,22,40,40,28,66,33,13,80],
		  [24,47,32,60,99,03,45,02,44,75,33,53,78,36,84,20,35,17,12,50],
		  [32,98,81,28,64,23,67,10,26,38,40,67,59,54,70,66,18,38,64,70],
		  [67,26,20,68,02,62,12,20,95,63,94,39,63,08,40,91,66,49,94,21],
		  [24,55,58,05,66,73,99,26,97,17,78,78,96,83,14,88,34,89,63,72],
		  [21,36,23,09,75,00,76,44,20,45,35,14,00,61,33,97,34,31,33,95],
		  [78,17,53,28,22,75,31,67,15,94,03,80,04,62,16,14,09,53,56,92],
		  [16,39,05,42,96,35,31,47,55,58,88,24,00,17,54,24,36,29,85,57],
		  [86,56,00,48,35,71,89,07,05,44,44,37,44,60,21,58,51,54,17,58],
		  [19,80,81,68,05,94,47,69,28,73,92,13,86,52,17,77,04,89,55,40],
		  [04,52,08,83,97,35,99,16,07,97,57,32,16,26,26,79,33,27,98,66],
		  [88,36,68,87,57,62,20,72,03,46,33,67,46,55,12,32,63,93,53,69],
		  [04,42,16,73,38,25,39,11,24,94,72,18,08,46,29,32,40,62,76,36],
		  [20,69,36,41,72,30,23,88,34,62,99,69,82,67,59,85,74,04,36,16],
		  [20,73,35,29,78,31,90,01,74,31,49,71,48,86,81,16,23,57,05,54],
		  [01,70,54,71,83,51,54,69,16,92,33,48,61,43,52,01,89,19,67,48]]
	max := 0
	every f := gen_4(t) do {
		prod := 1
		every prod *:= !f
		max := max < prod
	}
	write(max)
end

procedure count_factors(k)
	local d,f
	d := 0
	every f := 1 to sqrt(k) & k % f = 0 & d +:= 2
	return if f = sqrt(k) then d+1 else d
end

# generate all factors
procedure gen_factors(k)
	local f,t
	t := 0
	every f := 1 to sqrt(k) & k % f = 0 do suspend t := k/f | t ~= f
end

# What is the value of the first triangle number to have over five hundred divisors?
procedure problem12(n)
	local i,tri
	i := 0
	until count_factors(tri := sum_series(i +:= 1)) > n
	write(tri)
end

procedure main()
	problem1(1000)
	#problem2(4000000)
	#problem3(600851475143)
	#problem4()
	#problem5(20)
	#problem6(100)
	#problem7(10001)
	#problem8(7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450)
	#problem9(1000)
	#problem10(2000000)
	#problem11()
	#problem12(500)
end