Help us understand the problem. What is going on with this article?

有理数とか有限体とかのはなし

More than 3 years have passed since last update.

一部間違いがあります、本文中でも明記しましたが最初の終了条件はInt を覆っていません。
最後に訂正したバージョンとより厳しいquickCheck の結果も載せておきました。

Haskell Advent Calendar 2016 の10日目です。

去年は眺めているだけだったので今回枠取れたの嬉しいです!
レベル分け的にはAdvanced Beginner の一人が同じくらいのレベルの人に向けて書いてるつもりです、やや内容に初等整数論が含まれます。
あわよくばより詳しい人にトンテンカンカンと直していただきたい感じですが。

モチベーション

電卓などで(1/3) の計算をした後、答えに3 を掛けたことがある人、そしてそのとき1 に戻らなかった経験ある人もきっと多いことだと思います。
今回の話はそれにちょっと関係している、かもしれません。

今回は体を取り上げます、体と言うのはいわゆる四則演算に閉じている代数体系のことですね。
やや詳しく言うと環のうち、零元(0)以外が掛け算に対する逆元を持っているようなもの、それを体と呼びます。
割り算が安全に行える環、ということですね。
特に興味があるのが

\mathbb{Q}

で表される有理数体や、それを連続化した

\mathbb{R}

なる実数体、そして代数閉体である複素数体

\mathbb{C}

というのがアリます。
グレブナー基底で最近有名なこの本では色々な証明の為の$\mathbb{C}$、(連続的な)絵やグラフを描いたりする為の$\mathbb{R}$、そして計算機の中では$\mathbb{Q}$という書き方をしていたのが印象的です。

有理数の計算はコストが高い?

例えば計算機の中で素朴に割り算をすると大きく情報を損なう可能性が高いのは皆さんご存知だと思います、上であげた電卓の例がそれに当たります。
言うまでも無いことですが、計算機上で無限の精度で割り算を行うことは素朴には無理で、例えば浮動小数点数に丸めて計算すると普通は望む結果では無いものになってしまいます。
素朴な解決策は、Int を二つ取ってきた新しいデータ型で有理数を表現するという方法です。
これはまぁうまくいきます、しかし計算の途中で分母分子のInt が桁落ちしてないことは一般には保障されません。
たとえ入力と出力の型がInt であっても、中間表現が桁落ちしていれば結果は望んでいたものではないものが得られることになる可能性があるわけで、これは困るわけです。
もちろん任意制度のInteger を使うという手もありますが、、、遅くコストがかさむことは容易に想像できます。

私自身が定量的に確認したわけではないのでアレですが、割り算、さらに言えば(掛け算に対する)逆元の計算を真面目に取り扱おうと思うとスペースを喰うと言うのはなんとなく感じてもらえると思います。

というわけで

ここで取り上げるのは以前の記事で取り上げたアルゴリズムの紹介とそれを上手に使った有理数復元方法です。
詳細はP. S. Wang, "A p-adic algorithm for univariate partial fractions", (1981)あるいはP. S. Wang, M. J. T. Guy, and J. H. Davenport, "p-adic reconstruction of rational numbers" (1982)などを参照してください。

それに付随して有限体をメインに取り上げます。

かいつまんで言いますと、全ての演算をInt の上(実はさらにもっと「安全」な領域)で行って、元の有理数を復元する方法のご紹介です。
勝手にコントロールされた桁落ちと呼んでいますが、ある意味桁落ちしない環境とも言えます。

環と体

いきなり脱線話ですが、体と聞いて皆さん何を想像されますか?
これは書き言葉だから出来ることですが、ここではこの漢字は「たい」と読まれます。
(残念ながらタイトルのように環とくっつけてしまうと誤読していただける機会は減るかと思いますが)
困ったことに英語ではfields となり物理の界隈ではなかなか体として相手に伝わらずもどかしい気持ちになるのです、ちなみに物理屋がfields というとまぁ間違いなく「場」だと思われるので有限体を指したくてfinite fields というと大体皆なにかエキゾチックなものだと誤解してくれます。

集合と演算(代数構造)

