Yuki Nakata's Blog

天才 中ちゃん♪

カテゴリ: プログラミング

linux glibcのrand関数のソースコードを読んでみましょう。

どうやって乱数を生成しているのか調べてみましょう。全く怖くありません。

まず乱数についてですが、真性乱数と擬似乱数に分類されます。glibcの乱数は周期性があるため擬似乱数となります。優秀な擬似乱数としてはメルセンヌ・ツイスタ乱数があります。暗号に使われるのは真性乱数の方ですね。難しそうな説明はここまで。w

glibc ダウンロードディレクトリ glibc-2.22.tar.gz

rand関数は解凍してstdlibディレクトリのrandom.cとrandom_r.cに書かれています。

ソースコードを読む、というよりパクって実装して動かしてみると理解できます。

まず今は使われていない昔のrand関数を書いていきます。
oldrand.c

#include <stdio.h>
#include <stdlib.h>

/* x**31 + x**3 + 1.  */
#define    TYPE_3        3
#define    BREAK_3        128
#define    DEG_3        31
#define    SEP_3        3

//実はこの乱数テーブルをいじっているだけ 昔の乱数は2番目の-1726662223しか使っていない
static int32_t  randtbl[DEG_3 + 1] =
  {
    TYPE_3,

    -1726662223, 379960547, 1735697613, 1040273694, 1313901226,
    1627687941, -179304937, -2073333483, 1780058412, -1989503057,
    -615974602, 344556628, 939512070, -1249116260, 1507946756,
    -812545463, 154635395, 1388815473, -1926676823, 525320961,
    -1009028674, 968117788, -123449607, 1284210865, 435012392,
    -2017506339, -911064859, -370259173, 1132637927, 1398500161,
    -205601318,
  };

static struct random_data unsafe_state =
  {
    .fptr = &randtbl[SEP_3 + 1],
    .rptr = &randtbl[1],

    .state = &randtbl[1],

    .rand_type = TYPE_3,
    .rand_deg = DEG_3,
    .rand_sep = SEP_3,

    .end_ptr = &randtbl[sizeof (randtbl) / sizeof (randtbl[0])]
};

void oldrand(struct random_data *buf, int *result)
{
  int *state;

  state = buf->state;
  int val = state[0];
  val = ((state[0] * 1103515245) + 12345) & 0x7fffffff;  //randtbl[1]をごちゃごちゃいじっているだけ
  state[0] = val;
  *result = val;
}

int main()
{
  int i, result;

  for (i = 0; i < 10; i++) {
    oldrand(&unsafe_state, &result);
    printf("oldrand = %d\n", result);
  }

  return 0;
}

$ gcc oldrand.c -o oldrand
./oldrand
oldrand = 897822358
oldrand = 105331223
oldrand = 617672196
oldrand = 1888665581
oldrand = 1128537634
oldrand = 1434149043
oldrand = 980477552
oldrand = 1927835113
oldrand = 1772747374
oldrand = 1723946255

これが昔のrand関数です。randtblという乱数テーブルを用意していますが、たった一つしか使っていません。
val = ((state[0] * 1103515245) + 12345) & 0x7fffffff;
ここで数字をいじっているだけです。

次は現在使われている乱数
myrand.c

#include <stdio.h>
#include <stdlib.h>

/* x**31 + x**3 + 1.  */
#define    TYPE_0        0
#define    TYPE_3        3
#define    BREAK_3        128
#define    DEG_3        31
#define    SEP_3        3

//実はこの乱数テーブルをいじっているだけ ポインタを使って順送りしているだけ。最後まで行ったら最初に戻る
static int32_t randtbl[DEG_3 + 1] =
  {
    TYPE_3,

    -1726662223, 379960547, 1735697613, 1040273694, 1313901226,
    1627687941, -179304937, -2073333483, 1780058412, -1989503057,
    -615974602, 344556628, 939512070, -1249116260, 1507946756,
    -812545463, 154635395, 1388815473, -1926676823, 525320961,
    -1009028674, 968117788, -123449607, 1284210865, 435012392,
    -2017506339, -911064859, -370259173, 1132637927, 1398500161,
    -205601318,
  };

static struct random_data unsafe_state =
  {
    .fptr = &randtbl[SEP_3 + 1],
    .rptr = &randtbl[1],

    .state = &randtbl[1],

    .rand_type = TYPE_3,
    .rand_deg = DEG_3,
    .rand_sep = SEP_3,

    .end_ptr = &randtbl[sizeof (randtbl) / sizeof (randtbl[0])]
};

void myrand(struct random_data *buf, int *result)
{
  int *state;


      int32_t *fptr = buf->fptr;
      int32_t *rptr = buf->rptr;
      int32_t *end_ptr = buf->end_ptr;
      int32_t val;

      val = *fptr += *rptr;
      /* Chucking least random bit.  */
      *result = (val >> 1) & 0x7fffffff;  //ここで数字をいじる。他はポインタの書き換えでしかない
      ++fptr;
      if (fptr >= end_ptr) {
           fptr = state;
            ++rptr;
        } else {
           ++rptr;
          if (rptr >= end_ptr)
           rptr = state;
        }
      buf->fptr = fptr;
      buf->rptr = rptr;

}

int main()
{
  int i, result;

  for (i = 0; i < 10; i++) {
    myrand(&unsafe_state, &result);
    printf("oldrand = %d\n", result);
  }

  return 0;
}
$ gcc myrand.c -o myrand
$ ./myrand
oldrand = 1804289383
oldrand = 846930886
oldrand = 1681692777
oldrand = 1714636915
oldrand = 1957747793
oldrand = 424238335
oldrand = 719885386
oldrand = 1649760492
oldrand = 596516649
oldrand = 1189641421

難しそうですがやっていることはポインタを使って読み込むrandtblのアドレスを書き換えているだけです。

もちろん、これだけだと毎回同じ乱数が生成されます。

srandは与えられた種を元にrandtblのデータを書き換えているだけです。

それで毎回異なる乱数を生成することができるようになります。

仕組みは恐ろしく簡単です。

それよりも注意が必要なのは関数の実態がわかりづらいことです。

random.cの中でsrandが関連付けられています。

weak_alias (__srandom, srand)  とあります。

srandを呼び出すと__srandomに変換されますよ、という意味。

さらに__srandom_r関数へと変換されて、random_r.cファイルの中で定義されています。

つまりsrand 関数は内部では __srandom_r 関数になります。

仕組みはこんなもんです。この調子でglibcのソースコードを読んでいきましょう。

OpenMP 並列プログラミング2の続きです。

実際にOpenMPを使いこなしていきましょう。

グレゴリー・ライプニッツ級数でπを求めます。

π / 4 = arctan1 
 = (1 / 1) - (1 / 3) + (1 / 5) - (1 / 7) + (1 /9) - (1 / 11) ...無限に続く
 
グレゴリー・ライプニッツ級数は収束が遅いので実際にπを求める際には使われていません。

式がシンプルなforループで記述できるので並列プログラミング向きなのです。

並列化なし pi-leibniz.c
#include <stdio.h>
#include <time.h>

#define MAX 1000000000

int main()
{
int i;
double is = 1, a = 0;

clock_t start,end;
  start = clock();
for (i = 0; i < MAX; i++) {
if ((i % 2) == 0) {
is = 1;
} else {
is = -1;
}
a += is / (2 * i + 1); 
}
printf("%.15f\n", 4 *a);
end = clock();
  printf(" time [sec.] = %f \n", (double)(end - start) / CLOCKS_PER_SEC);
return 0;
}

$ gcc pi-leibniz.c -o pi-leibniz
$ ./pi-leibniz
3.141592652588050
 time [sec.] = 6.853473 
6.8秒かかりました。

