ザネリは列車を見送った ブログという名の備忘録

今回は僕にとっては難度が高かった…。
結局 Problem 26 は解法のヒントを検索して何とか解き、
Problem 27 は自力で解けるには解けたが実行時間がクッソ遅いのをどうにもできなかった…。

Problem 25. 1000-digit Fibonacci number

The Fibonacci sequence is defined by the recurrence relation:

Fn = Fn-1 + Fn-2, where F1 = 1 and F2 = 1.

Hence the first 12 terms will be:
F1 = 1
F2 = 1
F3 = 2
F4 = 3
F5 = 5
F6 = 8
F7 = 13
F8 = 21
F9 = 34
F10 = 55
F11 = 89
F12 = 144
The 12th term, F12, is the first term to contain three digits.
What is the first term in the Fibonacci sequence to contain 1000 digits?
Problem 25. 1000-digit Fibonacci number

---

フィボナッチ数列は以下の漸化式で定義される:

Fn = Fn-1 + Fn-2, ただし F1 = 1, F2 = 1.

最初の12項は以下である.

• F1 = 1
• F2 = 1
• F3 = 2
• F4 = 3
• F5 = 5
• F6 = 8
• F7 = 13
• F8 = 21
• F9 = 34
• F10 = 55
• F11 = 89
• F12 = 144

12番目の項, F12が3桁になる最初の項である.
1000桁になる最初の項の番号を答えよ.
Problem 25 「1000桁のフィボナッチ数」

最初に思いついたのがこちら。
まず、1から始まる項番号とフィボナッチ数列の無限リストをzipし、
n桁になる最初の値を 10^(n-1) で求め、それより小さいフィボナッチ数を dropWhile で落としている。
dropWhileしたリストの先頭がn桁になる最初のフィボナッチ数とその項なので、fstで項を取る。
少しスマートにしようと試みたのがこちら。
似たようなことをやっているが、1から始まる無限リストのフィボナッチ数列を算出していき、
findIndex で最初にn桁になる数のインデックスを取る。
インデックスは0から始まるが、この問題では項1から数えているので+1した値を取り出している。

Problem 26. Reciprocal cycles

A unit fraction contains 1 in the numerator. The decimal representation of the unit fractions with denominators 2 to 10 are given:
1/2 = 0.5
1/3 = 0.(3)
1/4 = 0.25
1/5 = 0.2
1/6 = 0.1(6)
1/7 = 0.(142857)
1/8 = 0.125
1/9 = 0.(1)
1/10 = 0.1
Where 0.1(6) means 0.166666..., and has a 1-digit recurring cycle. It can be seen that 1/7 has a 6-digit recurring cycle.
Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its decimal fraction part.
Problem 26. Reciprocal cycles

---

単位分数とは分子が1の分数である. 分母が2から10の単位分数を10進数で表記すると次のようになる.
1/2 = 0.5
1/3 = 0.(3)
1/4 = 0.25
1/5 = 0.2
1/6 = 0.1(6)
1/7 = 0.(142857)
1/8 = 0.125
1/9 = 0.(1)
1/10 = 0.1
0.1(6)は 0.166666... という数字であり, 1桁の循環節を持つ. 1/7 の循環節は6桁ある.
d < 1000 なる 1/d の中で小数部の循環節が最も長くなるような d を求めよ.
Problem 26 「逆数の循環節 その1」

苦戦した…結局自力では解法までたどり着けなかった。
まず、最初に Data.Fixed モジュールを使用して、精度の高い小数を扱おうと試みた。

import Data.Fixed (Fixed, HasResolution, resolution, showFixed)
import Data.List (elemIndex, maximumBy)
import Data.Ord (comparing)

data E300 = E300
instance HasResolution E300 where
    resolution _ = 10 ^ 300
type BigNum = Fixed E300

--main = print $ showFixed True $ fst $ maximumBy (comparing snd) $ map (\d -> (d, recurCycle $ dStr d)) [1,2..999]

-- 単位分数を10進数で表した小数部分の文字列
dStr :: BigNum -> String
dStr = reverse . dropWhile (=='0') . reverse . tail . dropWhile (/='.') . show . snd . properFraction . (1/)

-- 循環節の長さを返す
recurCycle :: Eq a => [a] -> Maybe Int
recurCycle xs = recurCycle' 1 xs
  where
    recurCycle' n xs = recurCycle'' $ (splited n xs) >>= cycleLength
      where
        recurCycle'' Nothing | n <= (length xs - 1) = recurCycle' (n+1) xs    -- 先頭要素と同じ要素が複数ある循環節に対応
                             | length xs > 1        = recurCycle' 1 (tail xs) -- リストの次の要素から始まる循環節に対応
                             | otherwise            = Nothing
        recurCycle'' x = x

