Yuki Nakata's Blog

スーパー中ちゃん♪

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

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

GMPライブラリを使って素数判定プログラムを書いていきます。

gmpには素数を判定する関数が用意されています。

int mpz_probab_prime_p(mpz_t n, int reps);

戻り値で素数かどうかを判定します。

2が帰ってきたときは絶対素数 (definitely prime)
1のときは素数っぽい、(probably prime)
0のときは合成数つまり素数ではない、です。

repsはミラーラビン法のパラメータですが、通常25を指定しておけばいいそうです。

実行してみてください。一瞬で計算が終わりますよ。

$ gcc prime4.c -o prime4 -lgmp

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

int main(void)
{
mpz_t number;
int answer = 0;

  mpz_init(number);

mpz_set_str(number, "940461086986004843694934910131104208392131088686023657900173332902199657733778583489", 10);
 
answer =  mpz_probab_prime_p(number, 25); 

mpz_out_str(stdout, 10, number);
if (answer) {
printf(" は素数だよ\n");
} else {
printf(" は素数じゃないよ\n");
}
return 0;
}

素数判定プログラム C言語の続きです。

今度はGMPライブラリを使って素数判定プログラムを書いていきます。

GMPプログラミング初めてです。

ぶっつけ本番でも大丈夫でしょう。

変数をそれぞれ保存しないといけないのでプログラムが汚くなってしまいました。

C言語風のコメントも書いておきました。見づらくなったかも。w

GMPライブラリーもマスター出来ました。

//$ gcc prime3.c -o prime3 -lgmp

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

int isPrime(mpz_t number)
{
mpz_t i, four, twomod, threemod, sqrt, tmpmod, tmpmod2, iplus2;

mpz_init(i);
mpz_set_str(i, "5", 10); //10進数 5で初期化 i = 5;

mpz_init(sqrt);
mpz_sqrt(sqrt, number); 

mpz_init(four);
mpz_set_str(four, "4", 10);//10進数 4で初期化 four = 4;
mpz_init(twomod);
mpz_set_str(twomod, "0", 10);//10進数 0で初期化 twomod = 0;

mpz_init(threemod);
mpz_set_str(threemod, "0", 10);//10進数 0で初期化 threemod = 0;

mpz_init(tmpmod);
mpz_init(tmpmod2);
mpz_init(iplus2);

mpz_mod_ui(twomod, number, 2); //twomod = number % 2;
mpz_mod_ui(threemod, number, 3); //threemod = number % 3;

if (mpz_cmp_ui(number, 1) > 0 && 0 < mpz_cmp(four, number)) {
return 1;
} else if (mpz_get_ui(twomod) == 0 || mpz_get_ui(threemod) == 0) {
return 0;
} else {
for (i; 0 <= mpz_cmp(sqrt, i); mpz_add_ui(i, i, 6)) {
mpz_mod(tmpmod, number, i); //tmpmod = number % i;
mpz_add_ui(iplus2, i, 2); //iplus2 = i + 2;
mpz_mod(tmpmod2, number, iplus2); //tmpmod2 = number % iplus2;
if (mpz_get_ui(tmpmod) == 0 || mpz_get_ui(tmpmod2) == 0) {
return 0;
}
}
}

mpz_clear(i);
mpz_clear(sqrt);
mpz_clear(four);
mpz_clear(twomod);
mpz_clear(threemod);
mpz_clear(tmpmod);
mpz_clear(tmpmod2);
mpz_clear(iplus2);
return 1;
}

int main(void)
{
mpz_t number;
int answer = 0;

  mpz_init(number);

mpz_set_str(number, "150094772735952593", 10);

answer = isPrime(number);

mpz_out_str(stdout, 10, number);
if (answer) {
printf(" は素数だよ\n");
} else {
printf(" は素数じゃないよ\n");
}
mpz_clear(number);
return 0;
}

150094772735952593 は18桁の素数です。

ここに判定する数字を入力してください。
 
素数サンプル 楽勝レベル
31389448217
282437925089
1853028778786433
5559077746424707
150094772735952593

PCだと厳しい 一気に時間がかかります。終わらないかもね。
8862938260389989451257
26588814640432479998443

実はGMPライブラリには素数を判定する関数が用意されています。

それを呼べば一発で素数かどうかを判定することができます。

次回に続きます。 素数判定プログラム3 C言語GMPライブラリ

参考サイト
多倍長整数演算ライブラリ (GNU MP) 

素数一覧 1000万 

素数表 

素数を判定するプログラムを書いていきます。

単純にC言語で書いたのがコレ。

//$ gcc prime1.c -o prime1

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

int isPrime(unsigned long number)
{
int i;
for (i = 2; i < number; i++) {
if (number % i == 0) {
return 0;
}
}
return 1;
}

int main(void)
{
unsigned long number = 9999991;
int answer = 0;
answer = isPrime(number);
if (answer) {
printf("%lu は素数だよ\n", number);
} else {
printf("%lu は素数じゃないよ\n", number);
}
return 0;
}

numberのところに判定する数字を入力してください。

3箇所アルゴリズムを改良して早くしていきます。 

1. まず2や3で割り切れるなら素数ではありません。
 候補ではない数字を最初に除外します。これをエラトステネスの篩といいます。

2. 判定する数字を最後までチェックする必要はありません。
 判定する数字を2つの積で表します。
 たとえば100を2つの積で表してみましょう。
 2 * 50     2で素数かどうかチェック済み
    4 * 25 4で素数かどうかチェック済み
    5 * 20 5で素数かどうかチェック済み
   10 * 10 10で素数かどうかチェック済み

 つまりわざわざ100まで調べなくても良いのです。100のルートを取り10を上限とします。
 素数かどうかを調べるには10以下まで調べれば良いです。
  
3. 素数の数字に注目する。
 5以上100までの素数を書いてみます。
 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97
 素数というのは6の倍数の隣の数字になります。確認してください。
 6k ± 1で表せます。

 この3つを考慮してプログラムを書くと次になります。

//$ gcc prime2.c -o prime2

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

int isPrime(unsigned long number)
{
unsigned int i;
if (number <= 3 && number > 1) {
return 1;
} else if (number % 2 == 0 || number % 3 == 0) {
return 0;
} else {
for (i = 5; i * i <= number; i += 6) {
if (number % i == 0 || number % (i + 2) == 0) {
return 0;
}
}
}
return 1;
}

int main(void)
{
unsigned long number = 9999991;
int answer = 0;
answer = isPrime(number);
if (answer) {
printf("%lu は素数だよ\n", number);
} else {
printf("%lu は素数じゃないよ\n", number);
}
return 0;
}


これでだいぶ早くなったはずです。
for分がややこしいかもしれませんが、一度紙に書いてトレースして確認してみてください。
 
次回はGMPライブラリーを使って素数判定プログラムに改良していきます。

素数判定プログラム C言語 GMPライブラリー


 フィボナッチ数列を計算してみよう1 C言語 の続きです。

前回はC言語だけでフィボナッチ数列を計算していました。

それだけだとunsigned long型の上限があり、大きな数字が計算できません。

今回はgmpライブラリを使ってフィボナッチ数列 1000001番目を計算してみます。

fibo.cというファイル名で保存してください。 

//gcc fibo.c -o fibo -lgmp
#include <stdio.h>
#include <string.h>
#include <stdint.h>
#include <gmp.h>

void omg_i_love_leonardo_of_pisa(uint32_t num, mpz_t * result) 
mpz_t retval, last, tmp;
mpz_init(retval);
  mpz_init(last);
mpz_init(tmp);

uint32_t i = 1;
if(num == 0) {
return;
}
  mpz_set_ui(retval, 1U);
mpz_set_ui(last, 0U);

for(; i < num; i++) {
mpz_set(tmp, retval);
mpz_add(retval, retval, last);

mpz_set(last, tmp);
}

mpz_set(*result, retval);
 }

int main() 
{
uint32_t num;
mpz_t fibo;
mpz_init(fibo);

//1000001
omg_i_love_leonardo_of_pisa(1000001, &fibo);
mpz_out_str(stdout, 10, fibo);
printf("\n");
return 1;
}

$ gcc fibo.c -o fibo -lgmp
$ ./fibo

すぐに計算結果が出てきます。
科学技術すげー。

http://www.wolframalpha.com/input/?i=fibonacci+1000001


実はgmpにはmpz_fib_uiというフィボナッチ関数が用意されています。
mpz_fib_ui関数を実行するだけフィボナッチ数列を求めることができます。

fibo3.cというファイル名で保存してください。
#include <stdio.h>
#include <string.h>
#include <stdint.h>
#include <gmp.h>

int main()
{
  int n = 1000000;
  mpz_t fm2;
  mpz_inits(fm2, NULL);
  mpz_fib_ui(fm2, n);
  gmp_printf("fib(%d) = %Zd\n", n, fm2);
  return 1;
}
 
$ gcc fibo3.c -o fibo3 -lgmp
$ ./fibo3
 
参考サイト
1,000,000th Fibonacci Number One-Liner in C

フィボナッチ数列を計算してみよう C言語

フィボナッチ数列とは直前の2つの数字の和になっている数列のことです。

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89

数列を見たら分かりますね。C言語でプログラミングしていきましょう。

3通りの方法で書いています。

fibo1, fibo2, fibo3 と進むにつれて処理が早くなっています。

fibo1が最も単純で何度も再起呼び出しをしていて遅い方法です。

fibo2はデータを保存して再起呼び出しの回数を少なくしています。

fibo3が一番シンプルです。for文で回しているだけです。
初めからfibo3で書けよという感じですね。

関数の再起呼び出しは難しい。プログラムを読むだけで頭を使います。

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

#define MAX_VAL 47
unsigned long fibo2Arr[MAX_VAL];
unsigned long fibo3Arr[MAX_VAL];

unsigned long fibo1(int n)
{
unsigned long f;
switch (n) {
case 1:
case 2:
f = 1L;
break;
default:
f = fibo1(n - 1) + fibo1(n - 2);
break;
}
return (f);
}

unsigned long fibo2(int n)
{
unsigned long f;
if (fibo2Arr[n]) {
f = fibo2Arr[n];
return (f);
}
switch (n) {
case 1:
case 2:
f = 1L;
break;
default:
f = fibo2(n  - 1) + fibo2(n - 2);
break;
}
fibo2Arr[n] = f;
return (f);
}

unsigned long fibo3(int n)
{
int i;
fibo3Arr[1] = 1;
fibo3Arr[2] = 1;
for (i = 3; i <= n; i++) {
fibo3Arr[i] = fibo3Arr[i - 1] + fibo3Arr[i - 2];
}
return (fibo3Arr[n]);
}

int main(void)
{
int n = MAX_VAL;
clock_t start, end;
unsigned long result1, result2, result3;
start = clock();
result1 = fibo1(n);
end = clock();
printf("fibo1 %lu %f\n", result1, (double)end - start);
start = clock();
result2 = fibo2(n);
end = clock();
printf("fibo2 %lu %f\n", result2, (double)end - start);
start = clock();
result3 = fibo3(n);
end = clock();
printf("fibo3 %lu %f\n", result3, (double)end - start);
return 0;


フィボナッチ数列の47まで計算して、それぞれ掛かった時間を計測しています。

私の環境ではこんな感じ。fibo1だけがやたら遅いです。

fibo1 2971215073 15680165.000000

fibo2 2971215073 2.000000

fibo3 2971215073 1.000000

フィボナッチ数列の47は
2971215073です。

たった47番目で 約30億になりますよ、と。

Wolfram Alphaで検索して合っているか確認してみましょう。

http://www.wolframalpha.com/input/?i=fibonacci+47 

同様に100番目でも検索してみてください。 

http://www.wolframalpha.com/input/?i=fibonacci+100 

前回のlinuxで科学技術計算をしようの続きです。

円周率を1億桁計算しました!からのプログラムです。

pai.c


/* Calculate pi based on Chudnovsky algorithm (& Binary Splitting method), using GMP */
/* [1] Computation of 2700 billion decimal digits of Pi using a Desktop Computer,
       Fabrice Bellard, Feb 11 2010 (4th revision),
       http://bellard.org/pi/pi2700e9/pipcrecord.pdf */
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>

mpz_t A, B, C, C3over24;

void computePQT(mpz_t P, mpz_t Q, mpz_t T, int n1, int n2) {
  int m;
  if (n1 + 1 == n2) {
    mpz_t a;
    mpz_init(a);
    /* P(n2-1, n2) = - (2 * n2 - 1) * (6 * n2 - 5) * (6 * n2 - 1) */
    mpz_set_ui(P, (2 * n2 - 1));
    mpz_mul_ui(P, P, (6 * n2 - 5));
    mpz_mul_ui(P, P, (6 * n2 - 1));
    mpz_neg(P, P);
    /* Q(n2-1, n2) = C^3 * n2^3 / 24 */
    mpz_mul_ui(Q, C3over24, n2);
    mpz_mul_ui(Q, Q, n2);
    mpz_mul_ui(Q, Q, n2);
    /* a_n2 = (A + B * n2) */
    mpz_mul_ui(a, B, n2);
    mpz_add(a, a, A);
    /* T(n2-1, n2) = a_n2 * P(n2-1, n2) */
    mpz_mul(T, a, P);
    mpz_clear(a);
  } else {
    mpz_t P1, Q1, T1, P2, Q2, T2;
    mpz_init(P1);
    mpz_init(Q1);
    mpz_init(T1);
    mpz_init(P2);
    mpz_init(Q2);
    mpz_init(T2);
    m = (n1 + n2) / 2;
    computePQT(P1, Q1, T1, n1, m);
    computePQT(P2, Q2, T2, m, n2);
    /* P = P1 * P2 */
    mpz_mul(P, P1, P2);
    /* Q = Q1 * Q2 */
    mpz_mul(Q, Q1, Q2);
    /* T = T1 * Q2 + P1 * T2 */
    mpz_mul(T, T1, Q2);
    mpz_mul(P1, P1, T2);
    mpz_add(T, T, P1);
    mpz_clear(P1);
    mpz_clear(Q1);
    mpz_clear(T1);
    mpz_clear(P2);
    mpz_clear(Q2);
    mpz_clear(T2);
  }
}

int main (int argc, char* argv[]) {

  clock_t start, end;
  FILE* fp;

  int digits = 10000;
  int prec = digits * log2(10);
  int n = digits / 14;

  start = clock();

  mpf_t pi, temp;
  mpz_t P, Q, T;

  /* Initialize GMP numbers */
  mpf_set_default_prec(prec);
  mpf_init(pi);
  mpf_init(temp);
  mpz_init(P);
  mpz_init(Q);
  mpz_init(T);
  mpz_init(A);
  mpz_init(B);
  mpz_init(C);
  mpz_init(C3over24);

  /* Assignment */
  mpz_set_str(A, "13591409", 10);
  mpz_set_str(B, "545140134", 10);
  mpz_set_str(C, "640320", 10);
  mpz_mul(C3over24, C, C);
  mpz_mul(C3over24, C3over24, C);
  mpz_div_ui(C3over24, C3over24, 24);

  computePQT(P, Q, T, 0, n);

  /* pi = C ^ (1 / 2)     */
  mpf_set_z(temp, C);
  mpf_sqrt(pi, temp);
  /*      * C             */
  mpf_mul(pi, pi, temp);
  /*      * Q             */
  mpf_set_z(temp, Q);
  mpf_mul(pi, pi, temp);
  /*      / (T + A * Q)   */
  mpz_mul(Q, A, Q);
  mpz_add(Q, Q, T);
  mpf_set_z(temp, Q);
  mpf_div(pi, pi, temp);
  /*      / 12            */
  mpf_div_ui(pi, pi, 12);

  /* mpf_out_str(stdout, 10, digits, pi); */
  fp = fopen("output.txt", "w");
  if (fp == NULL) {
    printf("couldn't open output.txt");
    exit(EXIT_FAILURE);
  }
  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);
  mpf_out_str(fp, 10, digits, pi);

  mpf_clear(pi);
  mpf_clear(temp);
  mpz_clear(P);
  mpz_clear(Q);
  mpz_clear(T);
  mpz_clear(A);
  mpz_clear(B);
  mpz_clear(C);
  mpz_clear(C3over24);

  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

  return 0;
}

桁を変更するには
int digits = 10000;

この箇所で変更してください。
1億桁にするには
int digits = 100000000;

10万桁からやってみましょう。
$ #10万桁
$ gcc pai.c -o pai -lgmp
$ ./pai
0.001 s
0.002 s

$#1億桁
$ gcc pai.c -o pai -lgmp
$ ./pai
147.206 s
180.558 s

1億桁がたった3分で計算できてしまいました。
MacBook Pro 2015 メモリ8Gです。
凄いパソコンを全く有効利用できていませんね。
プログラムの中身は後で読みたいと思います。ちゃんと読みますよ。 

processinglifegame
以前、はむくんのライフゲームの世界を再現しようでライフゲームを再現しました。

今回はprocessingを使ってライフゲームを作成しましょう。

といっても他の人が作ってくれていますので、それを参考に改造していきます。

新米コンサルの日常 のライフゲームを改造しています。

このライフゲームでは2つのモードを設定しています。

セルを作成したり、削除するEditモードとライフゲームを動かすPlayモードです。

遊び方
起動すると最初にエディットモードになります。
マウスをクリックしてセルを作成してください。
スペースキーを押すとプレイモードとなります。
再びスペースキーを押すとプレイモードとエディットモードを切り替えることができます。

ゴスパーのグライダー銃も再現することができますよ。楽しめます。

final int SIZE_OF_CELL = 20;
final int FRAME_SIZE_OF_X = 800;
final int FRAME_SIZE_OF_Y = 800;
int num_of_row = FRAME_SIZE_OF_Y / SIZE_OF_CELL; //gyou  
int num_of_col = FRAME_SIZE_OF_X / SIZE_OF_CELL; //retsu
int num_of_cell = num_of_row * num_of_col;
Cell[][] cells = new Cell[num_of_row][num_of_col];
boolean game_state = false;


void setup(){
  size(FRAME_SIZE_OF_X, FRAME_SIZE_OF_Y);
  frameRate(20);
  for(int row=0; row<num_of_row; row++){
    for(int col=0; col<num_of_col; col++){
      int position_x = col * SIZE_OF_CELL;
      int position_y = row * SIZE_OF_CELL;
      cells[row][col] = new Cell(position_x, position_y, false);
    }
  }
  print("Edit Mode");
}

class Point{
  int row, col;
  Point(int my_row, int my_col){
    row = my_row;
    col = my_col;
  }
  int getRow(){
    return row;
  }
  int getCol(){
    return col;
  }
}

class Cell{
  int x, y;
  boolean state;

  Cell(int my_x, int my_y, boolean s){
    x = my_x;
    y = my_y;
    state = s;
  }
  
  void changeState(){
    state = !state;
  }
  
  void appear(){
    if(state){
      fill(34, 139, 34);
      //fill(50, 50, 50);
    }else{
      fill(50, 50, 50);
      //fill(255, 255, 255);
    }
    rect(x, y, SIZE_OF_CELL, SIZE_OF_CELL);
  }
}

void check_life_and_death(){
  ArrayList<Point> points = new ArrayList<Point>();
  for(int row=0; row<num_of_row; row++){
    for(int col=0; col<num_of_col; col++){
      //jibun no mawari check
      int num_of_alivers = 0;
      for(int r=-1; r<=1; r++){
        for(int c=-1; c<=1; c++){
          if(row+r<0 || row+r>=num_of_row || col+c<0 || col+c>=num_of_col){
            continue;
          } 
          if(cells[row+r][col+c].state){
            num_of_alivers++;
          }
        }
      }
      if(cells[row][col].state){ //jibun no bun wo hiku
        num_of_alivers--;
      }
      //joutai change no note
      if(cells[row][col].state){
        if(num_of_alivers<=1  || num_of_alivers>=4){
          points.add(new Point(row, col));
        }
      }else{
        if(num_of_alivers==3){
          points.add(new Point(row, col));
        }
      }
    }
  } 
  for(int i=0; i<points.size(); i++){
    Point point = points.get(i);
    int point_row = point.getRow();
    int point_col = point.getCol();
    cells[point_row][point_col].changeState();
  }
}

//random setup
void randomSetup(){
  //size(FRAME_SIZE_OF_X, FRAME_SIZE_OF_Y);
  frameRate(20);
  for(int row=0; row<num_of_row; row++){
    for(int col=0; col<num_of_col; col++){
      int position_x = col * SIZE_OF_CELL;
      int position_y = row * SIZE_OF_CELL;
      float coin_seed = random(1);
      boolean coin;
      if(coin_seed>0.8){
        coin = true;
      }else{
        coin = false;
      }
      cells[row][col] = new Cell(position_x, position_y, coin);
    }
  }  
}



void keyPressed(){
  
  
  if(key == ' '){
    game_state = !game_state;
  }
  
  if(game_state) {
    print("Play Mode  ");
  } else {
    print("Edit Mode  ");
  }
}

void mouseClicked(){
  int cell_col = mouseX/SIZE_OF_CELL;
  int cell_row = mouseY/SIZE_OF_CELL;
  cells[cell_row][cell_col].changeState();
  
  print("(" + cell_col + ",");
  print(cell_row + ") ");
}

void draw() {
  if(!game_state){
  }else{
    check_life_and_death();
  }
  for(int row=0; row<num_of_row; row++){
    for(int col=0; col<num_of_col; col++){
      cells[row][col].appear();
    }
  }
  //delay(500);
}

↑このページのトップヘ