[ create a new paste ] login | about

Link: http://codepad.org/VYXqTeQh    [ raw code | fork ]

Haskell, pasted on Oct 28:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
{-# LANGUAGE BangPatterns #-}

module ProjectEuler (
	euler
) where


import Data.List
import Data.Maybe
import Data.Char
import Data.Ord (comparing)
import Data.Ratio
import Control.Monad.State
import Control.Parallel.Strategies
import qualified Text.Parsec as Parsec
import qualified Debug.Trace as Debug
import qualified Data.Map as Map
import Data.Word
import qualified Data.Numbers.Primes as Primes -- cabal install primes
import qualified Test.QuickCheck as QC
import qualified Data.Permute as Permute
import qualified Data.Vector as Vec -- cabal install vector
import qualified Data.Array.IArray as A

euler :: Int -> IO Integer

-- All performance is measured in GHCi, i.e. without optimization.

-- Definitions used in multiple problems

-- Ceiled integer square root
iSqrt :: (Integral a) => a -> a
iSqrt = floor . sqrt . fromIntegral

-- Divisibility test
a `divisibleBy` b = a `mod` b == 0

-- List of all divisors, i.e. proper divisors plus the number itself
divisors n = map head . group . sort $ smaller ++ greater
	where
		smaller = filter (n `divisibleBy`) [1..iSqrt n]
		greater = map (n `quot`) smaller

-- List of proper divisors of n
properDivisors n
	| n == 1 = [1]
	| otherwise = init . divisors $ n


digitSum base = sum . reverseExplodeInt base

-- explodeInt with the result reversed. This is the natural outcome of the
-- unfoldr, and it is useful if the order of the digits obtained doens't matter,
-- e.g. for calculating the digit sum.
reverseExplodeInt base = unfoldr extractLastDigit
	where
		extractLastDigit n'
			| n' <= 0   = Nothing
			| otherwise = Just (n' `rem` base, n' `quot` base)

-- Explodes an int to a list of digits.
-- explodeInt 10 1234 -> [4,3,2,1]
explodeInt base = reverse . reverseExplodeInt base

-- Inverse of explodeInt.
-- implodeInt 10 [1,2,3,4] -> 1234
implodeInt base ns = sum $ zipWith (*) (reverse ns) exponents
	where
		exponents = map (base^) [0..]


-- Fibonacci numbers
fibo = 1 : 1 : zipWith (+) fibo (tail fibo)

-- | Faculty function
faculty n = foldl' (*) 1 [1..n] -- Using foldl1 crashes for n=0


-- Splits a list at certain elements.
-- explode ',' "1,23,456" -> ["1", "23", "456"]
explode divider s =
	let
		(l, s') = break (== divider) s
	in
		l : case s' of
			[]      -> []
			(_:s'') -> explode divider s''


-- | Binomial coefficient. Error for n or k < 0, and for n < k.
choose :: (Integral a) => a -> a -> a
choose n k
	| n < k     = error "choose n k with n < k"
	| n < 0     = error "choose n k with n < 0"
	| k < 0     = error "choose n k with k < 0"
	| k == 0    = 1
	| n == 0    = 0
	| n - k < k = choose n (n - k)
	| otherwise = (n - k + 1) * choose n (k-1) `quot` k



{-
	Problem 1
		Find the sum of all the multiples of 3 or 5 below 1000.

	Result
		.01 s

	Comment
		The formula below is based on sum [1..n] == n(n+1)/2.
-}
euler 1 = return $ (3 * n3 * (n3 + 1) + 5 * n5 * (n5 + 1)) `div` 2
	where
		n3 = 999 `div` 3
		n5 = 999 `div` 5
-- Old implementation
-- euler 1 = return $ sum [0,3..999] + sum [0,5..999]
-- Even older implementation
-- euler 1 = return $ sum [n | n <- [1..999], n `divisibleBy` 5 || n `divisibleBy` 3]



{-
	Problem 2
		Find the sum of the even-valued Fibonacci numbers <= 4e6.

	Performance:
		.00 s
-}
euler 2 = return . sum . filter even $ takeWhile (<= 4000000) fibo



{-
	Problem 3
		The prime factors of 13195 are 5, 7, 13 and 29.
		What is the largest prime factor of the number 600851475143?

	Performance:
		.01 s
-}
euler 3 = return . last . Primes.primeFactors $ 600851475143


{-
	Problem 4
		A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99.
		Find the largest palindrome made from the product of two 3-digit numbers.

	Performance:
		.67 s
-}
euler 4 = return $ maximum [i*j | i <- [1..999], j <- [1..i], isPalindromic (i*j)]
	where
		isPalindromic n =
			let
				showNumber = show n
			in
				showNumber == reverse showNumber




{-
	Problem 5
		What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

	Result
		.00 s

	Comment
		Haskell provides a lcm function in the Prelude, but I re-implemented
		it so that the problem is less of a stub.
-}
euler 5 = return $ foldl1' lcm' [1..20]
	where
		lcm' a b = (a * b) `div` (gcd' a b)
		gcd' a b
			| b == 0    = a
			| otherwise = gcd' b (a `mod` b)



{-
	Problem 6
		Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

	Result
		.00 s

	Alternative solution:
		The formula can be simplified, yielding
		return $ (n-1) * n * (1+n) * (2+3*n) `div` 12

-}
euler 6 = return $ squareSum - sumSquares
	where
		squareSum  = ((^2) . sum) range
		sumSquares = (sum . map (^2)) range
		range = [1..100]



{-
	Problem 7
		What is the 10001st prime number?

	Performance:
		.04 s
-}
euler 7 = return $ Primes.primes !! (10001-1) -- Lists start counting at 0



{-
	Problem 8
		Find the greatest product of five consecutive digits in the 1000-digit number.

	Performance:
		.03 s
-}
euler 8 = fmap (largestProduct 5 . explodeInt 10) largeNumber
	where
		-- Finds the largest product of `len` consecutive list elements in a list
		largestProduct _ [] = 1
		largestProduct len xx@(x:xs)
			| length xx <= len = product xx
			| otherwise = max (product $ take len xs) (largestProduct len xs)

		largeNumber = fmap read $ readFile "problem_8" :: IO Integer



{-
	Problem 9
		Find the only Pythagorean triplet, {a, b, c}, for which a + b + c = 1000.

	Result
		.58 s
-}
euler 9 = return . head $ do
	c <- [1..1000]
	b <- [1..c-1]
	guard $ (1000-b-c)^2 + b^2 == c^2
	guard $ b + c < 1000
	guard $ (1000-b-c) <= b
	return $ (1000-b-c)*b*c



{-
	Problem 10
		Find the sum of all primes <= 2 * 10^6

	Result
		139 s
-}
euler 10 = return . sum . takeWhile (<= 2 * 10^6) $ Primes.primes



{-
	Problem 12
		What is the value of the first triangle number to have over five hundred divisors?

	Result
		0.74 s
-}
euler 12 = return $ head [triangle n | n <- [1..], triangleD n > 500]
	where
		-- Calculates the n-th triangle number
		triangle n = n * (n+1) `div` 2

		-- Number of divisors of the triangle number n*(n+1)/2. Uses the
		-- fact that n and (n+1) are coprime, so the divisor count is a
		-- homomorphism on them.
		triangleD n
			| even n = d (n `div` 2) * d (n + 1) -- Make sure the even value
			|  odd n = d n * d ((n + 1) `div` 2) -- is divided by 2
		d = numDivisors

		-- Number of divisors.
		-- Based on the fact that a number with prime factors of multi-
		-- plicities a1, a2, ... has (a1+1)(a2+1)... divisors.
		-- Equivalent to length.divisors, but performs way better.
		numDivisors = product . map (+1) . pfm
		-- pfm = prime factor multiplicities
		pfm = map length . group . Primes.primeFactors



{-
	Problem 13
		Work out the first ten digits of the sum of one-hundred 50-digit numbers.

	Result
		.01 s
-}
euler 13 = fmap (first10 . sum) numbers
	where
		numbers = fmap (map read . explode '\n') $ readFile "problem_13"
		first10 n
			| n < 10^10 = n -- 10^10 is the first number with 11 digits
			| otherwise  = first10 $ n `div` 10


{-
	Problem 14
		Which starting number, under one million, produces the longest Collatz chain?

	Result
		10 s
-}
euler 14 = return $ fst maxTuple
	where
		maxN :: Integer
		maxN = 10^6-1

		-- How many numbers to memoize.
		-- Too large consumes too much space, too little and memoization
		-- is useless.
		memoN :: Integer
		memoN = maxN

		-- Vector of the type [(n, collatz steps n has to take to reach loop)],
		-- intended to memoize already calculated lengths below a million.
		-- The first element is (1,1), since 1 has to take 1 step to loop.
		collatzLengthList :: A.Array Integer Integer
		collatzLengthList = A.array (1, memoN) [(n, collatzLength collatzLengthList n 1) | n <- [1..memoN]]
		-- collatzLengthList = V.generate (maxN+1) (\n -> if n == 0 then (0,-1) else (n, collatzLength collatzLengthList n 1))

		collatzLength :: A.Array Integer Integer -> Integer -> Integer -> Integer
		collatzLength memo n c
			| n == 1                                           = 1
			| even n && A.bounds memo `A.inRange` (n `quot` 2) = c + (memo A.! (n `quot` 2))
			| odd  n && A.bounds memo `A.inRange`  (3 * n + 1) = c + (memo A.! (3 * n +  1))
			| even n                                           = collatzLength memo (n `quot` 2) (c+1)
			| otherwise                                        = collatzLength memo (3*n + 1) (c+1)

		maxTuple = foldl1 maxC . take (fromIntegral maxN) $ A.assocs collatzLengthList
			where
				maxC a@(!a1,!a2) x@(_,x2)
					| a2 >= x2  = a
					| otherwise = x

-- Old algorithm, using no memoization. Performance: 228 s (GHCi), 1.3 s (GHC).
--
-- euler 14 = return . fromIntegral . snd . maxTuple $ [1..10^6-1]
-- 	where

-- 		collatzLength :: Word32 -> Word16
-- 		collatzLength n = go (1, n)
-- 			where
-- 				-- 'go (1, n)' calculates the length of the collatz chain
-- 				-- starting at n.
-- 				go :: (Word16, Word32) -> Word16
-- 				go (!c, !n) -- c = current count so far, n = current number
-- 					| n == 1    = c
-- 					| even n    = go (c + 1, n `div` 2)
-- 					| otherwise = go (c + 1, 3 * n + 1)

-- 		-- 'collatzMax prev candid' returns the tuple with the longer chain,
-- 		-- which is either the unmodified old tuple, or the one generated
-- 		-- out of the starting value provided as the second argument.
-- 		collatzMax :: (Word16, Word32) -> Word32 -> (Word16, Word32)
-- 		collatzMax previous candidate =
-- 			max previous (collatzLength candidate, candidate)

-- 		-- The accumulator is of the form '(length, start_value)' of the
-- 		-- maximum entry so far.
-- 		maxTuple :: [Word32] -> (Word16, Word32)
-- 		maxTuple = foldl collatzMax (1, 1)



{-
	Problem 15
		How many routes are there through a 20*20 grid?

	Result
		Perfect

	Comment
		To traverse the 20*20 grid, one has to take 20 steps right and
		another 20 down. The problem thus reduces in "how many different
		ways are there to arrange right and down", which is simply the
		binomial coefficient.
-}
euler 15 = return $ (20+20) `choose` 20



{-
	Problem 16
		Digit sum of 2^1000?

	Result
		Perfect
-}
euler 16 = return . digitSum 10 $ 2^1000



{-
	Problem 18
		Maximum sum in a triangle of numbers

	Solution idea
		Reduce the triangle row by row, starting from the bottom. For
		example, the bottom right corner is the triangle [[31],[04,23]]. If
		one reaches the 31, the maximizing step is taking 23, 04 is not an
		option. Therefore, the cell with 31 can effectively be considered a
		cell on the bottom of the triangle with value 31+23=54. Doing this
		for all 3-digit triangles in the bottom row reduces the triangle by
		one row; iterating the process until the tip of the pyramid has been
		reached yields the desired result.

	Result
		.01 s
-}
euler 18 = do
		triangle <- triangleData
		return . head $ foldr1 rowReduce triangle
	where

		-- rowReduce compares all the small triangles in the bottom rows,
		-- and reduces the two lines to one effective line. Folding this
		-- over the pyramid yields a singleton list containing the result.
		-- Example: rowReduce [1,2] [7,1,4] -> [8,6]
		rowReduce' upper lower = zipWith max lShifted rShifted
			where
				lShifted = zipWith (+) upper (init lower) -- init unnecessary, but makes the code nicely symmetric :-)
				rShifted = zipWith (+) upper (tail lower)
		-- Golfed version:
		rowReduce u l=z max (p id)$p tail
			where
				p f=z (+) u$f l
				z=zipWith

		-- Parsec
		triangleData = fmap parseNumberFile $ readFile "problem_18"
		numberFile = Parsec.sepBy line eol
		line = Parsec.sepBy cell $ Parsec.char ' '
		cell = fmap read $ Parsec.many Parsec.digit
		eol = Parsec.char '\n'
		parseNumberFile = unEither . Parsec.parse numberFile ""
		unEither = either (error . show) id

