タブを空白文字に展開するプログラムを作ります。プログラムを簡単にするためファイルはアスキーコードだけで書かれていて、タブの値は 8 個に限定します。コマンドラインでファイル名が省略された場合は標準入出力を使うものとします。なお、Unix 系 OS には同等の機能を持つコマンド expand があります。
リスト : タブを空白に展開する (detab.rs)
use std::fs::File;
use std::io::prelude::*;
use std::io::{self, BufReader};
fn detab<R: Read>(reader: BufReader<R>) {
for xs in reader.lines() {
let mut n = 0;
for c in xs.unwrap().chars() {
if c == '\t' {
loop {
print!("{}", ' ');
n += 1;
if n % 8 == 0 { break; }
}
} else {
print!("{}", c);
n += 1;
}
}
println!("");
}
}
fn main() {
let args: Vec<_> = std::env::args().collect();
if args.len() == 1 {
detab(BufReader::new(io::stdin()));
} else {
match File::open(&args[1]) {
Ok(file) => detab(BufReader::new(file)),
Err(_) => println!("{} not found", &args[1])
}
}
}
空白文字をタブに置き換えるプログラムを作ります。プログラムを簡単にするためファイルはアスキーコードだけで書かれていて、タブの値は 8 個に限定します。コマンドラインでファイル名が省略された場合は標準入出力を使うものとします。なお、Unix 系 OS には同等の機能を持つコマンド unexpand があります。
リスト : 空白をタブに置換する (entab.rs)
use std::fs::File;
use std::io::prelude::*;
use std::io::{self, BufReader};
fn entab<R: Read>(reader: BufReader<R>) {
for xs in reader.lines() {
let s = xs.unwrap();
let mut s0 = s.as_str();
loop {
if s0.len() < 8 {
println!("{}", s0);
break;
}
let (s1, s2) = s0.split_at(8);
let s3 = s1.trim_end();
if s3.len() < 8 {
print!("{}\t", s3);
} else {
print!("{}", s3);
}
s0 = s2;
}
}
}
fn main() {
let args: Vec<_> = std::env::args().collect();
if args.len() == 1 {
entab(BufReader::new(io::stdin()));
} else {
match File::open(&args[1]) {
Ok(file) => entab(BufReader::new(file)),
Err(_) => println!("{} not found", &args[1])
}
}
}
fn split_at(&self, mid: usize) -> (&str, &str)
ファイルを「ランレングス」で符号化するプログラムと復号するプログラムを作ります。
ランレングスとは「連続して現れるものの長さ」という意味で、データ内で同じ値が並んでいる場合はその値と個数で符号化する方法のことを、「ランレングス圧縮」または「ランレングス符号化」といいます。
ランレングスはとても簡単な符号化方式ですが、それでもいくつかの方法が考えられます。いちばん簡単な方法は、データの値とデータの個数で表す方法です。たとえば、次の文字列を考えてみましょう。
文字列 abccddeeeeffffgggggggghhhhhhhh (30)
同じ文字が連続して並んでいますね。これを文字と個数で表すと、次のようになります。
=>1b1c2d2e4f4g8h8 (16)
元データ 30 文字を 16 文字で表すことができます。また、復号も簡単にできることはすぐにわかると思います。このように、ランレングスはとても単純な方法ですが、画像データには大きな効果を発揮する場合があります。たとえば、白黒の画像ではデータが 2 種類しかないので、最初のデータが白か黒のどちらであるか区別できれば、あとは個数だけの情報でデータを圧縮することができます。
なお、ランレングス符号にはいろいろな種類があります。興味のある方は拙作のページ「Algorithms with Python: 連長圧縮」をお読みくださいませ。
リスト : ランレングス符号化 (rle.rs)
use std::fs::File;
use std::io::prelude::*;
use std::io::{BufReader, BufWriter};
fn rle_file(src: &String, dst: &String) -> std::io::Result<()> {
let rd = BufReader::new(File::open(src)?);
let mut wr = BufWriter::new(File::create(dst)?);
let mut n = 0;
let mut c = 0;
for r in rd.bytes() {
let x = r?;
if n == 0 {
c = x;
n = 1;
} else if c == x {
n += 1;
if n == 256 {
wr.write(&[c, 255u8])?;
n = 0;
}
} else {
wr.write(&[c, (n - 1) as u8])?;
c = x;
n = 1;
}
}
if n > 0 {
wr.write(&[c, (n - 1) as u8])?;
}
Ok(())
}
fn main() {
let args: Vec<_> = std::env::args().collect();
if args.len() < 3 {
println!("usage: rle input_file output_file");
} else {
match rle_file(&args[1], &args[2]) {
Ok(_) => (),
Err(err) => println!("{}", err.to_string())
}
}
}
リスト : ランレングス復号 (rld.rs)
use std::fs::File;
use std::io::prelude::*;
use std::io::{BufReader, BufWriter};
fn rld_file(src: &String, dst: &String) -> std::io::Result<()> {
let mut rd = BufReader::new(File::open(src)?).bytes();
let mut wr = BufWriter::new(File::create(dst)?);
while let Some(c) = rd.next() {
if let Some(n) = rd.next() {
let x = c?;
let m = n? as usize;
for _ in 0 .. m + 1 { wr.write(&[x])?; }
} else {
panic!("broken input_file");
}
}
Ok(())
}
fn main() {
let args: Vec<_> = std::env::args().collect();
if args.len() < 3 {
println!("usage: rld input_file output_file");
} else {
match rld_file(&args[1], &args[2]) {
Ok(_) => (),
Err(err) => println!("{}", err.to_string())
}
}
}
有理数 (分数) を計算する簡単なプログラムです。分子と分母には i64 を使っているので、四則演算でオーバーフローすることがあります。実用的なプログラムを作るのであれば、多倍長整数が必要になります。なお、Rust のクレート num には多倍長整数 BigInt や有理数 BigRational が用意されているので、私たちが有理数を自作する必要はありません。Rust のお勉強ということで、あえて簡単なプログラムを作ってみましょう。
リスト : 簡単な有理数 (ratio.rs)
use std::ops::{Add, Div, Mul, Sub};
use std::fmt;
#[derive(PartialEq, Debug, Copy, Clone)]
struct Ratio {
numerator: i64, // 分子
denominator: i64 // 分母
}
// ユークリッドの互除法
fn gcd(a: i64, b: i64) -> i64 {
if b == 0 {
a
} else {
gcd(b, a % b)
}
}
// 有理数の表示
// トレイト std::fmt::Display を実装すると、
// 書式文字列の {} で Ratio を表示することができる
impl fmt::Display for Ratio {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
if self.denominator == 1 {
write!(f, "{}", self.numerator)
} else {
write!(f, "{}/{}", self.numerator, self.denominator)
}
}
}
// メソッドの実装
impl Ratio {
fn new(p: i64, q: i64) -> Ratio {
if q == 0 { panic!("Ratio: divide by zero"); }
let g = gcd(p.abs(), q.abs());
let s = if q < 0 { -1 } else { 1 };
Ratio { numerator: s * p / g, denominator: s * q / g}
}
fn from_integer(n: i64) -> Ratio {
Ratio { numerator: n, denominator: 1 }
}
fn to_integer(&self) -> i64 {
self.numerator / self.denominator
}
fn to_float(&self) -> f64 {
self.numerator as f64 / self.denominator as f64
}
fn numer(&self) -> i64 { self.numerator }
fn denom(&self) -> i64 { self.denominator }
fn is_integer(&self) -> bool {
self.denominator == 1
}
}
// 演算子の多重定義
impl Add for Ratio {
type Output = Ratio;
fn add(self, other: Ratio) -> Ratio {
let p = self.numerator * other.denominator + other.numerator * self.denominator;
let q = self.denominator * other.denominator;
Ratio::new(p, q)
}
}
impl Sub for Ratio {
type Output = Ratio;
fn sub(self, other: Ratio) -> Ratio {
let p = self.numerator * other.denominator - other.numerator * self.denominator;
let q = self.denominator * other.denominator;
Ratio::new(p, q)
}
}
impl Mul for Ratio {
type Output = Ratio;
fn mul(self, other: Ratio) -> Ratio {
let p = self.numerator * other.numerator;
let q = self.denominator * other.denominator;
Ratio::new(p, q)
}
}
impl Div for Ratio {
type Output = Ratio;
fn div(self, other: Ratio) -> Ratio {
let p = self.numerator * other.denominator;
let q = self.denominator * other.numerator;
Ratio::new(p, q)
}
}
// 簡単なテスト
fn main() {
let a = Ratio::new(1, 2);
println!("{}", a);
println!("{}", a.numer());
println!("{}", a.denom());
println!("{}", a.is_integer());
println!("{}", a.to_integer());
println!("{}", a.to_float());
let b = Ratio::from_integer(7);
println!("{}", b);
println!("{}", b.numer());
println!("{}", b.denom());
println!("{}", b.is_integer());
println!("{}", b.to_integer());
println!("{}", b.to_float());
let c = Ratio::new(1, 7);
println!("{}", c + a);
println!("{}", c - a);
println!("{}", c * a);
println!("{}", c / a);
// 1 + 1/2 + 1/4 + 1/8 + ... + 1/2n + ... = 2
println!("----- 2 に収束する級数 -----");
let mut d = Ratio::from_integer(1);
let mut n = 2;
for _ in 0 .. 20 {
d = d + Ratio::new(1, n);
println!("{}, {}", d, d.to_float());
n *= 2;
}
}
1/2 1 2 false 0 0.5 7 7 1 true 7 7 9/14 -5/14 1/14 2/7 ----- 2 に収束する級数 ----- 3/2, 1.5 7/4, 1.75 15/8, 1.875 31/16, 1.9375 63/32, 1.96875 127/64, 1.984375 255/128, 1.9921875 511/256, 1.99609375 1023/512, 1.998046875 2047/1024, 1.9990234375 4095/2048, 1.99951171875 8191/4096, 1.999755859375 16383/8192, 1.9998779296875 32767/16384, 1.99993896484375 65535/32768, 1.999969482421875 131071/65536, 1.9999847412109375 262143/131072, 1.9999923706054688 524287/262144, 1.9999961853027344 1048575/524288, 1.9999980926513672 2097151/1048576, 1.9999990463256836
Ratio を使って簡単なパズルを解いてみましょう。
下図の A から I の場所に 1 から 9 までの数字をひとつずつ配置します。3 つの分数を足すと 1 / N になる配置を求めてください。ただし、N は 2, 3, 4, 6, 10 とします。
このパズルの元ネタは N = 1 の場合で、参考文献『超々難問数理パズル 解けるものなら解いてごらん』に掲載されています。ちなみに、3 つの分数の和が整数になる場合、その値は 1 しかありません。また、値が 1 / N (N は整数) になる場合は 2, 3, 4, 6, 10 の 5 通りです。
リスト : 小町分数
//
// Ratio の定義は省略
//
// 順列の生成
fn permutations(n: i64, m: i64, func: fn(&Vec<i64>) -> ()) {
fn perm(n: i64, m: i64, func: fn(&Vec<i64>) -> (), xs: &mut Vec<i64>) {
if m == 0 {
func(xs);
} else {
for x in 1 .. n + 1 {
if !xs.contains(&x) {
xs.push(x);
perm(n, m - 1, func, xs);
xs.pop();
}
}
}
}
perm(n, m, func, &mut vec![]);
}
fn komachi_solver(xs: &Vec<i64>) {
let a = Ratio::new(xs[0], xs[1] * 10 + xs[2]);
let b = Ratio::new(xs[3], xs[4] * 10 + xs[5]);
let c = Ratio::new(xs[6], xs[7] * 10 + xs[8]);
let d = a + b + c;
if d.numer() == 1 && xs[0] < xs[3] && xs[3] < xs[6] {
println!("{}/{}{} + {}/{}{} + {}/{}{} = {}",
xs[0],xs[1],xs[2],xs[3],xs[4],xs[5],xs[6],xs[7],xs[8], d);
}
}
fn main() {
println!("----- 小町分数 -----");
permutations(9, 9, komachi_solver);
}
----- 小町分数 ----- 1/24 + 3/56 + 7/98 = 1/6 1/26 + 5/39 + 7/84 = 1/4 1/32 + 5/96 + 7/84 = 1/6 1/38 + 2/95 + 4/76 = 1/10 1/48 + 5/32 + 7/96 = 1/4 1/56 + 3/72 + 9/84 = 1/6 1/96 + 5/32 + 7/84 = 1/4 1/96 + 5/48 + 7/32 = 1/3 2/18 + 5/63 + 7/49 = 1/3 2/19 + 4/57 + 6/38 = 1/3 3/27 + 6/54 + 9/81 = 1/3 3/48 + 5/16 + 9/72 = 1/2 3/54 + 6/72 + 9/81 = 1/4 5/34 + 7/68 + 9/12 = 1
分子が 1 の分数を「単位分数」といいます。単位分数の和は古代エジプト人がよく研究していたそうです。参考文献『Ｃ言語による最新アルゴリズム事典』に「分数を単位分数の和で表す方法」がありましたので紹介します。
0 < m / n < 1 の分数を単位分数の和で表します。まず、n / m の商 q を求めます。もし、割り切れれば単位分数になりますね。そうでなければ、m / n から 1 / (q + 1) を引き算して M / N を求めます。あとは、M / N に対して同じ処理を繰り返すだけです。次の式を見てください。
M / N = m / n - 1 / (q + 1) M / N = (m(q + 1) - n) / n(q + 1) M = m(q + 1) - n = m - (n - mq) = m - (n mod m)
0 < (n mod m) < m ですから、M は必ず m より小さな値になります。つまり、この処理を繰り返していくと m は必ず 1 になるので、分数を単位分数の和で表すことができる、というわけです。なるほど納得のアルゴリズムですね。たとえば、11 / 13 を単位分数の和で表してみましょう。
11 / 13 => q = 1, 11 / 13 - 1 / 2 = 9 / 26 9 / 26 => q = 2, 9 / 26 - 1 / 3 = 1 / 78 11 / 13 = 1 / 2 + 1 / 3 + 1 / 78
このように、分子 m の値は減少していきます。このアルゴリズムを Rust でプログラムすると、次のようになります。
リスト : 単位分数の和
//
// Ratio の定義は省略
//
// 分数を単位分数の和で表す
fn egy_frac(mut m: i64, mut n: i64) -> Vec<Ratio> {
let mut result: Vec<Ratio> = vec![];
while n % m != 0 {
let q = n / m + 1;
result.push(Ratio::new(1, q));
m = m * q - n;
n *= q;
}
result.push(Ratio::new(1, n / m));
result
}
// 単位分数の和を表示
fn print_egy_frac(buff: &[Ratio]) {
let mut a = Ratio::from_integer(0);
for x in buff {
print!("{} ", x);
a = a + *x;
}
println!("=> {}", a);
}
fn main() {
println!("----- 単位分数の和 -----");
print_egy_frac(&egy_frac(11, 13));
print_egy_frac(&egy_frac(12, 13));
print_egy_frac(&egy_frac(19, 23));
}
1/2 1/3 1/78 => 11/13 1/2 1/3 1/12 1/156 => 12/13 1/2 1/4 1/14 1/215 1/138460 => 19/23
実数 a の平方根 \(\sqrt a\) の値を求める場合、方程式 \(x^2 - a = 0\) を Newton (ニュートン) 法で解くことが多いと思います。方程式を \(f(x)\), その導関数を \(f'(x)\) とすると、ニュートン法は次の漸化式の値が収束するまで繰り返す方法です。
平方根を求める場合、導関数は f'(x) = 2x になるので、漸化式は次のようになります。
参考文献『Ｃ言語による最新アルゴリズム事典』によると、\(\sqrt a\) より大きめの初期値から出発し、置き換え x <- (x + a / x) / 2 を減少が止まるまで繰り返すことで \(\sqrt a\) の正確な値を求めることができるそうです。
Rust でプログラムすると、次のようになります。
リスト : ニュートン法により平方根を求める
fn sqrt(x: f64) -> f64 {
// √x よりも大きめの初期値を求める
fn init(mut x: f64) -> f64 {
let mut s = 1.0;
while s < x {
s *= 2.0;
x /= 2.0;
}
s
}
if x < 0.0 { panic!("sqrt: domain error"); }
let mut p = if x > 1.0 { init(x) } else { 1.0 };
loop {
let q = (p + x / p) / 2.0;
if q >= p { break; }
p = q;
}
p
}
fn main() {
for n in 2 .. 10 {
println!("{}, {}", sqrt(n as f64), (n as f64).sqrt());
}
}
1.414213562373095, 1.4142135623730951 1.7320508075688772, 1.7320508075688772 2, 2 2.23606797749979, 2.23606797749979 2.449489742783178, 2.449489742783178 2.6457513110645907, 2.6457513110645907 2.82842712474619, 2.8284271247461903 3, 3
平方根の整数部もニュートン法を使って求めることができますが、次の公式を使って平方根の整数部分を求めることもできます。
式 (1) は、奇数 \(1\) から \(2n - 1\) の総和は \(n^2\) になることを表しています。式 (2) のように、整数 m の値が \(n^2\) より大きくて \((n + 1)^2\) より小さいのであれば、m の平方根の整数部分は n であることがわかります。これは m から奇数 \(1, 3, 5, \ldots, (2n - 1), (2n + 1)\) を順番に引き算していき、引き算できなくなった時点の (2n + 1) / 2 = n が m の平方根になります。参考文献『平方根計算法 (PDF)』によると、この方法を「めのこ平方」と呼ぶそうです。
プログラムは次のようになります。
リスト : めのこ平方
fn isqrt(mut n: i64) -> i64 {
let mut m = 1;
while n >= m {
n -= m;
m += 2;
}
m / 2
}
この方法はとても簡単ですが、数が大きくなると時間がかかるようになります。そこで、整数を 2 桁ずつ分けて計算していくことにします。次の図を見てください。
整数 6789 を 67 と 89 に分ける 1 + 3 + ... + 15 = 82 < 67 両辺を 100 倍すると 802 < 6700 < 6789 802 = 1 + 3 + ... + 159 (= 2 * 80 - 1) 161 + 163 < (6789 - 6400 = 389) < 161 + 163 + 165
整数 6789 を 67 と 89 に分けます。最初に 67 の平方根を求めます。この場合は 8 になり、82 < 67 を満たします。次に、この式を 100 倍します。すると、802 < 6700 になり、6700 に 89 を足した 6789 も 802 より大きくなります。802 は 1 から 159 までの奇数の総和であることはすぐにわかるので、6789 - 6400 = 389 から奇数 161, 163, ... を順番に引き算していけば 6789 の平方根を求めることができます。
プログラムは次のようになります。
リスト : めのこ平方 (2)
// isqrt の定義は省略
// 高速版
fn isqrt_fast(n: i64) -> i64 {
fn isqrt(mut n: i64, mut m: i64) -> i64 {
while n >= m {
n -= m;
m += 2;
}
m / 2
}
if n < 100 {
isqrt(n, 1)
} else {
let m = 10 * isqrt_fast(n / 100);
isqrt(n - m * m, 2 * m + 1)
}
}
fn main() {
println!("{}", isqrt(4));
println!("{}", isqrt(16));
println!("{}", isqrt(64));
println!("{}", isqrt(80));
println!("{}", isqrt(81));
println!("{}", isqrt(82));
println!("{}", isqrt(100));
println!("{}", isqrt_fast(6789));
println!("{}", isqrt_fast(123456789));
println!("{}", isqrt_fast(1234567890));
}
2 4 8 8 9 9 10 82 11111 35136
配列の「べき集合」を求める高階関数 power_set() を作ります。たとえば配列 [1, 2, 3] のべき集合は [], [1], [2], [3], [1, 2], [1, 3], [2, 3], [1, 2, 3] になります。
リスト : べき集合
fn power_set<T: Copy>(func: fn(&[T]) -> (), table: &[T]) {
fn iter<T: Copy>(func: fn(&[T]) -> (), table: &[T], buff: &mut Vec<T>) {
if table.len() == 0 {
func(buff);
} else {
iter(func, &table[1..], buff);
buff.push(table[0]);
iter(func, &table[1..], buff);
buff.pop();
}
}
iter(func, table, &mut vec![])
}
fn main() {
fn print_vec(xs: &[i32]) {
println!("{:?}", xs);
}
power_set(print_vec, &[1,2,3,4]);
}
power_set() は簡単です。実際の処理は局所関数 iter() で行います。table の先頭要素を選択する場合は、その要素をベクタ buff に追加して iter() を再帰呼び出しします。選択しない場合は、ベクタ buff に要素を追加せずに iter() を再帰呼び出しするだけです。これでべき集合の要素をすべて求めることができます。
[] [4] [3] [3, 4] [2] [2, 4] [2, 3] [2, 3, 4] [1] [1, 4] [1, 3] [1, 3, 4] [1, 2] [1, 2, 4] [1, 2, 3] [1, 2, 3, 4]
m 個の整数 1, 2, ..., m の順列を考えます。このとき、n 番目 (先頭要素が 1 番目) の要素が整数 n ではない順列を「完全順列 (derangement)」といいます。今回は 1 から m までの整数値で完全順列を生成する関数を作ります。
リスト : 完全順列
fn derangement(m: i32, func: fn(&[i32]) -> ()) {
fn perm(m: i32, func: fn(&[i32]) -> (), xs: &mut Vec<i32>) {
let n = xs.len() as i32;
if n == m {
func(xs);
} else {
for x in 1 .. m + 1 {
if n + 1 != x && !xs.contains(&x) {
xs.push(x);
perm(m, func, xs);
xs.pop();
}
}
}
}
perm(m, func, &mut vec![]);
}
fn main() {
fn print_vec(xs: &[i32]) {
println!("{:?}", xs);
}
derangement(3, print_vec);
derangement(4, print_vec);
derangement(5, print_vec);
}
実際の処理は局所関数 perm() で行います。基本的には 1 から m までの数字を m 個選ぶ順列を生成する処理と同じです。n 番目の数字を選ぶとき、数字 x が n + 1 と等しい場合は x を選択しません。n が m と等しい場合は m 個の数字を選んだので xs の内容を表示します。これで完全順列を生成することができます。
[2, 3, 1] [3, 1, 2] [2, 1, 4, 3] [2, 3, 4, 1] [2, 4, 1, 3] [3, 1, 4, 2] [3, 4, 1, 2] [3, 4, 2, 1] [4, 1, 2, 3] [4, 3, 1, 2] [4, 3, 2, 1] [2, 1, 4, 5, 3] [2, 1, 5, 3, 4] [2, 3, 1, 5, 4] [2, 3, 4, 5, 1] [2, 3, 5, 1, 4] [2, 4, 1, 5, 3] [2, 4, 5, 1, 3] [2, 4, 5, 3, 1] [2, 5, 1, 3, 4] [2, 5, 4, 1, 3] [2, 5, 4, 3, 1] [3, 1, 2, 5, 4] [3, 1, 4, 5, 2] [3, 1, 5, 2, 4] [3, 4, 1, 5, 2] [3, 4, 2, 5, 1] [3, 4, 5, 1, 2] [3, 4, 5, 2, 1] [3, 5, 1, 2, 4] [3, 5, 2, 1, 4] [3, 5, 4, 1, 2] [3, 5, 4, 2, 1] [4, 1, 2, 5, 3] [4, 1, 5, 2, 3] [4, 1, 5, 3, 2] [4, 3, 1, 5, 2] [4, 3, 2, 5, 1] [4, 3, 5, 1, 2] [4, 3, 5, 2, 1] [4, 5, 1, 2, 3] [4, 5, 1, 3, 2] [4, 5, 2, 1, 3] [4, 5, 2, 3, 1] [5, 1, 2, 3, 4] [5, 1, 4, 2, 3] [5, 1, 4, 3, 2] [5, 3, 1, 2, 4] [5, 3, 2, 1, 4] [5, 3, 4, 1, 2] [5, 3, 4, 2, 1] [5, 4, 1, 2, 3] [5, 4, 1, 3, 2] [5, 4, 2, 1, 3] [5, 4, 2, 3, 1]
完全順列の総数を「モンモール数 (Montmort number)」といいます。モンモール数は次の漸化式で求めることができます。
モンモール数を求める関数 montmort_number() を作ります。
リスト : モンモール数
fn montmort_number(n: i64) -> i64 {
if n <= 1 {
0
} else if n == 2 {
1
} else {
(n - 1) * (montmort_number(n - 1) + montmort_number(n - 2))
}
}
// 別解
fn montmort_number1(n: i64) -> i64 {
let mut a = 0;
let mut b = 1;
for i in 1 .. n {
let c = (i + 1) * (a + b);
a = b;
b = c;
}
a
}
fn main() {
for n in 1 .. 20 {
println!("{}, {}", n, montmort_number(n));
// println!("{}, {}", n, montmort_number1(n));
}
}
1, 0 2, 1 3, 2 4, 9 5, 44 6, 265 7, 1854 8, 14833 9, 133496 10, 1334961 11, 14684570 12, 176214841 13, 2290792932 14, 32071101049 15, 481066515734 16, 7697064251745 17, 130850092279664 18, 2355301661033953 19, 44750731559645106
関数 montmort_number() は公式をそのままプログラムしただけです。二重再帰になっているので、実行速度はとても遅くなります。これを繰り返しに変換すると別解のようになります。考え方はフィボナッチ数列と同じです。変数 a に i 番目の値を、b に i + 1 番目の値を保存しておきます。すると、i + 2 番目の値は (i + 1) * (a + b) で計算することができます。あとは、b の値を a に、新しい値を b にセットして処理を繰り返すだけです。
なお、i64 では 19 番目のモンモール数までしか求めることができません。大きなモンモール数を求めたい場合は多倍長整数を使ってください。
カッコ列は ( と ) からなる列のことで、バランスが取れているカッコ列は、右カッコで閉じることができる、つまり右カッコに対応する左カッコがある状態のことをいいます。たとえば n = 1 の場合、( ) はバランスの取れたカッコ列ですが、) ( はバランスが取れていません。今回はカッコ列を生成する高階関数 kakko() を作ります。
リスト : カッコ列
fn kakko(func: fn(&str) -> (), n: i32) {
fn iter(x: i32, y: i32, n: i32, func: fn(&str) -> (), a: &mut String) {
if x == y && y == n {
func(&a)
} else {
if x < n {
a.push('(');
iter(x + 1, y, n, func, a);
a.pop();
}
if y < x {
a.push(')');
iter(x, y + 1, n, func, a);
a.pop();
}
}
}
iter(0, 0, n, func, &mut String::new());
}
fn main() {
fn print_str(s: &str) {
println!("{}", s);
}
kakko(print_str, 3);
kakko(print_str, 4);
}
((())) (()()) (())() ()(()) ()()() (((()))) ((()())) ((())()) ((()))() (()(())) (()()()) (()())() (())(()) (())()() ()((())) ()(()()) ()(())() ()()(()) ()()()()
カッコ列の生成は簡単です。局所関数 iter() の引数 x が左カッコの個数、引数 y が右カッコの個数を表します。引数 a は累積変数でカッコ列を表す文字列 (String 型) です。
バランスの取れたカッコ列の場合、x, y, n には y <= x && x <= n の関係が成り立ちます。x == y && y == n の場合、カッコ列がひとつ完成しました。引数の関数 func を呼び出します。そうでなければ、iter() を再帰呼び出しします。x < n であれば左カッコを追加し、y < x であれば右カッコを追加します。String のメソッド push() は文字列の末尾に文字を追加します。pop() は文字列の末尾から文字を取り除きます。これでカッコ列を生成することができます。
カタラン数 - Wikipedia によると、カッコ列の総数は「カタラン数 (Catalan number)」になるとのことです。カタラン数 \(\mathrm{C}_n\) は次に示す公式で求めることができます。
これをそのままプログラムしてもいいのですが、それではちょっと面白くないので別な方法でプログラムを作ってみましょう。カタラン数は次に示す経路図において、A から B までの最短距離の道順を求めるとき、対角線を超えないものの総数に一致します。
Ｂ Ｂ
┌─┬─┬─┬─┐ ┌─┬─┬─０─14
│ │ │ │ │ │ │ │ │ │
├─┼─┼─┼─┤ ├─┼─０─５─14
│ │ │ │ │ │ │ │ │ │
├─┼─┼─┼─┤ ├─０─２─５─９
│ │ │ │ │ │ │ │ │ │
├─┼─┼─┼─┤ ０─１─２─３─４
│ │ │ │ │ │ │ │ │ │
└─┴─┴─┴─┘ １─１─１─１─１
Ａ Ａ
図 : 道順の総数の求め方
A からある地点にいたる最短距離の道順の総数は、左隣と真下の地点の値を足したものになります。一番下の地点は 1 で、対角線を越えた地点は 0 になります。あとは下から順番に足し算していけば、A から B までの道順の総数を求めることができます。上図の場合はカラタン数 Ｃ4 に相当し、その値は 14 となります。
プログラムは配列を使うと簡単です。次の図を見てください。
0 : [1, 1, 1, 1, 1] 1 : [1, 1, 1, 1, 1,] 2 : [1, 1, 1+1=2, 2+1=3, 3+1=4] => [1, 1, 2, 3, 4] 3 : [1, 1, 2, 3+2=5, 5+4=9] => [1, 1, 2, 5, 9] 4 : [1, 1, 2, 5, 5+9=14] => [1, 1, 2, 5, 14]
上図は Ｃn (n = 4) を求める場合です。大きさが n + 1, 要素の値が 1 のベクタを用意します。n = 0, 1 の場合は n 番目の要素をそのまま返します。n が 2 よりも大きい場合、変数 i を 2 に初期化して、i - 1 番目以降の要素の累積和を求めます。
たとえば i = 2 の場合、2 番目の要素は 1 番目の要素と自分自身を加算した値 2 になります。3 番目の要素は 2 番目の要素と自分自身を足した値 3 になり、4 番目の要素は 3 + 1 = 4 になります。次に i を +1 して同じことを繰り返します。3 番目の要素は 2 + 3 = 5 になり、4 番目の要素は 5 + 4 = 9 になります。i = 4 のとき、4 番目の要素は 5 + 9 = 14 となり、Ｃ4 の値を求めることができました。
リスト : カッコ列の総数
// カタラン数
fn catalan_number(n: usize) -> usize {
let mut table: Vec<usize> = vec![1; n + 1];
for i in 2 .. n + 1 {
for j in i .. n + 1 {
table[j] += table[j - 1];
}
}
table[n]
}
fn main() {
for n in 1 .. 20 {
println!("{}, {}", n, catalan_number(n));
}
}
1, 1 2, 2 3, 5 4, 14 5, 42 6, 132 7, 429 8, 1430 9, 4862 10, 16796 11, 58786 12, 208012 13, 742900 14, 2674440 15, 9694845 16, 35357670 17, 129644790 18, 477638700 19, 1767263190
リスト : たらいまわし関数
use std::time::Instant;
fn tarai(x: i32, y: i32, z: i32) -> i32 {
if x <= y {
y
} else {
tarai(tarai(x - 1, y, z), tarai(y - 1, z, x), tarai(z - 1, x, y))
}
}
fn tak(x: i32, y: i32, z: i32) -> i32 {
if x <= y {
z
} else {
tak(tak(x - 1, y, z), tak(y - 1, z, x), tak(z - 1, x, y))
}
}
fn main() {
{
let start = Instant::now();
println!("{}", tarai(14, 7, 0));
let end = start.elapsed();
println!("{}.{:03}秒経過しました。", end.as_secs(), end.subsec_nanos() / 1_000_000);
}
{
let start = Instant::now();
println!("{}", tak(22, 11, 0));
let end = start.elapsed();
println!("{}.{:03}秒経過しました。", end.as_secs(), end.subsec_nanos() / 1_000_000);
}
}
関数 tarai() や tak() は「たらいまわし関数」といい、再帰的に定義されています。これらの関数は、引数の与え方によっては実行に時間がかかるため、Lisp などのベンチマークに利用されることがあります。tarai() は通称「竹内関数」と呼ばれていて、日本の代表的な Lisper である竹内郁雄先生によって考案されたそうです。そして、tak() は tarai() のバリエーションで、John Macarthy 先生によって作成されたそうです。
それでは、さっそく実行してみましょう。
14 1.059秒経過しました。 11 1.169秒経過しました。 実行環境 : Rust 1.83.0, Ubunts 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
このように、たらいまわし関数は引数の値が小さくても実行に時間がかかります。
tarai() は「遅延評価」を行う処理系、たとえば関数型言語の Haskell では高速に実行することができます。Scheme でも delay と force を使って遅延評価を行うことができます。tarai() のプログラムを見てください。x <= y のときに y を返しますが、このとき引数 z の値は必要ありませんね。引数 z の値は x > y のときに計算するようにすれば、無駄な計算を省略することができます。なお、tak() は x <= y のときに z を返しているため、遅延評価で高速化することはできません。ご注意ください。
クロージャをサポートしている言語の場合、完全ではありませんがクロージャを使って遅延評価を行うことができます。Rust にもクロージャがあるので、次のようにプログラムしたのですが、コンパイルエラーになりました。
リスト : 遅延評価による高速化 (コンパイルエラー)
fn tarai_lazy<F>(x: i32, y: i32, z: F) -> i32 where F: Fn() -> i32 {
if x <= y {
y
} else {
let zz = z();
tarai_lazy(tarai_lazy(x - 1, y, z),
tarai_lazy(y - 1, zz, || x),
|| tarai_lazy(zz - 1, x, || y))
}
}
reached the recursion limit while instantiating `tarai_lazy::<[closure@tarai.rs:25:42: 25:46]>`
再帰関数にクロージャを渡す場合、Rust には制限があるようです。具体的には、再帰関数の中でクロージャを生成し、それを再帰関数に渡すことはできないようです。簡単な例を示しましょう。
リスト : 整数の和 (正常に動作)
fn sum<F>(n: i32, a: i32, func: F) where F: Fn(i32) -> () {
if n <= 0 {
func(a);
} else {
sum(n - 1, a + n, func);
}
}
fn main() {
sum(1000, 0, |x| println!("{}", x));
}
500500
関数 sum は引数 func にクロージャを受け取ります。sum の中で再帰呼び出しをしていますが、新しいクロージャは生成せず、func をそのまま再帰呼び出しの引数に渡しています。これは正常にコンパイルすることができます。ただし、func を再帰呼び出しの sum() に渡すとき、所有権が移動することに注意してください。二重再帰のとき、このことを忘れるとコンパイルエラーになります。
ところが、次のように sum の中でクロージャを生成し、それを sum の引数に渡すと、先ほどと同じコンパイルエラーになります。
リスト : 整数の和 (コンパイルエラー)
fn sum<F>(n: i32, a: i32, func: F) where F: Fn(i32) -> () {
if n <= 0 {
func(a);
} else {
sum(n - 1, a + n, |x| println!("{}", x));
}
}
fn main() {
sum(1000, 0, |x| println!("{}", x));
}
reached the recursion limit while instantiating `sum::<[closure@sum2.rs:5:27: 5:48]>`
この場合、トレイトオブジェクトを使うと上手くいきます。
リスト : 整数の和 (トレイトオブジェクト)
fn sum(n: i32, a: i32, func: &dyn Fn(i32) -> ()) {
if n <= 0 {
func(a);
} else {
sum(n - 1, a + n, &|x| println!("{}", x));
}
}
fn main() {
sum(1000, 0, &|x| println!("{}", x));
}
500500
トレイトオブジェクトを使うと、たらいまわし関数は次のようになります。
リスト : 遅延評価による高速化 (トレイトオブジェクト)
fn tarai_lazy(x: i32, y: i32, z: &dyn Fn() -> i32) -> i32 {
if x <= y {
y
} else {
let zz = z();
tarai_lazy(tarai_lazy(x - 1, y, z),
tarai_lazy(y - 1, zz, &|| x),
&|| tarai_lazy(zz - 1, x, &|| y))
}
}
実際に試してみると、計測不能なほど高速になりました。興味のある方はいろいろ試してみてください。
たらいまわし関数が遅いのは、同じ値を何度も計算しているためです。この場合、表 (table) を使って処理を高速化することができます。同じ値を何度も計算することがないように、計算した値は表に格納しておいて、2 回目以降は表から計算結果を求めるようにします。このような手法を「表計算法」とか「メモ化」といいます。
Rust の場合、HashMap を使うと簡単です。次のリストを見てください。
リスト : メモ化
fn tak_memo(x: i32, y: i32, z: i32, m: &mut HashMap<(i32, i32, i32), i32>) -> i32 {
let key = (x, y, z);
match m.get(&key) {
Some(&v) => v,
None => {
if x <= y {
m.insert(key, z);
z
} else {
let v = tak_memo(tak_memo(x - 1, y, z, m),
tak_memo(y - 1, z, x, m),
tak_memo(z - 1, x, y, m), m);
m.insert(key, v);
v
}
}
}
}
tak_memo() の値を格納する HashMap を引数 m に渡します。tak_memo() では、引数 x, y, z を要素とするタプルを作り、それをキーとしてHashMap を検索します。HashMap に key があれば、その値を返します。そうでなければ、値を計算して HashMap にセットして、その値を返します。実際に試してみると、計測不能なほど高速になりました。興味のある方はいろいろ試してみてください。