整数をインクリメントするスレッド vs. その整数を表示するスレッド

1. 動機

以下のプログラムを考える：
メインスレッドで共有の整数nを用意する（初期値は0）。そして以下の２つのスレッドを立ち上げる。
・1000回だけnを標準出力に書き出す。
・1000回だけnをインクリメントする。

前者のスレッドは出力行を書式設定して書き出さなければならない（つまり、後者のスレッドに比べて多くの機械命令を必要とする）ので、前者のスレッドが表示するnの値はすぐに1000になるはずである。それを手を動かして確認したくなった。言語はRustを使用。

2. ロック等を使わず各々のスレッドが同期なしで動くとき

タイトルの通りのものをRustで書いた。Rustみたいに言語レベルで並行プログラミングの安全性を保障してる言語だと、簡単な並行プログラムなら同期なしで書くほうがかえって面倒だ。

use std::thread;
use std::ptr::read;

fn main() {
    let mut n = 0;
    let p = UnsafePtr::new(&mut n);
    let x = thread::spawn(move || {
        for _ in 0..1000 {
            println!("The value of n is {}.",unsafe {read(p.0) });
        }
    });
    let y = thread::spawn(move || { 
        for _ in 0..1000 {
            p.inc();
        }
    });
    let _ = x.join();
    let _ = y.join();
}

#[derive(Copy,Clone)]
struct UnsafePtr(*mut usize);

impl UnsafePtr {
    fn new(n: &mut usize) -> Self {
        Self(n as *mut usize)
    }
    fn inc(self) {
        let old = unsafe {self.0.read()};
        unsafe {self.0.write(old+1)};
    }
}

unsafe impl Sync for UnsafePtr {}
unsafe impl Send for UnsafePtr {}

上記のプログラムを3回実行した。

実行結果その１

The value of n is 100.
The value of n is 1000.
The value of n is 1000.
The value of n is 1000.
The value of n is 1000.
The value of n is 1000.
.
.
.
The value of n is 1000.

実行結果その２

The value of n is 187.
The value of n is 972.
The value of n is 1000.
The value of n is 1000.
The value of n is 1000.
The value of n is 1000.
.
.
.
The value of n is 1000.

実行結果その３

The value of n is 0.
The value of n is 342.
The value of n is 380.
The value of n is 405.
The value of n is 429.
The value of n is 452.
The value of n is 474.
The value of n is 497.
The value of n is 520.
The value of n is 541.
The value of n is 563.
The value of n is 586.
The value of n is 607.
The value of n is 630.
The value of n is 653.
The value of n is 674.
The value of n is 696.
The value of n is 718.
The value of n is 739.
The value of n is 761.
The value of n is 782.
The value of n is 804.
The value of n is 826.
The value of n is 848.
The value of n is 870.
The value of n is 891.
The value of n is 913.
The value of n is 935.
The value of n is 957.
The value of n is 980.
The value of n is 1000.
The value of n is 1000.
The value of n is 1000.
.
.
.
The value of n is 1000.

上述の通り、0から1000まで均等に表示されるのではなく、すぐに1000までインクリメントされる。

3. ロックを使って各々のスレッドが排他的にnにアクセスして動くとき

今度はロックを使って各々のスレッドは排他的にnにアクセスするようにする。

use std::{sync::{Arc,Mutex}, thread};

fn main() {
    let n = Arc::new(Mutex::new(0));
    let a = n.clone();
    let b = n.clone();
    let x = thread::spawn(move || {
        for _ in 0..1000 {
            println!("The value of n is {}.",a.lock().unwrap());
        }
    });
    let y = thread::spawn(move || { 
        for _ in 0..1000 {
            *b.lock().unwrap() += 1;
        }
    });
    let _ = x.join();
    let _ = y.join();
}

上記のプログラムを3回実行した。

実行結果は・・・nの値が1000になるのは大体500行目以降となったので、それを載せると無駄に記事が長くなってしまうのでやらない。ともかく、nの増加の勾配は2. のときに比べて緩やかになった。これは、println!マクロを実行している間はスレッドxがロックを取得しているため、スレッドyはprintln!が完了するまでnにアクセスすることが出来ないからだと思う。

4. Haskellでやってみる

ついでにHaskellでもやってみた。Haskellで2.と同等のプログラムの書き方が分からない（この辺使ったらできるのかな）ので、3. と同等のプログラムを書いてみた。

import Control.Concurrent.Async
import Control.Monad
import Control.Concurrent
import Text.Printf

main :: IO ()
main = do
    m <- newMVar (0 :: Int)
    x <- async . replicateM_ 1000 $ do
        n <- takeMVar m
        printf "The value of n is %d.\n" n 
        putMVar m n
    y <- async . replicateM_ 1000 $ do
        n <- takeMVar m
        putMVar m (n+1)
    wait x
    wait y

実行結果はこちら

The value of n is 0.
The value of n is 1.
The value of n is 2.
The value of n is 3.
The value of n is 4.
The value of n is 5.
The value of n is 6.
The value of n is 7.
The value of n is 8.
The value of n is 9.
The value of n is 10.
The value of n is 11.
The value of n is 12.
The value of n is 13.
.
.
.
The value of n is 994.
The value of n is 995.
The value of n is 996.
The value of n is 997.
The value of n is 998.
The value of n is 999.

あれ・・・ものすごく綺麗にインクリメントしてる（交互にスレッドが実行されている）。10回ほど実行しても毎回上記の結果となった。理由がわかったら追記します。

Rustのbreakとラベル付きループのメモ

たまに使いたくなるたびに調べているのでメモ

fn main() {
    // loopは式なのでbreakで値を返すことができる. forやwhileではできない. 
    let n = {
        let mut n = 3;
        loop {
            if 2u32.pow(n)+ n >= n.pow(3)-n {
                break n;
            }
            n += 1;
        }
    };
    println!("The minimum integer n >= 3 satisfying 2^n +n >= n^3-n is {}.",n);

    // loop, while, forにはラベルを付けてbreak, continueすることができる. 
    let mut x = (0,0,0);
    'main: for i in  2.. {
        'sub: for j in 2..300 {
            for k in 2..300 {
                if 2*j+k <= 200 {
                    continue 'sub;
                }
                if i+j+k >= 1000 {
                    x = (i,j,k);
                    break 'main;
                }
            }
        }
    }
    println!("An example of (i,j,k) satisfying 2*j+k > 200 and i+j+k <= 1000 is ({},{},{}).",x.0,x.1,x.2)
}

実行結果

The minimum integer n >= 3 satisfying 2^n +n >= n^3-n is 10.
An example of (i,j,k) satisfying 2*j+k > 200 and i+j+k <= 1000 is (402,299,299).

Rustで多倍長整数と素数生成

暗号のお勉強をしているときにRustで何かしらプログラムを書こうと思ったので、Rustでの多倍長整数の取扱についての軽いメモ。以下のソースコードはNバイトの素数をMiller–Rabinの素数判定法を用いて確率的に生成するもの（まともにデバッグしていないので間違えているかも）。多倍長整数演算はnum_bigintを使用。

use num_bigint_dig::BigUint;
use rand::{RngCore, rngs::ThreadRng};
use competitive::prime::*;
use num_traits::identities::Zero;
use num_traits::One;
use num_traits::Pow;

pub fn gen_prime<const N: usize>(r: usize) -> BigUint {
    let primes =  primes(1_000_000);
    let mut rng = ThreadRng::default();
    let mut vv = [0u8; N];
    rng.fill_bytes(&mut vv);
    vv[0] |= 1 << 7;
    vv[N-1] |= 1;
    let mut ret = BigUint::from_bytes_be(&vv);
    loop {
        if primes.iter().copied().all(|p| &ret % p != BigUint::zero()) {
            if miller_rabin::<N>(&ret, r, &mut rng) {
                return ret
            }
        }
        rng.fill_bytes(&mut vv);
        vv[0] |= 1 << 7;
        vv[N-1] |= 1;
        ret = BigUint::from_bytes_be(&vv);
    }
} 

fn miller_rabin<const N: usize>(n: &BigUint, r: usize, rng: &mut ThreadRng) -> bool {
    let n1 = n.clone()-1u32;
    let s = {
        let mut s = 0;
        for i in 1.. {
            let bb: BigUint = 2u32.into();
            if &n1 % bb.pow(i) != BigUint::zero() {
                s = i-1;
                break;
            }
        }
        s
    };
    let d = &n1/2u32.pow(s);
    let mut vv = [0u8; N];
    rng.fill_bytes(&mut vv);
    let mut a = BigUint::from_bytes_be(&vv);
    while &a >= &n1 {
        rng.fill_bytes(&mut vv);
        a = BigUint::from_bytes_be(&vv);
    }
    for _ in 0..r {
        let mut y = a.modpow(&d,&n);
        if y == BigUint::one() || &y == &n1 {
            continue;
        } else {
            for _ in 1..s {
                y = &y*&y;
                if &y == &n1 {
                    continue;
                } else {
                    return false
                }
            }
        }
        rng.fill_bytes(&mut vv);
        a = BigUint::from_bytes_be(&vv);
        while &a >= &n1 {
            rng.fill_bytes(&mut vv);
            a = BigUint::from_bytes_be(&vv);
        }
    }
    true
}

#[test]
fn test() {
    for _ in 0..5 {
        println!("{}",gen_prime::<64>(3))
    }
}

テストの結果:

running 1 test
8384442391023910089559381911244353428872487399318179273396015886507100438953585665678365906489416840722213922620210237392111730172861086168993311923836111
7117174289562103609173135325238483107034344612555783222340008053564899447108584024808296706859414511376712644654207288058426605121566231239674354322869571
7563023309084596682599574034479919223499374791081086371863943325349323174848426538472921938956966836116894178088405681548677619092250925604872140829006499
6845191684205106300741531169084809632281172184322263496188751375542125818145681444532252331962842985727123170223256606669735938766457801060667384315729663
11140042094397453071834468857251093913894980235722953128972718807724718935294988778661504441980549669968722894429691286957748143827788668127013507759578059
test prime::test ... ok

軽い説明:
まず、1,000,000未満の素数を生成する。続いて、Nバイトの配列をランダムに生成する（ThreadRngはCryptoRngを実装している）。このとき、最上位・最下位のビットは1とする。続いて、その配列を多倍長符号なし整数に変化する。あとは先程生成した1,000,000未満の素数で割り切れないかつMiller–Rabinの素数判定法でacceptされたら終了、そうでなければNバイトの配列をランダムに生成してやり直す。

このプログラムの素数生成ではp-1法などで簡単に因数分解される場合を考慮していないので、間違ってもRSA等の素数生成でこのプログラムを使ってはならない（そもそもまともにデバッグをしていない）