この記事は ひとりアドベントカレンダーRosettaCodeで楽しむプログラミング Advent Calendar 2025の5日めの記事です。
英文タイトル長いのよ。
双子素数とは、$p$ と $p + 2$ がどちらも素数であるような $p$ と $p+2$ 両方のこと。
この問題では、双子を揃えて使う必要はないので、3,5 ゆえに双子素数である 3 と、5,7 ゆえに双子素数である 7 を使って 3 + 7 = 10 とするのもよし、双子素数である 5 を使って 5 + 5 = 10 とするのもよし、という条件。
この条件で、どうやっても作れない偶数を見つけよ、というのがお題。
タスク
- 5,000以下の 双子素数ふたつの和で表せない偶数 を全て挙げよ。
- 1を素数に加える（そして3は素数なので1は双子素数にも追加される）としたらどうなるか。
- （ボーナス）5,000超にはそのような数はないとされていることを、10,000まで、あるいはそれ以上好きなところまでの範囲で確認せよ。（少なくとも $10^9$ までは確認されている。）
考える
とりあえず双子素数の無限リストを作ろう。
import Data.Numbers.Primes
twinPrimes :: [Int]
twinPrimes =
[ q
| (p,q,r) <- zip3 (1:primes) primes (tail primes)
, p + 2 == q || q + 2 == r]
2列でやろうとするとプログラムがバタバタするので3列でやるのが正解。
ghci> take 20 twinPrimes
[3,5,7,11,13,17,19,29,31,41,43,59,61,71,73,101,103,107,109,137]
指定の上限までについて、双子素数2つの全ての組み合わせで足してみて、チェックの付かなかった偶数を取り出すのがタスク共通の仕事。
表には配列を用いるが、奇数の添え字の空間が無駄なので半分にして節約しよう。
import Data.Array.Unboxed
import Data.List
-- 上限ubまでの偶数で、双子素数リストpsを使って作れないものを見つける
engine :: Int -> [Int] -> [Int]
engine ub ps = map ((2 *) . fst) $ filter snd $ assocs checks
where
-- tailsを使い、足して ub を超えない組み合わせでの和を全て作る
pqss = takeWhile (not . null)
[ takeWhile (ub >=) $ map (p +) pqs
| pqs@(p:_) <- tails ps ]
-- pqssに出現しなかった偶数を見つける
checks :: UArray Int Bool
checks = accumArray (flip const) True (1, div ub 2)
[(div pq 2, False) | pqs <- pqss, pq <- pqs]
task1 = engine 5000 twinPrimes
task2 = engine 5000 (1 : twinPrimes)
task3 = dropWhile (5000 >) $ engine 10_000 twinPrimes
ghci> task1
[2,4,94,96,98,400,402,404,514,516,518,784,786,788,904,906,908,1114,1116,1118
,1144,1146,1148,1264,1266,1268,1354,1356,1358,3244,3246,3248,4204,4206,4208]
(0.04 secs, 15,478,656 bytes)
ghci> task2
[94,96,98,400,402,404,514,516,518,784,786,788,904,906,908,1114,1116,1118,1144
,1146,1148,1264,1266,1268,1354,1356,1358,3244,3246,3248,4204,4206,4208]
(0.03 secs, 10,730,200 bytes)
ghci> task3
[]
(0.05 secs, 32,684,624 bytes)
ghci> dropWhile (5000 >) $ engine 100_000 twinPrimes
[]
(2.07 secs, 895,809,368 bytes)
タスク3は上限を10倍にするたびにかかる時間は100倍になる。
OEIS
Rosetta CodeにはHaskell版はないが、OEISにはある。
-- A007534 Positive even numbers that are not the sum of a pair of twin primes.
import qualified Data.Set as Set (map, null)
import Data.Set (empty, insert, intersection)
a007534 n = a007534_list !! (n-1)
a007534_list = f [2, 4..] empty 1 a001097_list where
f xs'@(x:xs) s m ps'@(p:ps)
| x > m = f xs' (insert p s) p ps
| Set.null (s `intersection` Set.map (x -) s) = x : f xs s m ps'
| otherwise = f xs s m ps'
-- A001097 Twin primes.
a001097 n = a001097_list !! (n-1)
a001097_list = filter ((== 1) . a164292) [1..]
-- A164292 Binary sequence identifying the twin primes (characteristic function of twin primes: 1 if n is a twin prime else 0).
a164292 1 = 0
a164292 2 = 0
a164292 n = signum (a010051' n * (a010051' (n - 2) + a010051' (n + 2))) -- ダッシュ付きは意味不明
-- A010051 Characteristic function of primes: 1 if n is prime, else 0.
import Data.List (unfoldr)
a010051 :: Integer -> Int
a010051 n = a010051_list !! (fromInteger n-1)
a010051_list = unfoldr ch (1, a000040_list) where
ch (i, ps'@(p:ps)) = Just (fromEnum (i == p),
(i + 1, if i == p then ps else ps'))
-- A000040 The prime numbers.
import Data.List (genericIndex)
a000040 n = genericIndex a000040_list (n - 1)
a000040_list = base ++ larger where
base = [2, 3, 5, 7, 11, 13, 17]
larger = p : filter prime more
prime n = all ((> 0) . mod n) $ takeWhile (\x -> x*x <= n) larger
_ : p : more = roll $ makeWheels base
roll (Wheel n rs) = [n * k + r | k <- [0..], r <- rs]
makeWheels = foldl nextSize (Wheel 1 [1])
nextSize (Wheel size bs) p = Wheel (size * p)
[r | k <- [0..p-1], b <- bs, let r = size*k+b, mod r p > 0]
data Wheel = Wheel Integer [Integer]
書いたのは、まぎれもないヤツさ。
ghci> take 35 a007534_list
[2,4,94,96,98,400,402,404,514,516,518,784,786,788,904,906,908,1114,1116,1118
,1144,1146,1148,1264,1266,1268,1354,1356,1358,3244,3246,3248,4204,4206,4208]
(0.38 secs, 163,941,504 bytes)
何をしているのだろう？コードを整理して読み解く。
a007534_list = f [2, 4..] empty 1 twinPrimes
where
f xxs@(x:xs) s m pps@(p:ps)
| x > m = f xxs (insert p s) p ps
| cond = x : f xs s m pps
| otherwise = f xs s m pps
where
cond = Set.null (s `intersection` Set.map (x -) s)
xxs 偶数リスト 双子素数の足し合わせで作れるか判定されてリストに送出される情報源
s $S$ これから検討する $x$ 以下の双子素数全ての集合(Set Integer)。そうなるように最初の等式で維持する。
m s を管理するための補助変数、s の最大値（厳密には、x以上の最小の素数となる）
pps 双子素数リスト 使うタイミングになったら s に順に移される
cond 今回の $x$ に関して、$p \in S$ で $x - p \in S$ を満たすものがない、つまり、（$x$ 以下の）双子素数の足し合わせで表現できないとき真
真なら x は2つめの等式で出力され、さもなくば3つめの等式でスルーされる
個人的には xxs は最後の引数にして欲しいな、とかあるけど、まぁこんなものか。
cond も「共通部分が空集合」より「条件を満たす要素が存在しない」 $\forall p \in S; x - p \not\in S$
cond = all (\p -> notMember (x - p) s) (Set.elems s)
の方が判りやすくないかしら。
今日はここまで。