集合に構造を入れたもののうち、割合便利なものに大体名前が付いています、例えば群、環、体など。
平たく言えば群とは集合に(性質の良い)掛け算の構造が入っているもので、環は掛け算と足し算、体は四則演算が入ったものと言う感じに。
Haskeller の皆さんはさらにモノイドもご存知かもしれません、あとは圏とかも。
上で挙げた有理数体、実数体そして複素数体は(濃度は異なりますが)無限集合に体の構造が入っているものになります。
誤解を恐れずに言えば、電卓の桁落ち問題は有理数体の濃度が無限だということに起因します、従って有限な体があればもしかしたら、、、というのが希望であり、今回の話のモチベーションにつながるわけです。

剰余環

有限な環というのは割合簡単に見つかるものです、たとえばある正の整数に着目した余りの計算が環の構造を取ることは知られています。
例えば6 で割った余りでもって無限個の整数を分類するとなんと要素が6 個しかない有限の環が出来ます、これが剰余環(あるいは商環)と言われるもので

\mathbb{Z}_6

とここでは表します。
集合としては

\{0,1,2,3,4,5\}

だと思ったら良く、足し算引き算そして掛け算はすべて普通に実行した後$\mod 6$ をとるという風に約束するものとします、例えば

1+2 == 3 \mod 6 \\
4+5 == 3 \mod 6 \\
1-5 == 2 \mod 6 \\
3*4 == 0 \mod 6

など。
さてここでこの剰余環を体にすることが出来るかどうか?というのが浮かびうる疑問です。
それが出来れば第一ステップクリアですがその前に道具を一つ用意します。

最大公約数(gcd)

そう最大公約数、語呂が良くて覚えるのは覚えたのですが最小公倍数といつもごっちゃになってしまあれです、私だけかもしれませんが。
最大公約数を求めるアルゴリズムは古くから知られていてwikiには大変直観的な素敵なアニメーションも置いてありました。
実装も簡単で自前でも書けます、よく演習問題などになるのでは無いでしょうか?
細かいことを言うと私は大小判定、等号判定、そして引き算だけでの再帰的な定義の方が見通しが良くて好きです、これが上記のアニメーションに相当しているやつです。
以前にここで解いたように最大公約数は以下のように

myGCD.lhs
> myGCD :: Integral a => a -> a -> a
> myGCD a b
>   | b < 0 = myGCD a (-b)
> myGCD a b
>  | a == b = a
>  | b >  a = myGCD b a
>  | b <  a = myGCD (a-b) b 

大きい方から小さい方を引き続けて行ったらよいことになります、終了判定は両者が等しくなったときです、あるいはもう一度適用させて小さいほうが零でも大丈夫。

Bézout の補題

さてここで重要な補題を一つ、任意の整数a, b とそれらの最大公約数の間に以下の等号を満たすx, y が存在します。
すなわち

\forall a,b \in \mathbb{Z}, \exists x,y \in \mathbb{Z}, s.t. \\
a*x+b*y == \gcd(a,b)

以前ここで取り上げたこのアルゴリズムがこの二整数x, y を構成的に作る方法として知られています。

拡張されたユークリッドのアルゴリズム

Haskeller の方には(型には?)実物のコードの方が可読性が高いことだと思いますので、まず実物をば。

exGCD.lhs
> exGCD' :: Integral n => n -> n -> ([n], [n], [n], [n])
> exGCD' a b = (qs, rs, ss, ts)
>   where
>     qs = zipWith quot rs (tail rs)
>     rs = takeUntil (==0) r'
>     r' = steps a b
>     ss = steps 1 0
>     ts = steps 0 1
>     steps a b = rr
>       where rr@(_:rs) = a:b: zipWith (-) rr (zipWith (*) qs rs)
>         
> takeUntil :: (a -> Bool) -> [a] -> [a]
> takeUntil p = foldr func []
>   where
>     func x xs 
>       | p x       = []
>       | otherwise = x : xs
>
> -- a*x + b*y = gcd a b
> exGcd a b = (g, x, y)
>   where
>     (_,r,s,t) = exGCD' a b
>     g = last r
>     x = last . init $ s
>     y = last . init $ t

今回はこの「補助関数」と思われてたexGCD'とその「副産物」である整数リストたちが今回の話の道具になります。

証明、の前に

どういうアルゴリズムかというと、初期条件として

(r_0,s_0,t_0) := (a,1,0) \\
(r_1,s_1,t_1) := (b,0,1)

さらに

q_i := quot(r_{i-2}, r_{i-1})

とします、つまり余りです、例えば

quot(12,5) = 2

あとは再帰的に

r_i := r_{i-2} - q_i * r_{i-1} \\
s_i := s_{i-2} - q_i * s_{i-1} \\
t_i := t_{i-2} - q_i * t_{i-1} \\

と定義していきます。
停止条件は

r_k = 0

となることとします。

証明

何を証明するかというと、そう、Bézout の補題です。

\forall a,b \in \mathbb{Z}, \exists x,y \in \mathbb{Z}, s.t. \\
a*x+b*y == \gcd(a,b)

すなわち任意の整数a,b に対して$ a*x+b*y == \gcd(a,b) $ なる二整数x,y と$gcd(a,b)$ を、拡張されたユークリッドのアルゴリズムが構成的に作ることが出来ますよ!というのが主張です。

ここでは$a>b$ を仮定します、これは特に意味のある仮定ではありませんが一般性は失わない、単に簡便性のためです。
まず始めにこのアルゴリズムの停止性を述べねばなりません。
定義より

r_i := r_{i-2} - q_i * r_{i-1}

この数列$r_i$は単調に減少することが分かります、ここで$a>b$が生きるわけです。
停止条件は$r=0$なので、いつか必ずヒットするので停止性は大丈夫ということになります、すなわち有限回のステップで止まります。
このステップを$k$とします、すなわち$r_k = 0$。

次に各ステップでの$r$の挙動を見ましょう。
割り算の過程で

\gcd(r_{i-1}, r_i) = \gcd(r_{i-1}, r_{i-2} - q_i*r_{i-1}) \\
= \gcd(r_{i-1}, r_{i-2})

となることから、

\gcd(a,b) = \gcd(r_0, r_1) = ... = \gcd(r_{k-1},r_k=0)

が分かります、すなわち

\gcd(a,b) = r_{k-1}

次に各ステップでBézout 恒等式が成り立つこと、は読者の方への演習としましょう:

a*s_i + b * t_i = r_i

これらを合わせて

\gcd(a,b) = r_{k-1} = a*s_{k-1} + b*t_{k-1}

が成り立ちます、すなわち最大公約数と$x=s_{k-1}, y=t_{k-1}$が構成的に得られたことが分かりました。

逆元

ここからは$\gcd$ を使う立場に代わります。

もし二整数a, b が互いに素ならば$\gcd(a,b) == 1$を満たします、さらにこの時Bézout の補題から(あるいは拡張されたユークリッドのアルゴリズムから)

a*x+b*y == \gcd(a,b) == 1

なる二整数x, yがあります。
このとき例えば剰余環$\mathbb{Z}_b$の上で考えてみると

1 == a*x+b*y \mod b \\
1 == a*x \mod b

となり、なんと剰余環$\mathbb{Z}_b$の上では$x \mod b$ が$a \mod b$ の(掛け算に対する)逆元となるわけです。

例えば

Prelude> gcd 6 1
1
Prelude> gcd 6 5
1

などから$\mathbb{Z}_6$では1 は 1の逆元で5 は5 の逆元であることが分かります。

1*1 == 1 \\
5*5 == 4*6+1  

しかし残りの2,3,4 は逆元を持たないことが分かるので$\mathbb{Z}_6$ は体にならないことが分かります。

素数体

二整数a,b が互いに素ならば$a \mod b$ の逆元が存在し、拡張されたユークリッドのアルゴリズムに則ればその候補x が得られ、あとは$x \mod b$を計算すれば良いと言うのがミソです。

素数p に対して$\mathbb{Z}_p$を考えると、この要素

\{0,1 \cdots (p-1)\}

の内0 を除く全ての元a に関して

\gcd(a,p) == 1, a=1,2\cdots (p-1)

が成り立つので、この素数環は体になる!ということが分かります。
例えば$\mathbb{Z}_5$や$\mathbb{Z}_7$が体になることは

Prelude> map (gcd 5) [1..4]
[1,1,1,1]
Prelude> map (gcd 7) [1..6]
[1,1,1,1,1,1]

から分かると思います、問題は$\gcd$だけでは逆元はさっぱり分からないということですが、我々には拡張されたユークリッドのアルゴリズムがあります!したがって

*Ffield> map (exGCD 5) [1..4]
[(1,0,1),(1,1,-2),(1,-1,2),(1,1,-1)]
*Ffield> (-2) `mod` 5
3
*Ffield> 2 `mod` 5
2
*Ffield> (-1) `mod` 5
4

などから

1^{-1} == 1 \mod 5 \\
2^{-1} == 3 \mod 5 \\
3^{-1} == 2 \mod 5 \\
4^{-1} == 4 \mod 5 

が分かります、真ん中の二組は

2*3 == 5 + 1

から理解していただけると思います。

$\mathbb{Z}_7$では

*Ffield> map (exGCD 7) [1..6]
[(1,0,1),(1,1,-3),(1,1,-2),(1,-1,2),(1,-2,3),(1,1,-1)]
*Ffield> map (\(_,_,c) -> c) it
[1,-3,-2,2,3,-1]
*Ffield> map (`mod` 7) it
[1,4,5,2,3,6]

から

1^{-1} == 1 \mod 7 \\
2^{-1} == 4 \mod 7 \\
3^{-1} == 5 \mod 7 \\
4^{-1} == 2 \mod 7 \\
5^{-1} == 3 \mod 7 \\
6^{-1} == 6 \mod 7 

が分かります。

拡張された$\mod$

普通の$\mod p$ は整数を取って整数を返す関数です、すなわち

- \mod p :: \mathbb{Z} \to \mathbb{Z}_p

という"型"を取ると考えられます。
ここでこの定義域を有理数まで拡張しておきます、すなわち

q=\frac{n}{d} \mapsto (q \mod p) := (a \times ?(b^{-1} \mod p)) \mod p

で有理数q の$\mod p$ を定義します。

ここで

? :: \mathbb{Z}_p \to \mathbb{Z}

は自然な包含あるいは忘却関数で、$\mathbb{Z}_p$の構造を忘れてただの整数だと思いましょう、という関数です。

今回のミッションはこの$\mathbb{Z}_p$ の元である、すなわち

0 \leq a \leq (p-1)

を満たす、「有理数q の$\mod p$ での像」から元のq を推測しましょう、という話です。

有理数の復元(素朴な方法)

さて拡張ユークリッドの副産物として示すことが出来るものとして、各ステップ$i$ においてBezout の恒等式が 成り立つというのがあります。

a* s_i + b*t_i = r_i

素朴な推測

ここで我々の問題設定として$b=p$ を素数としてとり、かつ$a$ はターゲットの有理数の$Z_p$ の上での像だとします、すなわち

a = (q \mod p)

この時各ステップでのBezout 恒等式の両辺を$Z_p$ で考えると

a * s_i = r_i \mod p

従って

a = (\frac{n}{d} \mod p )

は$Z_p$ の上で

(\frac{n}{d} \mod p) = r_i * {s_i}^{-1} \mod p

を満たすことから $\frac{n}{d}$は各ステップで

\frac{s_i}{r_i}

と"推測される"というわけです。
ここでWang が示したのは

s_i^2,r_i^2 \lessapprox p

なるステップ$i$ において元の有理数は

q = \frac{s_i}{r_i}

である、ということです。

問題点

このアルゴリズムの問題点は推測の精度が

s_i^2,r_i^2 \lessapprox p

という条件のせいで平方根で抑えられているという点です。

例えば素数として10007を取ったとすれば推測できる有理数は分母分子ともに$10^2$程度までだということです。

中国の剰余定理と有理数の復元

主張と証明はwikiの該当ページに任せるとします。
ミソは2つの互いに素な$p_1,p_2$(素数でなくても、互いに素ならok)があって、かつ$Z_{p_1}$と$Z_{p_2}$の上での値

a_1 = (q \mod p_1) \\
a_2 = (q \mod p_1)

が与えられていれば

a = (q \mod p_1*p_2)

なる$a$ をただ1つ決めてあげられるという点につきます。

問題設定

まず比較的大きな素数を幾つか用意しておきます、例えば

Prelude Data.Numbers.Primes> take 10 $ dropWhile (<10^4) primes
[10007,10009,10037,10039,10061,10067,10069,10079,10091,10093]

それぞれの素数における商体$Z_p$ に対して

(a_1, p_1)\\
(a_2, p_2)\\
\vdots

なるデータがあるとします、なんでこんなもんを考えるの?それは最後に少しだけ書きますがまぁそういう問題設定だと思ってください。
単なる有理数の計算なら整数を2つ分とるデータ型に、色々ふさわしい演算を導入してあげたらいいだけなのでこんな質面倒くさいことを考える必要はありません。

ここで中国の剰余定理の定理により最初の2つの像と素数から新たに

(a,p_1*p_2)

なるデータを作ります、$a$ のレシピはwikiの該当ページ参照のこと。
この新しいデータを拡張ユークリッドにかけますと今度は分母分子ともに

s_i^2,r_i^2 \lessapprox (p_1*p_2)

という範囲で決めることができます。
これで足らなければ更にもう一つ素数を取ってきて

s_i^2,r_i^2 \lessapprox (p_1*p_2*p_3)

という風にどんどん制限をゆるくしていくことができます。

例えば私の環境(すなわち素数として10^4 程度のものから取った場合)ならば上記の素数なら4つか5つか取ればInt において充分なことが分かります。

Prelude> maxBound :: Int
9223372036854775807
Prelude> 2^63 -1
9223372036854775807
Prelude> logBase 10 it
18.964889726830812

従って例えば3つ取って復元した値と4つの値、そして5つの値が一致していれば実用上ユニークに元の有理数が復元できてるとして良いということです。

これは嘘です、最後に追記で訂正されたバージョンを載せておきます。
オーダーの読みとしては

s_i^2,r_i^2 \lessapprox (p_1*\cdots)

の右辺がInt の最大値の二乗を上から抑えたら良いので、すなわち10 個が目安になります。
なぜなら

p_1, \cdots \in o(10^4)

と取った素数と

Prelude> logBase 10 (2^63-1)
18.964889726830812

を比べたら良いので、大体10個の積でInt をカバーできるはずです。

有理数の復元

以下が今回のコードの全景になります、無事quickCheck も通りました。
中国の剰余定理の計算のところだけ$Integer$ を使ってますが、それ以外の演算は基本的にマシンサイズの$Int$ で出来るのがミソです。

> module Ffield where

> import Data.Ratio 
> import Data.Maybe
> import Data.Numbers.Primes
> import Test.QuickCheck

> coprime :: Integral a => a -> a -> Bool
> coprime a b = gcd a b == 1

Consider a finite ring
  Z_n := [0..(n-1)]
of some Int number.
If any non-zero element has its multiplication inverse, then the ring is a field:

> -- Our target should be in Int.
> isField :: Int -> Bool
> isField = isPrime

> exGCD' :: (Integral n) => n -> n -> ([n], [n], [n], [n])
> exGCD' a b = (qs, rs, ss, ts)
>   where
>     qs = zipWith quot rs (tail rs)
>     rs = takeUntil (==0) r'
>     r' = steps a b
>     ss = steps 1 0
>     ts = steps 0 1
>     steps a b = rr
>       where 
>         rr@(_:rs) = a:b: zipWith (-) rr (zipWith (*) qs rs)
>
> takeUntil :: (a -> Bool) -> [a] -> [a]
> takeUntil p = foldr func []
>   where
>     func x xs 
>       | p x       = []
>       | otherwise = x : xs
>
> -- a*x + b*y = gcd a b
> exGCD :: Integral t => t -> t -> (t, t, t)
> exGCD a b = (g, x, y)
>   where
>     (_,r,s,t) = exGCD' a b
>     g = last r
>     x = last . init $ s
>     y = last . init $ t
>
> -- a^{-1} (in Z_p) == a `inversep` p
> inversep :: Integral a => a -> a -> Maybe a -- We also use in CRT.
> a `inversep` p = let (g,x,_) = exGCD a p in
>   if (g == 1) 
>     then Just (x `mod` p) -- g==1 <=> coprime a p
>     else Nothing
>
> inversesp :: Int -> [Maybe Int]
> inversesp p = map (`inversep` p) [1..(p-1)]
>
> -- A map from Q to Z_p.
> -- p should be prime.
> modp :: Ratio Int -> Int -> Maybe Int
> q `modp` p 
>   | coprime b p = Just $ (a * (bi `mod` p)) `mod` p
>   | otherwise   = Nothing
>   where
>     (a,b)   = (numerator q, denominator q)
>     Just bi = b `inversep` p
>
> -- This is guess function without Chinese Reminder Theorem.
> guess :: Integral t => 
>          (Maybe t, t)       -- (q `modp` p, p)
>       -> Maybe (Ratio t, t)
> guess (Nothing, _) = Nothing
> guess (Just a, p) = let (_,rs,ss,_) = exGCD' a p in
>   Just (select rs ss p, p)
>     where
>       select :: Integral t => [t] -> [t] -> t -> Ratio t
>       select [] _ _ = 0%1
>       select (r:rs) (s:ss) p
>         | s /= 0 && r*r <= p && s*s <= p = r%s
>         | otherwise = select rs ss p
>
> -- Hard code of big primes
> bigPrimes :: [Int]
> bigPrimes = take 10 $ dropWhile (<10^4) primes
> -- bigPrimes = take 10 $ dropWhile (<10^6) primes
> -- 10 is just for "practical" reason.

  *Ffield> let knownData q = zip (map (modp q) bigPrimes) bigPrimes
  *Ffield> let ds = knownData (12%13)
  *Ffield> map guess ds
  [Just (12 % 13,10007)
  ,Just (12 % 13,10009)
  ,Just (12 % 13,10037)
  ,Just (12 % 13,10039) ..

  *Ffield> let ds = knownData (112%113)
  *Ffield> map guess ds
  [Just ((-39) % 50,10007)
  ,Just ((-41) % 48,10009)
  ,Just ((-69) % 20,10037)
  ,Just ((-71) % 18,10039) ..

--
Chinese Remainder Theorem, and its usage

> imagesAndPrimes :: Ratio Int -> [(Maybe Int, Int)]
> imagesAndPrimes q = zip (map (modp q) bigPrimes) bigPrimes

Our data is a list of the type
  [(Maybe Int, Int)]
In order to use CRT, we should cast its type.

> toInteger2 :: [(Maybe Int, Int)] -> [(Maybe Integer, Integer)]
> toInteger2 = map helper
>   where 
>     helper (x,y) = (fmap toInteger x, toInteger y)
>
> crtRec':: Integral a => (Maybe a, a) -> (Maybe a, a) -> (Maybe a, a)
> crtRec' (Nothing,p) (_,q)       = (Nothing, p*q)
> crtRec' (_,p)       (Nothing,q) = (Nothing, p*q)
> crtRec' (Just a1,p1) (Just a2,p2) = (Just a,p)
>   where
>     a = (a1*p2*m2 + a2*p1*m1) `mod` p
>     Just m1 = p1 `inversep` p2 
>     Just m2 = p2 `inversep` p1
>     p = p1*p2
>
> matches3 :: Eq a => [Maybe (a,b)] -> Maybe (a,b)
> matches3  (b1@(Just (q1,p1)):bb@((Just (q2,_)):(Just (q3,_)):_))
>   | q1==q2 && q2==q3 = b1
>   | otherwise        = matches3 bb
> matches3 _ = Nothing

  *Ffield> let ds = imagesAndPrimes (1123%1135)
  *Ffield> map guess ds
  [Just (25 % 52,10007)
  ,Just ((-81) % 34,10009)
  ,Just ((-88) % 63,10037) ..

  *Ffield> matches3 it
  Nothing

  *Ffield> scanl1 crtRec' ds

  *Ffield> scanl1 crtRec' . toInteger2 $ ds
  [(Just 3272,10007)
  ,(Just 14913702,100160063)
  ,(Just 298491901442,1005306552331) ..

  *Ffield> map guess it
  [Just (25 % 52,10007)
  ,Just (1123 % 1135,100160063)
  ,Just (1123 % 1135,1005306552331)
  ,Just (1123 % 1135,10092272478850909) ..

  *Ffield> matches3 it
  Just (1123 % 1135,100160063)

The final reconstruction function takes a list of Z_p values and returns the three times matched guess.

> reconstruct :: [(Maybe Int, Int)] -> Maybe (Ratio Integer, Integer)
> reconstruct = matches3 . map guess . scanl1 crtRec' . toInteger2

--

> aux :: [(Maybe Int, Int)] -> Ratio Int
> aux = fromRational . fst . fromJust . reconstruct
>
> prop_rec :: Ratio Int -> Bool
> prop_rec q = q == aux ds
>   where
>     ds = imagesAndPrimes q

  *Ffield> quickCheck prop_rec 
  +++ OK, passed 100 tests.

終わりに

最初に述べたように、Int で与えられた数字でも演算の途中で桁落ちする可能性があります。

Prelude> let a = maxBound :: Int
Prelude> a
9223372036854775807
Prelude> a*2 :: Int
-2

このように、桁落ちした場合本質的な情報は完全に消えてしまうわけですが、素数体と中国の剰余定理を組み合わせ、かつ二項演算の際に毎回$\mod$ を施すことによって、コントロールされた桁落ち環境が実現されるわけです。
あるいは元から要素が有限個なので桁落ち「しない」環境といっても良いと思います。

ここで取り上げた方法は演算は全て安価な有限精度のInt の上で行い、最後復元作業の時だけ任意精度のInteger で行うという方法になっています。
したがって、全演算を任意制度の有理数Ratio Integer で行うより安価でありそうな雰囲気を感じていただけると思います。

一体どうして?

ここでどういう問題を考えたかというと、実は有理数係数の多項式

f :: Q \to Q

などを(我々のプロジェクトでは)考えています。
なお、このとき$f$ がブラックボックスで与えられていて係数などを知り得ない状況に置かれてるとします。
(例えば離散点における$f$ の値だけ知ってる(グラフが与えられてる)とか)
多項式ならニュートン補間とか、有理関数でもThiele の方法などが知られていて、十分な点を与えさえすれば係数を決めることができます。
この補完法をInt できまる有理数でやると途中で桁落ちとかが怖いので、いろんな有限体Z_p$ の上で行い、結果各係数の

(a_1, p_1)\\
(a_2, p_2)\\
\vdots

なるデータが得られている、という状況を想定してます。
$Z_p$ でももちろん $ \mod $ のせいで桁落ちしているといえばしているわけですが、きちんとコントロールされているので、中国の剰余定理と組み合わせることできちんと有理数係数を復元することが出来る、という寸法です。

追記

何の気なしにquickCheck を何度かタイプして悦に入ってましたら、時折こけることがあることが分かりました。
例えば

> quickCheck prop_rec 
*** Failed! Exception: 'Maybe.fromJust: Nothing' (after 9 tests and 1 shrink): 
56307372244017 % 5226552907916

これは

> 2^2*10009*130546331
5226552907916

ということからわかるように$Z_ {10009}$ の上では$\frac{1}{0}$ の型をした不定形になってしまいます、いわばフェイクの無限大です。
したがって下処理としてNothing を取り除いておけば大丈夫。。。と思いましたが、これは嘘でした。
もちろんNothing を取り除く処理も大事ですが一番大事なのはInt を覆うことで、このためには10回の合致が必要です。
(この10 という数字は私の素数のチョイスとInt の最大値によって決まっています。)

> reconstruct :: [(Maybe Int, Int)] -> Maybe (Ratio Integer)
> reconstruct = matches 10 . makeList
>   where
>     matches n (a:as)
>       | all (a==) $ take (n-1) as = a
>       | otherwise                     = matches n as
>     makeList = map (fmap fst . guess) . scanl1 crtRec' . toInteger2 . filter (isJust . fst)

--

> prop_rec :: Ratio Int -> Bool
> prop_rec q = Just q == answer
>   where
>    answer :: Maybe (Ratio Int)
>    answer = fmap fromRational . reconstruct $ ds
>    ds = imagesAndPrimes q

  *Ffield> quickCheckWith stdArgs { maxSuccess = 100000 } prop_rec 
  +++ OK, passed 100000 tests.

念のため多めに取りましたが、今度こそ大丈夫なはずです。

bra_cat_ket
A Haskell, Maxima lover, PhD in physics. About to leave qiita, perhaps.
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした