Skip to content

Instantly share code, notes, and snippets.

@genkuroki

genkuroki/find pi by julia.ipynb

Last active February 21, 2023 08:22
Show Gist options
  • Save genkuroki/45e47f56799aee57c3dc11cda6df869a to your computer and use it in GitHub Desktop.
Save genkuroki/45e47f56799aee57c3dc11cda6df869a to your computer and use it in GitHub Desktop.
Download ZIP
Julia/find pi by julia.ipynb
Display the source blob
Display the rendered blob
Raw
find pi by julia.ipynb
{
"cells": [
{
"metadata": {},
"cell_type": "markdown",
"source": "# Juliaにおける円周率のモンテカルロ法による計算"
},
{
"metadata": {},
"cell_type": "markdown",
"source": "黒木玄\n\n2017年7月17日更新 (2017年7月作成)"
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## 要約\n\n円周率を求めるモンテカルロ法のコードを Julia では `sum()` 函数を使うことによって1行で書けるが、`for` ループを使った素朴なコードよりも遅いし、 `@parallel` マクロによる並列化の高速化もできなくなる。しかし、Juliaのマクロ機能を使えば、`for` ループを使った高速化と `@parallel` マクロを使った並列化と「1行プログラミング」をすべて同時に実現することができる。"
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## 実行環境\n\nこのノートブックの実行環境は以下の通り。\n\n後で `@parallel` マクロを使うので、julia を -p auto オプション付きで起動するようにしたカーネルを Jupyter に登録して利用している。"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "versioninfo()",
"execution_count": 1,
"outputs": [
{
"output_type": "stream",
"text": "Julia Version 0.6.0\nCommit 903644385b* (2017-06-19 13:05 UTC)\nPlatform Info:\n OS: Windows (x86_64-w64-mingw32)\n CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz\n WORD_SIZE: 64\n BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)\n LAPACK: libopenblas64_\n LIBM: libopenlibm\n LLVM: libLLVM-3.9.1 (ORCJIT, haswell)\n",
"name": "stdout"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## `sum()` 函数を使った1行プログラミング\n\n次のように `sum()` 函数を使って書けばモンテカルロ法が1行で書ける。\n\n`sum(f, 1:n)` は `f(i)` の `i=1` から `n` までの和を意味する。\n\n`ifelse(condition, a, b)` の値は `condition` が真ならば `a` に偽ならば `b` になる。\n\n`i::Int -> f(i)` は整数 `i` を `f(i)` に対応させる函数を意味する。\n\n次の数式をそのまま Julia のコードに直しただけだと考えると分かり易い：\n$$\\displaystyle\nf(n)=\n\\frac4n \\sum_{i=1}^n\n\\left\\{\n\\begin{array}{c}\n\\text{if}\\;\\; X_i^2+Y_i^2\\leqq 1 \\\\\n \\text{then}\\; 1, \\\\\n \\text{else}\\; 0 \\\\\n\\end{array}\n\\right\\}\n$$\nここで $X_i,Y_i$ は区間 $[0,1]$ 上の一様分布に従う独立な確率変数達である。\n\nしかし、後で示すように、 `for` ループを使うよりも遅いし、並列化による高速化も容易ではない。"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "findpi_sum(n) = (4.0/n) * sum(i::Int -> ifelse(rand()^2 + rand()^2 <= 1.0,1,0),1:n)\nfindpi_sum(10^5)",
"execution_count": 1,
"outputs": [
{
"data": {
"text/plain": "3.14488"
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@time findpi_sum(10^9)",
"execution_count": 10,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": " 8.804012 seconds (6 allocations: 192 bytes)\n"
},
{
"data": {
"text/plain": "3.141682668"
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## `for` ループによる高速化\n\n`for` ループを使って書くとこうなる。\n\n上の `sum()` 函数を使った場合よりも倍以上速い。\n\nしかし、上の `sum()` 函数を使った場合と違って数式通りのコードになっていない。"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function findpi(n)\n count = 0\n for i in 1:n\n count += ifelse(rand()^2+rand()^2<=1.0,1,0)\n end\n return 4.0 * count / n\nend\nfindpi(10^5)",
"execution_count": 3,
"outputs": [
{
"data": {
"text/plain": "3.14496"
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@time findpi(10^9)",
"execution_count": 11,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": " 3.967847 seconds (6 allocations: 192 bytes)\n"
},
{
"data": {
"text/plain": "3.141732276"
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## `@parallel` マクロによる並列化\n\n`@parallel` マクロで並列化すると次のようになる。\n\nこの実行環境では並列化すると4倍近く速くなる。\n\nしかし、並列化のための記述が直接的に見えるせいで、コードが読み難くなった。"
},
{
"metadata": {
"code_folding": [],
"trusted": true
},
"cell_type": "code",
"source": "function findpi_parallel(n)\n count = @parallel (+) for i in 1:n\n ifelse(rand()^2+rand()^2<=1.0,1,0)\n end\n return 4.0 * count / n\nend\nfindpi_parallel(10^5)",
"execution_count": 5,
"outputs": [
{
"data": {
"text/plain": "3.149"
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@time findpi_parallel(10^9)",
"execution_count": 12,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": " 1.072139 seconds (1.20 k allocations: 99.750 KiB)\n"
},
{
"data": {
"text/plain": "3.14148082"
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## マクロによる「高速化・並列化」と「1行プログラミング」の同時実現\n\nそこで `@sum` マクロを次のように定義しよう。\n\n`@sum i 1:n f(i)` は `f(i)` を `i=1` から `n` まで足し上げた結果を意味する。\n\n次の数式をそのまま Julia のコードに直しただけだと考えると分かり易い：\n$$\\displaystyle\nf(n)=\n\\frac4n \\sum_{i=1}^n\n\\left\\{\n\\begin{array}{c}\n\\text{if}\\;\\; X_i^2+Y_i^2\\leqq 1 \\\\\n \\text{then}\\; 1, \\\\\n \\text{else}\\; 0 \\\\\n\\end{array}\n\\right\\}\n$$\nここで $X_i,Y_i$ は区間 $[0,1]$ 上の一様分布に従う独立な確率変数達である。\n\n`@sum` マクロは上の `for` ループによる高速化と `@parallel` マクロによる並列化を `sum()` 函数を使った場合と同じように1行で実現するための手段を与える！"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "macro sum(i, iter, expr)\n quote\n begin\n s = @parallel (+) for $(esc(i)) in $(esc(iter))\n $(esc(expr))\n end\n end\n end\nend",
"execution_count": 5,
"outputs": [
{
"data": {
"text/plain": "@sum (macro with 1 method)"
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "findpi_macrosum(n::Int) = (4.0/n) * (@sum i 1:n ifelse(rand()^2 + rand()^2 <= 1.0, 1, 0))\nfindpi_macrosum(10^5)",
"execution_count": 6,
"outputs": [
{
"data": {
"text/plain": "3.1436400000000004"
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@time findpi_macrosum(10^9)",
"execution_count": 33,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": " 1.067360 seconds (1.21 k allocations: 99.484 KiB)\n"
},
{
"data": {
"text/plain": "3.141621944"
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@time findpi_macrosum(4*10^9)",
"execution_count": 35,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": " 4.280931 seconds (1.28 k allocations: 101.734 KiB)\n"
},
{
"data": {
"text/plain": "3.1415888770000002"
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"_draft": {
"nbviewer_url": "https://gist.github.com/45e47f56799aee57c3dc11cda6df869a"
},
"gist": {
"id": "45e47f56799aee57c3dc11cda6df869a",
"data": {
"description": "Julia/find pi by julia.ipynb",
"public": true
}
},
"kernelspec": {
"name": "julia-0.6-pauto",
"display_name": "Julia 0.6.0 procs auto",
"language": "julia"
},
"language_info": {
"file_extension": ".jl",
"name": "julia",
"mimetype": "application/julia",
"version": "0.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Display the source blob
Display the rendered blob
Raw
find pi by numpy.ipynb
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment