複素数

Pythonは複素数を直接書くことができます。

Python 2.7.12 (default, Dec  4 2017, 14:50:18)
[GCC 5.4.0 20160609] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> a=0+1j
>>> a
1j

実数部が０のとき省略できます。

>>> a=1j

>>> a

1j

虚数部が０のときこれも省略できます。

>>> b=1

>>> b

1

これはただの実数ですが複素数と組み合わせるとちゃんと複素数の計算をしてくれます。

>>> a-b

(-1+1j)

複素数は距離の計算をとても楽にしてくれます。

a,b　２点間の距離は差をとって絶対値を計算すればよい。

>>> d=abs(a-b)

>>> d

1.4142135623730951

これは√２です。

これをｂに代入します。

>>> b=d
>>> b
1.4142135623730951
>>> d=abs(a-b)
>>> d
1.7320508075688774

これは√３です。

このように計算結果を次々に代入していけば全ての整数の平方根がもとめられます。

(複素数の計算の内部で　sqrt を使っているのでは？　というツッコミはあると思いますが無視してください）

Haskellではこのような記法はできません。

Data.Complex　モジュールの

(:+)

というコンストラクターで複素数のインスタンスをつくります。

Haskellではコンストラクターは関数であり、（：＋）は演算子として定義されているので次のようにします。

GHCi, version 7.10.3: http://www.haskell.org/ghc/  :? for help

Prelude> import Data.Complex

Prelude Data.Complex> let a=0:+1

Prelude Data.Complex> a

0 :+ 1

Prelude Data.Complex> let b=1:+0

Prelude Data.Complex> b

1 :+ 0

ちょっと不格好ですが・・・

Haskell　は、複素数の絶対値を複素数で返します。

Prelude Data.Complex> abs(a-b)

1.4142135623730951 :+ 0.0

これは abs が

abs :: Num a => a -> a

と定義されているからです。

かわりにmagnitudeをつかいます。

Prelude Data.Complex> :t magnitude

magnitude :: RealFloat a => Complex a -> a

Prelude Data.Complex> magnitude(a-b)

1.4142135623730951

ちなみに　Data.Complex　では　abs を

abs z               =  magnitude z :+ 0

と再定義しているようです。(^_^;)

N−１　の平方根がわかれば　N　の平方根が計算できることから再帰的に定義できます。

--ｒ．ｈｓ

import Data.Complex

dist a b = magnitude (a-b)

a = 0 :+ 1

r 0 = 0

r n = dist a b

    where b = r (n-1) :+ 0

Prelude> :l ｒ.hs
[1 of 1] Compiling Main             ( r.hs, interpreted )
Ok, modules loaded: Main.
*Main> map r [0..9]
[0.0,1.0,1.4142135623730951,1.7320508075688774,2.0,2.23606797749979,2.4494897427831783,2.6457513110645907,2.8284271247461903,3.0000000000000004]

容易に想像できると思いますがNが大きくなると時間がかかるようになります。

*Main> :set +s
*Main> r 1000000
1000.0000000000299
(7.74 secs, 1,673,363,512 bytes)

このような再帰は fold　で置き換えることができます。

*Main> let r' n = foldl (\x _ -> dist a (x :+ 0)) 0 [1..n]

*Main> r' 1000000

1000.0000000000299

(5.28 secs, 1,474,441,216 bytes)

さらに　Data.List　モジュールの　foldl' 　を使うと

*Main> import Data.List

*Main Data.List> let r'' n = foldl' (\x _ -> dist a (x :+ 0)) 0 [1..n]

*Main Data.List> r'' 1000000

1000.0000000000299

(2.52 secs, 1,293,766,728 bytes)

magnitude　自体が重い処理なのでそこそこかかります。

fold　を使ったコードは最初に説明した手続き的処理をそのまま置き換えている点に注目です。

手続き的なscanlの説明

http://hhg2the-haskell.blogspot.com/2017/10/scanl.html

参照

余談ですが　（：＋）　の型は　

Prelude Data.Complex> :t (:+)
(:+) :: a -> a -> Complex a

です。引数の型は何でもいいようです。ですので

Prelude Data.Complex> "hello" :+ "world"
"hello" :+ "world"
Prelude Data.Complex> :t it
it :: Complex [Char]

もちろん演算とかはできません。

ランダム法

「Haskell98標準のrandomパッケージは、1.長くない周期、2.Haskell98由来の古いAPIの制約で内部状態を毎回コピーする、3.そのため遅い、4.CPU時間と現在時間で初期化(ポケモンかな？？？)など、過去に標準だったため多くの環境で使用出来ると考えられる、という点以外に余りお勧めできる所が有りません。」
（Haskellの乱数事情より）
https://qiita.com/philopon/items/8f647fc8dafe66b7381b

だそうですがここでは便利さをとり
 System.Random
モジュールを使います。

 System.Randomモジュールはcabalコマンドでパッケージを導入しなければなりません。
まず

sudo apt install cabal-install

でcabalをインストール。

$ cabal　update　-ｖ
$ cabal install random　-ｖ

でImportできるようになります。
結構時間がかかるので　-v オプションをつけて経過を見ながら。

なお表題の「ランダム法」はランダムを使っているという意味であり「モンテカルロ法」とは関係ありません。

基本のコードは

---

これで、未知数が一つ、解の範囲があらかじめ推測できるという条件の方程式を解くことができます。

上のコードに

root n x = evalStateT (try (pf n x)) (1,x) >>= return
pf n x g = g^n < x
main = root 2 2 >>= print

を付け加えて実行しますと

$ runghc -XFlexibleContexts  randomMethod.hs
(1.4142135623728276,1.4142135623731693)

---

3角形の３辺の長さがわかっていれば、その形と大きさが決まることから面積を計算できるはずです。
ヘロンの公式です。
ちなみにヘロンの公式はHaskellでは極めて簡潔にこう書けます。
heron a b c = sqrt (s*(s-a)*(s-b)*(s-c)) where s = (a+b+c) / 2

余談ですがどんな複雑な形の土地でも３角形に分割することで巻き尺１本で面積計算ができます。
（ただし平面地）

ヘロンの公式が思い出せなくてもランダム法で直感的に解くことができます。

三角形の３辺の長さが　a,b,c で

ｂを複素平面上の実数軸においたとき、

三角形の頂点は極座標で

（a,θ）

で表せます。

このθをランダム法で探ればよいのです。

---

$ runghc -XFlexibleContexts  heron.hs
6.0

円に内接する正６角形と、同じ円に外接する正４角形を考えると、
円周率は３と４の間に存在するとわかります。
３と４の間の実数のどれかはπの近似値になるはずです。
では、それをランダム法で求められるでしょうか。
残念ながらできません。円周率を求める方程式が存在しないからです。

複素数を使えば間接的にπを求めることはできます。
複素数　ｃ＝-１＋0i　を極座標で表しますと　（1、π）
位相を未知数θとして　θを　３と　4　の間で、さまざまに変化させ　ｃ　と近似する値をみつければよい。

でもこれはチートなのでは？

フィボナッチ数

---

f6 の実行例

*Main> last $ take 10 f6
[55,34,21,13,8,5,3,2,1,1,0]

これは

foldl (\x@(a:b:_) _ -> (a+b):x) [1,0] [0..８]

と同じです。

フィボナッチ数列の隣り合う 2 項の比は黄金比に収束^Wします。

$ runghc fibo2.hs | runghc gr.hs

~
~
~
1.618033988749895
1.618033988749895
1.618033988749895
~
~
~
Ctrl-Cを押して終了。

πの計算

ライプニッツの公式

---

---

実行結果
＄　runghc ｌ.hs
3.140592653839793

マチンの公式によるπ

---

参考にしたページ
http://handasse.blogspot.com/2011/12/1.html

---

Parallel版

1万桁では差がでないので
１０万桁ほど計算します。

コンパイルします。

ghc  --make -O2 pm.hs -rtsopts -threaded -eventlog

CPUはCore2Duo。デュアルコアなので　-N2 オプションで実行します。
$./pm +RTS -N2 -l
threadscope で動作を確認します。
$ threadscope pm.eventlog

僅かですが速くなっています。

並列計算を効率よく利用するには最初から並列化を意識してコードを書く必要があるようです。

---

spigot アルゴリズムをLoopで書いたコード
kk62526　さんが１０進Basicで書いたコードをそのまま移植しています。

グローバル変数　「i(注）」に現在何桁目を計算しているかが保持されてますので（ただし＋２された数字）
(注：　本当は無名ですが、こう捉えたほうが分かりやすいです。）
if (i-2) `mod` 100 == 0　then ・・・・・・
のようなコードを挿めば１００桁ごとに区切ったりできます。
実行結果
 $ runghc p.hs > a
$ less a
3 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679()
100 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196()
200 4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273()
300
~
~
~
（略）
ちなみに１００００桁を計算し終えるのに２秒少々でした。
ギボンズ博士が spigot アルゴリズムに対し
"not competitive with state-of-the-art arithmetic-geometric mean algorithms"
と言うとおり決して高速とはいえませんが、いきなり１桁目から表示を始めますので、
体感的には超高速です。
(注：１００００桁の計算では８万数千桁のIntegerが３つ、作られます。）

ギボンズ博士のspigot アルゴリズムは
https://www.cs.ox.ac.uk/jeremy.gibbons/publications/spigot.pdf
コードだけ紹介しますと

平方根、立方根

---

2の平方根

$ runghc nroot.hs 2 2 | tee a

2の立方根

$ runghc nroot.hs 2 ３ | tee ｂ

猛列な勢いで計算しますので適当なファイルに保存もしています。（Ctrl-Cで停止）

---

なお√２に限ればこう書けます。

One liner:

loop a b n = return (a*2+1) >>= \c -> if c<b then loop (a+1)(b-c)(n+1) else putStr (show n) >> loop (a*10)(b*100)0; main = loop 1 1 1

 

---

平方根は正方形の色紙を百分割して並べることで、
立方根は正四角形のサイコロを千分割して貼り付けていくことで深い桁まで計算をすすめることができます。

参照ページ
http://senjounosantafe.blogspot.com/2018/05/2.html

stern-brocot 木

stern-brocot 木
Stern-Brocot 木の「親」の値は数直線の区間です。
頂点の親「根」の値は
０〜∞
左端が０で、右は無限大です。
それを分数で０/１、１/０とそれぞれ表します。
さらにHaskellではタプルで
(0,1) ,(1,0)
と表します。
つまりHaskellでStern-Brocot 木の根はタプルのタプルで
（(0,1) ,(1,0))
と表せます。
親はいつでも自動で左の子、右の子を生成することができます。
Haskellで書くとこのような関数になります。

sb ((x,y),(z,w)) = (a,b)  where  a = ((x,y),c) ; b = (c,(z,w)); c = (z+x,y+w)

aが左の子ｂが右の子です。左の子の右端と右の子の左端は同じです。分割点になります。

main = do
    let a = (0,1)
    let b = (1,0)
    let r = (a,b)
    print r
    print $ sb r

((0,1),(1,0))
(((0,1),(1,1)),((1,1),(1,0)))

数直線の区間
０〜∞
を
０〜１
と
１〜∞
の２つの区間に分割したのと同じです。

もう一世代増やしてみましょう。

main = do
    let a = (0,1)
    let b = (1,0)
    let r = (a,b)
   
    let (a,b) = sb r
    print $ sb a
    print $ sb b

(((0,1),(1,2)),((1,2),(1,1)))
(((1,1),(2,1)),((2,1),(1,0)))

これは
数直線の区間
０〜１
を
０〜０．５
と
０．５〜１
に
１〜∞
を
１〜２
と
２〜∞
に分割したのと同じです。
これを繰り返せば数直線を無限に分割できます。
Stern-Brocot 木の木構造としての性質等は
Spaghetti Source - Stern-Brocot 木
http://www.prefield.com/algorithm/math/stern_brocot_tree.html

を参照してください。
さてStern-Brocot木を２分探索してゆけば全ての少数を分数に変換できるのですが、
そのためにはStern-Brocot木を前もって作成しておく必要はありません。
与えられた少数と分割点とを比較しながら片側の子のみを辿ってゆけばよい。

この辺はコードにしたほうがわかりやすい。

---

---

実行結果

245850922 % 78256779

---