OpenMPで並列化 pi-leibniz2.c
#include <stdio.h>
#include <time.h>
#include <omp.h>

#define MAX 1000000000

int main()
{
int i;
double is = 1, a = 0;

double start,end;
  start = omp_get_wtime();

#pragma omp parallel for private(is) reduction(+:a)
for (i = 0; i < MAX; i++) {
if ((i % 2) == 0) {
is = 1;
} else {
is = -1;
}
a += is / (2 * i + 1); 
}
printf("%.15f\n", 4 *a);
end = omp_get_wtime();
  printf(" time [sec.] = %lf \n", end - start);

return 0;
}

$ gcc pi-leibniz2.c -o pi-leibniz2 -fopenmp
$ ./pi-leibniz2
3.141592652589210
 time [sec.] = 1.899626 
 
並列化をすると1.8秒。 すごい効果が出ています。コア4つなので、早く終わります。

ところがpiの結果が異なります。

並列なし  3.141592652588050
並列あり  3.141592652589210

微妙に異なります。MAX 10にして調べました。

並列なし ループ毎に加算しています。
0 : 1.000000000000000
1 : 0.666666666666667
2 : 0.866666666666667
3 : 0.723809523809524
4 : 0.834920634920635
5 : 0.744011544011544
6 : 0.820934620934621
7 : 0.754267954267954
8 : 0.813091483679719
9 : 0.760459904732351
3.041839618929403

並列化あり 
6 : 0.076923076923077
7 : -0.073777203188968
8 : 0.135746606334842
9 : -0.126408782136336
3 : -0.007110536522301
4 : -0.015297671025225
5 : -0.106206761934316
0 : 0.893793238065684
1 : 0.560459904732351
2 : 0.760459904732351
3.041839618929402
スレッド毎に加算して最後にまとめ上げています。
当たり前といえば当たり前か。そういう仕組みだからしょうがないのか。

その辺りに誤差の原因がありそうですね。

誤差のレベルなので無視してもいいのかもしれませんが、わかりません。

並列なしと並列ありで突き合せて、確かめた方がいいですね。注意が必要です。

OpenMP 並列プログラミング1の続きです。

for文で並列化してみましょう。

forループの前に#pragma omp parallel forを追加するだけです。

#include <stdio.h>
#include <omp.h>

int main()
{

int i;

#pragma omp parallel for
for (i = 0; i < 5; i++) {
printf("i = %d: thread number = %d\n", i, omp_get_thread_num());
}

return 0;
}
$ ./openmp1
i = 4: thread number = 3
i = 3: thread number = 2
i = 2: thread number = 1
i = 0: thread number = 0
i = 1: thread number = 0 

$ ./openmp1
i = 0: thread number = 0
i = 1: thread number = 0
i = 3: thread number = 2
i = 4: thread number = 3
i = 2: thread number = 1

iの並びがランダムになっています。 


並列プログラミングのサンプルソースを読んでみました。OpenMPのポイントです。

一時的な変数にはprivate, 結果を残す変数にはreductionを変数に付ける。

配列を加算するプログラム。
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#define MAX 50000

int *box;

int main()
{
int i, total = 0;

box = malloc(MAX * sizeof(int));
for (i = 0; i <= MAX; i++) {
box[i] = i;
}

#pragma omp parallel
{
printf("thread = %d, %d\n", omp_get_thread_num(), omp_get_num_threads());

#pragma omp for reduction(+:total)
for (i = 0; i <= MAX; i++) {
total += box[i];
}
}
printf("total = %d\n", total);
return 0;
}

$ gcc openmp2.c -o openmp2 -fopenmp
$ ./openmp2
thread = 0, 4
thread = 2, 4
thread = 3, 4
thread = 1, 4
total = 1250025000
 
結果を残すのでtotal変数にはreduction(+:total)を付けています。
変数totalをスレッドで共有する必要があります。

実験でreductionなしにして、 #pragma omp for に変更してみてください。
間違った結果が帰ってきます。

参考サイト
OpenMP チュートリアル 

OpenMP* を使用したワークシェアリング 

linuxでOpenMP並列プログラミング1

corei5 のcpuで動かしています。まずは環境設定

virtual boxで仮想化しているので、その設定
34
設定 -> システム -> プロセッサーから最大の4コアを指定します。

それ以上を選択しても無効となります。下のPAE / NXはメモリー関係です。関係ないのでスルー。

cpuを確認
$ cat /proc/cpuinfo 
processor : 0
vendor_id : GenuineIntel
cpu family : 6
model : 61
model name : Intel(R) Core(TM) i5-5257U CPU @ 2.70GHz

processor : 1
vendor_id : GenuineIntel
cpu family : 6
model : 61
model name : Intel(R) Core(TM) i5-5257U CPU @ 2.70GHz

processor : 2
vendor_id : GenuineIntel
cpu family : 6
model : 61
model name : Intel(R) Core(TM) i5-5257U CPU @ 2.70GHz

processor : 3
vendor_id : GenuineIntel
cpu family : 6
model : 61
model name : Intel(R) Core(TM) i5-5257U CPU @ 2.70GHz

こんな感じで4つ認識していたらOK

 //gcc openmp0.c -o openmp0 -fopenmp

#include <stdio.h>
#include <omp.h>
int main()
{
   #pragma omp parallel
   {
     printf("Hello World! \n");
   }
}

$ ./openmp0
Hello World! 
Hello World! 
Hello World! 
Hello World! 

コア数分の4回Hello World!が出た! 
ここから並列プログラミングが始まります。
 

 

素数判定法5 繰り返し二乗法の続きになります。

前回はフェルマーテストの変数をlong long型で宣言していました。

long long型でも巨大な数字になると溢れてしまいます。

巨大な素数に対応するためにフェルマーテストをgmpライブラリーを使ってプログラミングしました。

えらい大変やった。脳みそ使うわ。

fermatgmp.c

#include<stdio.h>
#include<gmp.h>

int main()
{
mpz_t m, n, s, a;
mpz_t one, tmp1, tmp2, tmp3;

printf("数字を入力してください\n");
mpz_init(n);
gmp_scanf("%Zd", n);

mpz_init(m);
mpz_sub_ui(m, n, 1); //m = n - 1;

mpz_init(s);
mpz_set_str(s, "1", 10); // s = 1;

mpz_init(a);
mpz_set_str(a, "2", 10); // a = 2;

mpz_init(one);
mpz_set_str(one, "1", 10);

mpz_init(tmp1);
mpz_init(tmp2);
mpz_init(tmp3);

while(mpz_sgn(m)) { //while(m != 0) {
mpz_and(tmp1, m, one); //if (m & 1) {
if (mpz_sgn(tmp1)) {
mpz_mul(tmp2, s, a); //s = (s * a) % n;
mpz_mod(s, tmp2, n);

}
printf("m = ");
mpz_out_str(stdout, 10, m);
printf(" s = ");
mpz_out_str(stdout, 10, s);
printf(" a = ");
mpz_out_str(stdout, 10, a);
printf("\n");

//printf("m = %lld s = %lld a = %lld\n", m, s, a);

mpz_mul(tmp3, a, a); //a = (a * a) % n;
mpz_mod(a, tmp3, n);

mpz_div_2exp(m, m, 1); //m = m>>1;
}

mpz_out_str(stdout, 10, n);
if(!mpz_cmp_ui(s, 1)) {
printf("は素数です\n");
} else {
printf("は素数じゃないよ\n");
}

mpz_clear(m);
mpz_clear(n);
mpz_clear(s);
mpz_clear(a);
mpz_clear(tmp1);
mpz_clear(tmp2);
mpz_clear(tmp3);
return 0;

}

$ gcc fermatgmp.c -o fermatgmp -lgmp
$ ./fermatgmp
数字を入力してください
26588814640432479998443
26588814640432479998443は素数です

私のMac BookProだと一瞬で計算が終わります。
対応するC言語をコメントで残しています。
検証用に計算過程を表示させています。
邪魔であればコメントアウトしてくださぁい。

フェルマーテストの場合、素数を誤判定する可能性があります。

それで改良されたのがミラーラビン法です。

ミラーラビン法を使うことで誤判定の確率が減ります。

ミラーラビン法へ続く予定。いや続かないかも。多分続きます。
 

素数判定法5 フェルマーテストの続きです。

n が素数であるかを調べるために2 ^ nを計算しなければならない、ということでした。

n = 13のとき本来であれば

2 ^ 13 = 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2

2を12回掛けなければいけません。

これを簡単にするのが繰り返し二乗法です。

乗数を分解して、

2 ^ 13 = (2 ^ 1) * (2 ^ 4) * (2 ^ 8)  = 2 * 16 * 256 に分解するのが繰り返し二乗法です。

繰り返し二乗法を使えば計算回数を減らすことができます。

フェルマーの小定理を貼り付けておきます。

2 ^ (n - 1) mod n = 1 

c言語プログラム 
#include<stdio.h>

int main()
{
long long m, n;
long long s = 1;
long long a = 2;
printf("数字を入力してください\n");
scanf("%lld", &n);
m = n -1;
while(m != 0) {
if (m & 1) {
s = (s * a) %n;
}
printf("m = %lld s = %lld a = %lld\n", m, s, a);
a = (a * a) % n;
m = m>>1;
}
if (s == 1) {
printf("%lldは素数です\n", n);
} else {
printf("%lldは素数ではありません\n", n);
}
}

$ gcc fermat.c -o fermat

$ ./fermat

数字を入力してください

9999637

9999637は素数です 

素数判定法6 フェルマーテスト GMPライブラリ に続きます。

参考サイト
繰り返し二乗法

素数判定プログラム3 C言語GMPライブラリの続きです。

今までとは違った素数判定プログラムを書いていきます。

前回のgmpプログラムではmpz_probab_prime_p 関数を使いました。

実行すると計算が一瞬で終わります。一体どういう仕組みなのでしょうか。

mpz_probab_prime_p 関数は内部でMiller_Rabin 素数テストを行っています。

では、Miller Rabin素数テストとは何かを調べていきましょう。

Miller Rabin素数テストはフェルマーテストを改良して作られたものです。

というわけでまずはフェルマーテストの説明から。


nが素数の場合、

a ^ (n - 1) ≡ 1 (mod n)

が成立します。 a については 1 <= a < n という条件がついています。 ま、そりゃそうです。

aについては単純にa = 2 と置き換えて問題ありません。

上の式を変換すると

2 ^ (n - 1 ) mod n = 1

計算結果が1となると素数と判定します。

素数3の場合 2 ^ 2= 4,  4 % 3 = 1

素数5の場合 2 ^ 4 = 16, 16 % 5 = 1

素数7の場合 2 ^ 6 = 64, 64 % 7 = 1

素数の場合は必ず計算結果の余りが1 となります。これがフェルマーテストです。

1にならない場合は素数ではないので弾くことができます。

ところが素数でない場合も1 となることがあります。

341 = 11 * 31 は素数ではありませんが2となります。 

(2 ^ 340) % 341 = 1

前回までのプログラムは素数かどうかを割り算をして割り切れるかどうかを調べました。

フェルマーテストの場合は乗数計算となります。

しかも数が爆発します。

41が素数かどうかを調べるには 2 ^ 40です。計算すると1099511627776 になります。

2を40回掛け算しなければいけません。

これは大変、ということで簡単にするアルゴリズムがあります。

繰り返し二乗法です。

素数判定法5 繰り返し二乗法に続きます。

↑このページのトップヘ