{-
	Problem 20
		Digit sum of 100!?

	Result
		0.00 s
-}
euler 20 = return . digitSum 10 $ faculty 100


{-
	Problem 21
		Find the sum of all the amicable numbers under 10000.
		(Amicable == d^2 = id; d = sum of proper divisors

	Solution
		31626

	Result
		1.84 s
-}
euler 21 = return . sum $ filter isAmicable [1..maxN]
	where
		maxN = 10^4
		d = sum . properDivisors
		isAmicable n =
			let dn = d n
			in  dn <= maxN && dn /= n && d dn == n


{-
	Problem 22
		Process a file with a lot of names

	Result
		.09 s
-}
euler 22 = fmap nameListScore names
	where
		-- Calculates the total score of a list of names.
		-- The score consists of the position of the name in the sorted list
		-- (starting with 1) multiplied by the name score.
		-- Example: "ABD" at position 8 yields (1+2+4)*8 = 56.
		nameListScore :: (Enum a, Num a) => [[Char]] -> a
		nameListScore names = sum $ zipWith (*) [1..] (fmap nameScore $ sort names)

		-- Calculates the name score of a single name.
		-- Example: "ABD" = 1+2+4 = 7
		nameScore :: (Num a) => String -> a
		nameScore name = fromIntegral . sum $ map (\x -> ord x - ord 'A' + 1) name

		-- Reads the names file, creating data of the type
		-- IO [name1, name2, ...]
		names :: IO [String]
		names = fmap (explode ',' . filter (/= '"')) $ readFile "problem_22"




{-
	Problem 25
		Index of the first Fibonacci number with 1000 digits in the
		Fibonacci sequence?

	Result
		4782
		0.11 s
-}
euler 25 = return . fst . fromJust $ find (has1000 . snd) fiboPairs
	where
		fiboPairs = zip [1..] fibo -- Pairs of (n, F_n)
		has1000 n = n >= 10^999    -- Does n have 1000 digits?



{-
	Problem 28
		Find the sum of corners in a spiral list

	Result
		669171001
		0.00 s
-}
euler 28 = return . sum $ spiralList
	where
		-- List of offsets to the next number.
		-- = [1, 2,2,2,2, 4,4,4,4, 6,6,6,6, ...]
		spiralOffsets = 1 : ([2,4..1000] >>= replicate 4)

		-- Accumulate offsets to generate the actual list of spiral elements
		spiralList = scanl1 (+) spiralOffsets


{-
	Problem 29
		How many distinct values of a^b for 2 <= a,b <= 100 are there?

	Result
		.94 s
-}
euler 29 = return . fromIntegral . numDistinct $ powers
	where
		numDistinct = length . group . sort
		powers = liftM2 (^) [2..100] [2..100]


{-
	Problem 30
		Numbers that can be written as powers of their digits

	Result
		6.28 s

	Comment
		The upper boundary can be estimated since 999... = 10^k - 1 has to
		be equal to 9^5 + 9^5 + ... = k 9^5, which yields the maximum
		condition k 9^5 = 10^k - 1. A numeric solution for this is 5.51257,
		which yields a maximum of 10^5.51257 = 325514.24.
-}
euler 30 = return . fromIntegral . sum $ filter f [2..325515]
	where
		-- f x = can x be written as the sum of fifth power of its digits?
		f x = x == (sum . map toTheFifth $ explodeInt 10 x)

		-- Memoize powers
		toTheFifth = Vec.unsafeIndex fifthPowers
		fifthPowers = Vec.generate 10 (^5)


{-
	Problem 34
		Numbers that can be written as factorials of their digits

	Result
		54 s

	Comment
		Similar estimate as in problem 30 results in a numerical upper
		boundary estimate of 2309192.
-}
euler 34 = return . fromIntegral . sum $ filter isSumOfFacDigits [3..2309192]
	where
		-- f x = can x be written as the sum of faculties of its digits?
		isSumOfFacDigits x = x == (sum . map faculty' $ explodeInt 10 x)

		-- Hashed faculty function
		faculty' = Vec.unsafeIndex faculties
		faculties = Vec.generate 10 faculty



{-
	Problem 35
		How many circular primes are there below a million?

	Result
		87 s
-}
euler 35 = return . genericLength $ filter isCircularPrime candidates
	where
		candidates = takeWhile (< 10^6) Primes.primes
		isCircularPrime p = all Primes.isPrime rotations
			where
				digits = explodeInt 10 p
				rotateList (x:xs) = xs ++ [x]
				rotations = map (implodeInt 10) .
				            take (length digits) $
				            iterate rotateList digits



{-
	Problem 36
		Find the sum of all numbers below 10^6 that are palindromic in
		base 2 and 10.

	Result
		18 s
-}
euler 36 = return . sum $ filter isDoublePali [1..10^6]
	where
		isDoublePali n = isBasePali 10 n && isBasePali 2 n
		isBasePali base number =
			let
				exploded = reverseExplodeInt base number
				-- reverseExplodeInt is slightly more effective than
				-- explodeInt
			in
				exploded == reverse exploded



{-
	Problem 37
		Find all 11 primes that are left- and right-truncatable

	Result
		2.7 s
-}
euler 37 = return . sum . take 11 . filter isTruncatable $ candidates
	where
		candidates = dropWhile (<= 7) Primes.primes -- exclude single digit primes
		isTruncatable n = allPrimes truncates
			where
				n' = explodeInt 10 n
				allPrimes = and . map (Primes.isPrime . implodeInt 10)
				truncates = (middle $ inits n') ++ (middle $ tails n')
				middle = init . tail -- Everything but head and last
				                     -- (since inits xs = [xs ... []])



{-
	Problem 40
		Product of decimals in a long number

	Result
		2.2 s
-}
euler 40 = return . product $ map (concatInts !!) d_n
	where
		-- Concatenated decimals
		-- = [1,2,3,4,5,6,7,8,9,1,0,1,1,1,2,1,3,1,4,1,5,...]
		concatInts = [1..] >>= explodeInt 10

		-- Indices of d to be considered
		-- = [1, 10, 100, 1000, 10000, 100000, 1000000]
		d_n = [10^n - 1 | n <- [0..6]] -- -1 because lists are 0-indexed



{-
	Problem 41
		Find the largest pandigital prime

	Result
		0.01 s

	Comment
		"Why does the search start with a 7-pandigital and not one of rank 9?
		All n-pandigitals of a certain order have the same digit sum:
			n = 9 => digit sum 45
			n = 8 => digit sum 36
			n = 7 => digit sum 28
			...
		This shows that n = 8,9 cannot be pandigital primes, as their digit
		sum is divisible by 3. It is also known that there is a pandigital
		prime because the problem gives 2143 as an example, which is a
		4-pandigital. The rest of the algorithm is brute force.
-}
euler 41 = return . fromIntegral . prevPanPrime $ Permute.listPermute 7 [6,5..0]
	where

		prevPanPrime p
			| Primes.isPrime n = n
			| otherwise = prevPanPrime . maybe prunedP id $ Permute.prev p
			where
				n = getNumber p
				-- When the algorithm reaches the smallest n-pandigital,
				-- this generates the permutation for the largest
				-- (n-1) pandigital. Example: p = 123 => prunedP = 21
				-- Note that this is not actually needed here, since
				-- there is a 7-digit pandigital, but you can't know
				-- that in advance :-)
				prunedP = Permute.listPermute (length ipe) $ reverse ipe
				ipe = init . Permute.elems $ p

				-- Extract current number out of a permutation
				getNumber = implodeInt 10 . map (+1) . Permute.elems



{-
	Problem 48
		Find the last ten digits of the series 1^1 + 2^2 + 3^3 + ... + 1000^1000

	Result
		.03 s
-}
euler 48 = return $ sum [n^n | n <- [1..1000]] `mod` 10^10



{-
	Problem 53
		How many binomial coefficients with 1 <= n <= 100 exceed 10^6?

	Result
		.42 s
-}
euler 53 = return $ genericLength [() | n <- [1..100],
                                        k <- [1..n],
                                        choose n k > 10^6
                                        ]



{-
	Problem 56
		Considering natural numbers of the form a^b, where a, b < 100, what
		is the maximum digital sum?

	Result
		1.8 s
-}
euler 56 = return . maximum . map (digitSum 10) $
           [a^b | a <- [1..100], b <- [1..100]]


-- Code copied from problem 18. Do not modify the code here. Performance: .05 s
euler 67 = do
		triangle <- triangleData
		return . head $ foldr1 rowReduce triangle
	where
		rowReduce upper lower = zipWith max lShifted rShifted
			where
				lShifted = zipWith (+) upper (init lower)
				rShifted = zipWith (+) upper (tail lower)
		triangleData = do
			contents <- readFile "problem_67"
			return . map (map (read :: String -> Integer)) . unRight $ parseNumberFile contents
			where
				numberFile = Parsec.endBy line eol
				line = Parsec.sepBy cell (Parsec.char ' ')
				cell = Parsec.many (Parsec.noneOf " \n")
				eol = Parsec.char '\n'
				parseNumberFile = Parsec.parse numberFile "(numberfile)"
				unRight (Left err)  = error $ show err
				unRight (Right x)   = x



{-
	Problem 71
		Find the numerator of the fraction before 3/7 in the Farey sequence
		for d <= 10^6.

	Result
		1.05 s

	Comment
		This algorithm was constructed using the small note from [1], which
		reads "For a method of computing a successive sequence from an
		existing one of n terms, insert the mediant fraction (a+b)/(c+d)
		between terms a/c and b/d when c+d<=n".

		This is *very* brief, so here's an elaboration on that:

		The mediant of two fractions a/b and c/d is (a+b)/(c+d). This
		fraction has two properties that are important here:
			1. a/b < mediant < c/d
			2. The mediant's denominator is the smallest denominator of
			   all the (reduced) fractions that lie between a/b and c/d.
		We're now looking for the fraction before 3/7 in the Farey Sequence.
		Calculating the mediant of these two has two possible outcomes:
			x) The mediant has a denominator greater than n, the index
			   of the Farey sequence. Since we know from 1. that it's
			   between the two fractions we put in and from 2. that it's
			   the one with the smallest possible denominator, we can
			   conclude that this mediant is not part of the Farey
			   Sequence. Therefore, the two input fractions are adjacent.
			y) The mediant is in the Farey Sequence and the denominator
			   is not too large. From 1. we can be sure that it's
			   strictly between the two input fractions. If we now change
			   the first input fraction to the generated mediant and
			   recurse, the algorithm is called again with a narrower
			   focus. Repeating this will eventually yield case x), and
			   the algorithm terminates.
		So to calculate the fraction before 3/7, the algorithm moves
		the left boundary towards 3/7, until it detects that there's no
		fraction between the boundaries, in which case two neighbours have
		been found.
-}
euler 71 = return . numerator $ fareyBefore (10^6) 0 (3/7)
	where fareyBefore n ac bd
		| c + d > n = ac
		| otherwise = fareyBefore n abcd bd
		where
			(a,c) = (numerator ac, denominator ac)
			(b,d) = (numerator bd, denominator bd)
			abcd = (a+b) % (c+d)

{-
Here's an algorithm to calculate the full Farey sequence F_n, deduced from
http://en.wikipedia.org/wiki/Farey_sequence#Next_term

fareySequence :: (Integral a) => a -> [Ratio a]
fareySequence n = 0 : unfoldr nextFarey (0 % 1, 1 % n)
	where
		-- Computes the next fraction out of the predecessors
		-- a%b and c%d.
		-- Reference: http://en.wikipedia.org/wiki/Farey_sequence#Next_term
		nextFarey (ab, cd)
			| ab == 1 = Nothing
			| otherwise = Just (cd, (cd, pq))
			where
				(a,b) = (numerator ab, denominator ab)
				(c,d) = (numerator cd, denominator cd)
				pq = p % q
				k = (n + b) `div` d
				p = k * c - a
				q = k * d - b
-}



{-
	Problem 73
		How many reduced fractions are between 1/3 and 1/2 with maximum
		denominator 12000?

	Result
		108 s

	Comment
		See problem 71 for more details on this, namely Farey Sequences.
-}
euler 73 = return $ fareyCount 12000 (1/3) (1/2)
	where
		mediant ab cd = (numerator ab + numerator cd)
		                %
		                (denominator ab + denominator cd)
		fareyCount n ab cd
			| denominator m > n = 0
			| otherwise = 1 + fareyCount n ab m + fareyCount n m cd
			where
				m = mediant ab cd



{-
	Problem 99
		Comparing large numbers of the form base^exponent

	Result
		0.04 s
-}
euler 99 = fmap (fst . maximumBy logCompare . zip [1..]) numbers
	where

		logCompare (_,[x,xExp]) (_,[y,yExp]) =
			compare
				(fromIntegral xExp * log (fromIntegral x))
				(fromIntegral yExp * log (fromIntegral y))

		-- Shorter comparison function, using Ord.comparing
		-- logCompare' = comparing (\(_, [x, xExp]) -> fromIntegral xExp * log (fromIntegral x))

		-- Parsec part (reading the numbers)
		numbers = do
			contents <- readFile "problem_99"
			let parsed = either (error.show) id $ parseNumberFile contents
			return . (map . map) (read :: String -> Integer) $ parsed

		numberFile = Parsec.endBy line eol
		line = Parsec.sepBy cell (Parsec.char ',')
		cell = Parsec.many (Parsec.noneOf ",\n")
		eol = Parsec.char '\n'
		parseNumberFile = Parsec.parse numberFile "(numberfile)"

euler n = error ("Problem " ++ show n ++ " has not yet been solved in Haskell.")



-- == PLAYGROUND ==


problem_26 = fst $ maximumBy (comparing snd)
                            [(n,recurringCycle n) | n <- [1..999]]
    where  recurringCycle d = remainders d 10 []
           remainders d 0 rs = 0
           remainders d r rs = let r' = r `mod` d
                               in case elemIndex r' rs of
                                    Just i  -> i + 1
                                    Nothing -> remainders d (10*r') (r':rs)



-- Recurring digit algorithm, raw and very buggy form
-- For example, 0 recurring digits should be in a Maybe type because there aren't always some
-- :: (digits after comma, length of the repeating decimals)
recurDiv n d rems
	| r == 0 = (q, 0) -- Division finishes
	| r `elem` rems = (q, ei) -- Recurring cycle finishes
	| otherwise = recurDiv (10*n) d (r : rems)
	where
		(q,r) = n `quotRem` d
		ei = fromJust (elemIndex r rems) + 1 -- Nothing shouldn't appear because this is only called when r`elem`rems

getRecur n d = case recurDiv n d [] of
	-- ~(digits, nonrep, rep) -> splitAt nonrep $ show digits
	~(digits, rep) -> digits `quotRem` (10^rep)

recurringDivide n d = (beforeComma, nonrecurring, recurring)
	where
		(beforeComma, afterComma) = n `quotRem` d
		(nonrecurring, recurring) = getRecur afterComma d


Create a new paste based on this one


Comments: