コラッツの問題

コラッツの問題とは、ある自然数nに対し

• nが偶数ならnを2で割る
• nが奇数ならnを3倍して1を加える

を繰り返すと最終的にnは1になるというものです。詳しいことはWikipediaを読んで下さい＾＾；（このブログ、難しいことはいつもWikipediaまかせなのです）

これ、何が問題かというと「初期値にどんな自然数を与えても最終的には1に到達するのか」という問題です。（未だ未解決だが、3×2の53乗までは確かめられている）

コラッツの問題で行う操作をC言語プログラム化したのが次のコードです。（2007年のパソコン甲子園本戦問題にもなっています）

#include 

int collatz(int n)
{
  int count = 0;

  while(n != 1){
    n = (n%2 == 0) ? (n / 2) : (n*3 + 1);
    printf("%d\n", n);
    count++;
  }

  return count;
}

int main(void)
{
  int n, count;

  printf("Input n:");
  scanf("%d", &n);
  count = collatz(n);
  printf("COUNT:%d\n", count);

  return 0;
}

例えば、nに5を入力すると5,16,8,4,2,1で繰り返した回数が5回となります。

フィボナッチ数を求める

フィボナッチ数、それは自然界の現象にも数多くみられる不思議な数です。（詳しくはWikipediaを）

n番目のフィボナッチ数をFnで表すと、

と定義されます。これをプログラムにすると

  i, n, f_old, f, f_new;
  f_old = 0; f_new = 1;
  Input n;
 
  for(i=2 ; i<=n ; i++){
    f_new = f + f_old;
    print f_new;
    f_old = f;
    f = f_new;
  }

注）これは、なんちゃってコードです。

このコードを見てもわかるとおり、フィボナッチ数を求めるプログラムは繰り返しと変数の更新だけで実現できます。

C言語を始めたばっかりの頃にフィボナッチ数を求めるプログラムを課題にされたときは難しく感じたのですが・・・。（しみじみです）

べき剰余を求める

RSA暗号プログラムでべき剰余を計算する必要がありましたが、べき剰余を計算するプログラムをまだ紹介していなかったのでここで。

#include <stdio.h>
#define N 28

int imod(int a, int n)
{
  int r;
  r = (a + n) % n;
  if(r < 0) r += n;
  return r;
}

int pow_mod(int a, int e, int n)
{
  int r = 1;
  a = imod(a, n);
  while(e > 0){
    if ((e & 1) == 1) r = imod(r*a, n); 
    e >>= 1;
    a = imod(a*a, n);
  }
  return r;
}

int main(void)
{
  int a, e;
  printf("input a:"); scanf("%d", &a);
  printf("input e:"); scanf("%d", &e);

  printf("%d^%d (mod %d) = %d\n"
          , a, e, N, pow_mod(a, e, N));

  return 0;
}

べき剰余を計算するのは関数pow_modです。

RSA暗号プログラム

本丸のRSA暗号プログラムです。RSA暗号/復号を処理する関数は次のようになります。

//---RSA Encryption & Decryption---
void rsa_angou(int xy[2], int key[2], int stu[3])
{
  int m, code;
  (mode == 1) ? (m = bind2(xy)) : (m = bind3(stu));
  code = pow_mod(m, key[0], key[1]);
  (mode == 1) ? (split3(code, stu)) : (split2(code, xy));
}

引数には（2文字を格納した配列、鍵を格納した配列、3文字を格納した配列）を与えます。関数内で使われているmodeという変数はmode == 1のときは暗号モード、mode != 1のときは復号モードで処理するために設定したものです。

また、この関数内にはpow_mod(a, b, c)やbind、splitといった関数が存在していますが、pow_mode(a, b, c)はaのb乗(mod n)を計算するもの（後述）、bindやsplitは次のようになっています。( #define N 28 )

//---{x, y} --> x*N + y---
int bind2(int xy[2])
{
  return (xy[0]*N + xy[1]);
}
//---{s, t, u} --> s*N^2 + t*N + u---
int bind3(int stu[3])
{
  return (stu[0]*N*N + stu[1]*N + stu[2]);
}

//---m --> {x, y}---
void split2(int m, int xy[2])
{
  if(m > N*N){
    fprintf(stderr, "Error Argument 'm'.\n");
    exit(1);
  }
  xy[0] = m / N;
  xy[1] = m - xy[0]*N;
}
//---m --> {s, t, u}---
void split3(int m, int stu[3])
{
  if(m > N*N*N){
    fprintf(stderr, "Error Argument 'm'.\n");
    exit(1);
  }
  stu[0] = m / (N*N); 
  stu[1] = (m - stu[0]*N*N) / N;
  stu[2] = m - stu[0]*N*N - stu[1]*N;
}

何故このような関数が必要かというと、今回私がつくったRSA暗号プログラムはブロック長2で暗号化しています。つまり、2文字ずつ読み込んで暗号化しているのですが、RSA暗号で2文字を暗号化すると3文字になって出力されてきます。（←このページの下の方にある例を読むとわかると思います）したがって、

暗号化：2文字をbind2で計算してmする→mを暗号化「m^key1 mod key2」（このとき、(key1, key2) = 公開鍵(e, n)である）→split3でmを3文字に変換して出力

復号化：3文字をbind3で計算してmする→mを復号化「m^key1 mod key2」（このとき、(key1, key2) = 秘密鍵(d, n)である）→split2でmを2文字に変換して出力

このような処理をするためにbindやsplitがあるわけです。

次にpow_modによるaのb乗(mod n)の計算についてですが、普通はa^bを計算した後にnで割った余りを求めればいいわけですが、これは’b’の値が小さいときのみに有効な方法です。RSA暗号では、’b’に該当する’e’や’d’に大きな値（今回は3000までの素数しか使っていませんが、それでも6000くらい）が代入されるので時間がかかります。（単純に言えばb=6000の場合、5999回の乗算を必要とする）

これを高速に計算するにはバイナリ法（このWikipediaのページにはコードもあるので参考になります）と呼ばれる手法を用います。

RSA暗号の公開鍵と秘密鍵を生成する

いよいよRSA暗号です。（最近忙しくて遅れました）今回はRSA暗号で用いる公開鍵と秘密鍵を生成するプログラムを紹介します。

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

#define RANDOMIZE() srand(time(NULL))
#define RANDOM(x)   (rand()%(x))
#define N   28
#define MAX 3000
#define PRM 500

//---Prime Number Search---
int prime_search(int ptable[])
{
  int i, j, count=0;
  int p[MAX];
  for(i=0 ; i<MAX ; i++){
    p[i] = 0;
  }
  p[0] = 1;
  for(i=2 ; i<=MAX/2 ; i++){
    for(j=2 ; i*j<=MAX ; j++){
      if(p[i*j-1] == 0){
       p[i*j-1] = 1;
      }
    }
  }
  for(i=0 ; i<MAX ; i++){
    if(p[i] == 0){
      ptable[count] = i+1;
      count++;
    }
  }
  return count;
}

int gcd(int x, int y)
{
  while(y != 0) {
    int temp = y;
    y = x%y;
    x = temp;
  }
  return x;
}

int add_gcd(int a, int b)
{
  int r_pre = a, r = b, r_next;
  int x_pre = 1, x = 0, x_next;
  int y_pre = 0, y = 1, y_next;
  int ans_x, ans_y, ans_gcd, ans, q, k = 1;
  while((r_next = r_pre%r) != 0){
   q = r_pre / r;
   k++;
   x_next = q*x + x_pre;
   y_next = q*y + y_pre;
   r_pre = r;
   r = r_next;
   x_pre = x;
   x = x_next;
   y_pre = y;
   y = y_next;
  }
  ans_x = pow(-1, k) * x;
  ans_y = pow(-1, k+1) * y;
  ans_gcd = (ans_x * a)+(ans_y * b);
  if(ans_gcd == 1) ans = (ans_x + b)%b;
  else ans = 0;
  return ans;
}

int main()
{
  int e, d, p, q, n, euler;
  int count, ptable[PRM];
  
  RANDOMIZE();
  count = prime_search(ptable);
  
  do{
    p = ptable[RANDOM(count)];
    q = ptable[RANDOM(count)];
    n = p * q;
  }while(n > N*N*N || n < N*N);
  euler = (p - 1) * (q - 1);
  do{
    e = RANDOM(euler);
  }while(gcd(e, euler) != 1);
  d = add_gcd(e, euler);

  printf("(e, n) = (%d, %d)\n(d, n) = (%d, %d)\n",
         e, n, d, n);
  
  return 0;
}

今回は簡単のために3000までの素数を扱うことにします。鍵生成までのステップは次のようになります。

1．素数表からp、qを選ぶ（疑似乱数を使用）
2．n = p * q を計算し、N^2 < n < N^3 （今回はN=28）であることを確かめる（満たさない場合はp、qを選びなおす）
3．オイラー関数φ(n) = (p-1)(q-1)を計算する
4．オイラー関数φ(n)より小さい正の整数eを選ぶ（eはオイラー関数と互いに素である）
5．φ(n)を法としたeの逆元を求め（拡張ユークリッドを使用）、この逆元をdとする
6．(e, n)を公開鍵、(d, n)を秘密鍵とする

エラトステネスの篩

エラトステネスの篩は素数判定法の代表的なアルゴリズムです。簡単に言えば、「１より大きい数の表から合成数を取り除き、残ったものを素数と判定する」ということです。

例えば、3000までの素数判定は次のようになります。

#include <stdio.h>
#define MAX 3000

int main(void)
{
  int i, j;
  int p[MAX];

  //---初期化---
  for(i=0 ; i<MAX ; i++) p[i] = 0;
  p[0] = 1;
  //---倍数を除外---
  for(i=2 ; i<=MAX/2 ; i++){
    for(j=2 ; i*j<=MAX ; j++){
      if(p[i*j-1] == 0)
        p[i*j-1] = 1;
    }
  }
  for(i=0 ; i<MAX ; i++){
    if(p[i] == 0)
      printf("%d\n", i+1);
  }

  return 0;
}

ヒル暗号プログラム解説

昨日は眠くてできなかったヒル暗号プログラムの解説をしたいと思います。前回に引き続きプログラムを表示します。

void hill(int u[2], int key[4], int hkey[4])
{
  int i, det_a, t, angou[2];
  det_a = (key[0]*key[3] - key[1]*key[2] + 28) % 28;
  if(det_a < 0) det_a += 28;
  t = add_gcd(det_a, 28);
  if(t == 0){
   fprintf(stderr, "This keys are incorrect!\n");
   exit(1);
  }
  angou[0] = ((key[0]*u[0] + key[1]*u[1])+28) % 28;
  angou[1] = ((key[2]*u[0] + key[3]*u[1])+28) % 28;
  for(i=0 ; i<2 ; i++){
    if(angou[i] < 0) angou[i] += 28;
  }
  //---Hukugou Kagi---
  hkey[0] = (t*key[3]+28) % 28;
  hkey[1] = (t*key[1]*(-1)+28) % 28;
  hkey[2] = (t*key[2]*(-1)+28) % 28;
  hkey[3] = (t*key[0]+28) % 28;
  for(i=0 ; i<4 ; i++){
    if(hkey[i] < 0) hkey[i] += 28;
  }
  u[0] = angou[0];
  u[1] = angou[1];
}

このヒル暗号プログラムはブロック長２で暗号化を行います。（平文から２文字ずつ読み込むということ）

このvoid関数hillの引数は（数値化した２文字、４つの暗号鍵、４つの復号鍵）となっていて、この関数に２文字と暗号鍵を渡し、暗号化された２文字と復号鍵を取り出すようにしました。（復号鍵を取り出すのはプログラム実行時に復号鍵を表示するためですが、この関数内で表示するようにすれば取り出す必要はないです）

この関数の処理は基本的にはヒル暗号のロジックどおりに行っています。

1．暗号鍵A=（a&b c&d)の行列式を計算し、GCD(detA、28）=1であることを確認（このときプログラム内では拡張ユークリッドを用いて、GCD(detA、28)=1であればdetAの逆元を、そうでなければ0を変数tに代入している）

2．暗号化する

3．復号鍵を求める

4．暗号化した文字を取り出すために変数を入れ替え

このプログラムで頻繁に行っている

(var + 28) % 28;
if(var < 0) var += 28;

という操作はmod 28を守るために行っています。つまり、0～27までの数字のみ使用できる。29なら1、-3なら25と同じになります。そのため、1行目では例のようにvarの値を変換します。では何故2行目が必要かというと、varに入っている値が-28より小さい場合には1行目で変換しても、-1～-27になってしまうので、2行目のようにしてmod 28に収めてやる必要があります。（varの値の絶対値をとってから1行目の計算をすれば2行目は必要ないかな？）

ヒル暗号完成

ヒル暗号プログラムが完成しました。

void hill(int u[2], int key[4], int hkey[4]) 
{ 
  int i, det_a, t, angou[2]; 
  det_a = (key[0]*key[3] - key[1]*key[2] + 28) % 28; 
  if(det_a < 0) det_a += 28; 
  t = add_gcd(det_a, 28); 
  if(t == 0){ 
    fprintf(stderr, "This keys are incorrect!\n"); 
    exit(1); 
  } 
  angou[0] = ((key[0]*u[0] + key[1]*u[1])+28) % 28; 
  angou[1] = ((key[2]*u[0] + key[3]*u[1])+28) % 28; 
  for(i=0 ; i<2 ; i++){ 
    if(angou[i] < 0) angou[i] += 28; 
  } 
  //---Hukugou Kagi--- 
  hkey[0] = (t*key[3]+28) % 28; 
  hkey[1] = (t*key[1]*(-1)+28) % 28; 
  hkey[2] = (t*key[2]*(-1)+28) % 28; 
  hkey[3] = (t*key[0]+28) % 28; 
  for(i=0 ; i<4 ; i++){ 
    if(hkey[i] < 0) hkey[i] += 28; 
  } 
  u[0] = angou[0]; 
  u[1] = angou[1]; 
}

解説はまた今度ということで・・・。（とっても眠い）

ところで、ヒル暗号の暗号鍵を考えるのは面倒なので鍵生成プログラムを作りました。

#include <stdio.h> 
#include <math.h> 

int add_gcd(int a, int b) 
{ 
  int r_pre = a, r = b, r_next; 
  int x_pre = 1, x = 0, x_next; 
  int y_pre = 0, y = 1, y_next; 
  int ans_x, ans_y, ans_gcd, ans, q, k = 1; 
  while((r_next = r_pre%r) != 0){ 
    q = r_pre / r; 
    k++; 
    x_next = q*x + x_pre; 
    y_next = q*y + y_pre; 
    r_pre = r; 
    r = r_next; 
    x_pre = x; 
    x = x_next; 
    y_pre = y; 
    y = y_next; 
  } 
  ans_x = pow(-1, k) * x; 
  ans_y = pow(-1, k+1) * y; 
  ans_gcd = (ans_x * a)+(ans_y * b); 
  if(ans_gcd == 1) ans = (ans_x + b)%b; 
  else ans = 0; 
  return ans; 
} 

int main(void) 
{ 
  int i, j, k, l; 
  int det, check; 

  for(l=1 ; l<=28 ; l++){ 
   for(k=1 ; k<=28 ; k++){ 
    for(j=1 ; j<=28 ; j++){ 
     for(i=1 ; i<=28 ; i++){ 
      det = (i*l - j*k + 28) % 28; 
      check = add_gcd(det, 28); 
      if(check != 0) 
       printf("key=(%d %d %d %d)\n", i, j, k, l); 
     } 
    } 
   } 
  } 

  return 0; 
} 

3乗根を求める

2乗根（平方根）を求めるならsqrt()（C言語のmathライブラリの関数）で求められますが、3乗根となると手が止まる人もいると思います。（私もその一人）

3乗根が整数になる値（8,27,64,125,216,...）ならばpow(n, 1.0/3.0)で求めてもなんら問題ないと思いますが、3乗根が実数になる値となると精度が気になります。

そんなときに私が思いつくのはNewton法くらいです。（Newton法の詳細はリンク先のWikipediaで）つまり、f(x) = x^3 - n = 0 とおいてNewton法により x の値を求めるわけです。

う～ん、どちらを使えばいいんだろう？

pow()のほうが汎用性高いけど、Newton法のほうが精度よさそうだし・・・

アドバイス欲しいなぁ

復号を考慮した暗号プログラム

このブログでは以前、シーザー（シフト）暗号とアフィン暗号プログラムを紹介しましたが、このとき作ったプログラムは復号のことを考慮していませんでした。（つまり、暗号化したはいいものの、それを復号できない代物）

復号化を考慮したものは次のようになります。

1．シフト暗号（使用する文字数28の場合）

int angouka(int num, int key)
{
  int angou;

  key = (key + 28) % 28;
  if(key < 0){
    key = key + 28;
  }
  angou = (num + key) % 28;

  return angou;
}

[解説]
この関数の引数のnumは入力した文字に割り振られた番号（別関数によって規格化されたもの、例えばa=0, b=1,.....,z=25）で、keyは（暗号/復号）鍵です。復号するには復号鍵の扱い方を考えなければなりません。（復号鍵=暗号鍵×(-1)としています）
つまり、鍵が－3なら28-3=25というようにします。このとき復号鍵の絶対値＞使用する文字数だと復号できないので上のコードのようにすることで鍵の絶対値が28より大きくならないようにする必要があります。

2．アフィン暗号（使用文字数28）

int affine(int num, int key1, int key2)
{
  int angou;
  int check = gcd(key1, 28);
  
  if(check == 0){
    fprintf(stderr, "This Key is incorrect!\n");
    exit(1);
  }else{
    if(mode == 1){
      angou = (num*key1 + key2) % 28;
    }else if(mode == -1){
      angou = ( key1*(num - key2) ) % 28;
      if(angou < 0) angou = angou + 28;
    }
  }
  
  return angou;
}

[解説]
このプログラムでは暗号モード(mode=1)と復号モード(mode=-1)を選んで動作させるようにしてあります。（復号鍵1は暗号鍵1の逆元（mod 28)、復号鍵2=暗号鍵2）ここでも復号するときに注意が必要で、（鍵1+鍵2）＞使用文字数になると正しく復号できません。そこで、このコードのように書いたら正しく復号できるようになりました。（動作確認はしましたが、思いつきなので、深く考察していません）