-- aの繰り返しがbと等しい場合、aの長さを循環節の長さとして返す
cycleLength :: Eq a => ([a], [a]) -> Maybe Int
cycleLength (a, b) | isCycle   = Just (length a)
                   | otherwise = Nothing
  where isCycle = (length a) <= (length b) && (take (length b) (cycle a)) == b

splited :: Eq a => Int -> [a] -> Maybe ([a], [a])
splited _ [] = Nothing
splited n xs = fmap (flip splitAt xs) $ sameIdx n
  where sameIdx n = fmap (+n) $ elemIndex (head xs) (drop n xs)

同じ数値の並びが2回以上出現することで循環節とみなしているため、小数点以下の桁は循環節の長さ×2ないといけない。
しかし、resolution _ = 10 ^ 300 の値を増やし桁数を延ばして行くと Exception: Negative exponent が発生し、
十分な長さを取ることができなかった。
(また、どの程度が十分な長さかも分かっていなかった)
色々調べて、フェルマーの小定理から「循環節の長さは最大で n-1」ということを導けるらしく、それを踏まえた上で、このようにしてみた。
さらに、高速化のために数の大きい方から再帰的に調べていって、
循環節の長さの最大値が求められれば再帰を打ち切るようにしたのがこちら。
最初、 dPart 関数を dPart :: (Integral a, Show a) => a -> String と定義していたが、
それでは最終的な答えが7となってしまった。
型の違いによる桁の大きさの問題だろうか…
なぜ dPart :: Integer -> String として、toInteger d を渡さないといけないのか理解できていないのだが、
このようにすると euler26-1 と同値の回答がより速く得られた。

Problem 27. Quadratic primes

Euler discovered the remarkable quadratic formula:

n^2 + n + 41

It turns out that the formula will produce 40 primes for the consecutive values n = 0 to 39.
However, when n = 40, 402 + 40 + 41 = 40(40 + 1) + 41 is divisible by 41,
and certainly when n = 41, 412 + 41 + 41 is clearly divisible by 41.
The incredible formula n^2 - 79n + 1601 was discovered,
which produces 80 primes for the consecutive values n = 0 to 79.
The product of the coefficients, -79 and 1601, is -126479.
Considering quadratics of the form:
n^2 + an + b, where |a| < 1000 and |b| < 1000
where |n| is the modulus/absolute value of n
e.g. |11| = 11 and | -4| = 4
Find the product of the coefficients, a and b,
for the quadratic expression that produces the maximum number of primes for consecutive values of n, starting with n = 0.
Problem 27. Quadratic primes

---

オイラーは以下の二次式を考案している:

n^2 + n + 41.

この式は, n を0から39までの連続する整数としたときに40個の素数を生成する.
しかし, n = 40 のとき 402 + 40 + 41 = 40(40 + 1) + 41 となり41で割り切れる.
また, n = 41 のときは 412 + 41 + 41 であり明らかに41で割り切れる.

計算機を用いて, 二次式 n2 - 79n + 1601 という式が発見できた.
これは n = 0 から 79 の連続する整数で80個の素数を生成する.
係数の積は, -79 × 1601 で -126479である.

さて, |a| < 1000, |b| < 1000 として以下の二次式を考える (ここで |a| は絶対値): 例えば |11| = 11 |-4| = 4である.

n^2 + an + b

n = 0 から始めて連続する整数で素数を生成したときに最長の長さとなる上の二次式の, 係数 a, b の積を答えよ.
Problem 27 「二次式素数」

こちらはゴリ押しです…なんのひねりもなく、力技で愚直に問題の通りの式をコードに落とし込んで解いて、
最適化オプション付きコンパイルしても5分かかる…。
何かプログラミングテクニックか数学的知識の利用かで高速化できるとは思うのだが…うーむ。

色々小細工してこんな感じにもしてみた。
Nothing も含めて maybe 0 id するよりも少しはやることの節約になるかな、と思ったが、
コードがゴチャゴチャしただけで特に処理時間短縮にはならなかったな…。

import Data.List (findIndex, maximumBy)
import Data.Maybe (isJust)
import Data.Ord (comparing)
import Zaneli.Euler (isPrime)

main = let list = [(a*b, i) | a <- [-999..999], b <- [-999..999], let i = f a b, isJust i] in
       print $ fst $ maximumBy (comparing snd) list
  where f a b = findIndex (\n -> not $ isPrime (n^2 + a*n + b)) [0..]