Since June 16, 2000
STUDIO KAMADAJapanese to English by Excite
フラクタル | サイトマップ | ホーム
フラクタルプログラム

mandelbmp.c…マンデルブロ集合を描画する
julia3bmp.c…3次関数のジュリア集合を描画する
julia3cbmp.c…3次関数のジュリア集合を関数の定数項について描画する
j3cxbmp.c…3次関数のジュリア集合を関数の定数項について描画する(拡張付き)
complex.h…複素数演算マクロ定義ファイル
output.c…彩色してBMPファイルを生成する


mandelbmp.c

説明
 mandelbmp.cは、マンデルブロ集合のBMPデータを生成するプログラムです。complex.hが必要です。output.cと一緒にgccでコンパイルします。

mandelbmp.c
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#include <stdio.h>
#include <setjmp.h>
#include <time.h>
#include "complex.h"

#define USAGE "\
マンデルブロ集合の画像をBMPファイルに出力します\n\
>mandelbmp BMPファイル名 実部最小 虚部最小 辺の長さ 縦横のピクセル数 回数制限\n\
省略値        なし       -2.0     -2.0      4.0          256          256\n\
"

/* 描画範囲のデフォルト */
#define RE_MIN (-2.0)
#define IM_MIN (-2.0)
#define SC_LEN (4.0)

/* 縦横のピクセル数のデフォルト */
#define SC_CNT 256

/* 回数制限のデフォルト */
#define MX_CNT 256

char *bmpName;  /* BMPファイル名の先頭 */

double reMin = RE_MIN;  /* 実部最小 */
double imMin= IM_MIN;  /* 虚部最小 */
double scLen = SC_LEN;  /* 辺の長さ */
int scCnt = SC_CNT;  /* 縦横のピクセル数(1〜32767) */
int mxCnt = MX_CNT;  /* 回数制限(1〜1000000) */

int *cntTable;  /* 結果を記録する領域へのポインタ */
int *paTable;  /* パレットテーブルへのポインタ */

/* パレットテーブルを作る */
int *make_palet_table(int mxCnt);

/* 24ビットのBMPファイルを出力する */
int output_bmp(char *bmpName, int *paTable, int *cntTable, int width, int height);

/* 初期値cでマンデルブロ集合の漸化式の反復回数を返す */
int iteration(complex c, int limit_cnt)
{
  complex z;
  int cnt;

  if (setjmp(cjmpenv)) {
    return limit_cnt;  /* 計算できなかった */
  }

  /* 漸化式の反復 */
  cnt = limit_cnt;
  z = c;
  while (cnt--) {
    if (cabs2(z) > 4.0) {
      break;
    }
    z = cadd(cpow2(z), c);  /* マンデルブロ集合の漸化式 z=z^2+c */
  }
  return limit_cnt - (cnt + 1);
}

/* 描画ルーチンの本体 */
void draw_body()
{
  complex c0 = tocomplex(reMin, imMin), c;
  double step = scLen / (double)scCnt;
  int *p = cntTable;
  int reCnt, imCnt;

  /* 虚軸方向のループ */
  imCnt = scCnt;
  while (imCnt--) {
    c = c0;

    /* 実軸方向のループ */
    reCnt = scCnt;
    while (reCnt--) {
      *p++ = iteration(c, mxCnt);  /* 漸化式の反復回数を数えて記録する */

      Re(c) += step;  /* 実軸方向に進む */
    }

    Im(c0) += step;  /* 虚軸方向に進む */
  }
}

int main(int argc, char *argv[])
{
  clock_t tm0, tm1;

  /* コマンドラインのパラメータを取得 */
  switch (argc) {
  case 7:
    sscanf(argv[6], "%d", &mxCnt);
  case 6:
    sscanf(argv[5], "%d", &scCnt);
  case 5:
    sscanf(argv[4], "%lf", &scLen);
    sscanf(argv[3], "%lf", &imMin);
    sscanf(argv[2], "%lf", &reMin);
  case 2:
    bmpName = argv[1];
    if (bmpName[0] != '-') {
      break;
    }
  default:
    fprintf(stderr, USAGE);
    return 1;
  }

  /* 結果を記録する領域を用意する */
  if ((cntTable = (int *)malloc(sizeof(int) * scCnt * scCnt)) == NULL) {
    fprintf(stderr, "メモリが不足しています\n");
    return 1;
  }

  /* 計測開始 */
  {
    clock_t t = clock();
    while ((tm0 = clock()) == t);
  }

  /* 描画ルーチンの本体 */
  draw_body();

  /* 計測終了 */
  tm1 = clock();

  /* 所要時間を表示 */
  printf("所要時間は %.8g 秒です.\n", ((double)(tm1 - tm0)) / CLK_TCK);

  /* パレットテーブルを作る */
  if ((paTable = make_palet_table(mxCnt)) == NULL) {
    return 1;
  }

  /* BMPファイルを出力する */
  if (output_bmp(bmpName, paTable, cntTable, scCnt, scCnt)) {
    return 1;
  }

  return 0;
}


julia3bmp.c

説明
 julia3bmp.cは、3次関数のジュリア集合のBMPデータを生成するプログラムです。complex.hが必要です。output.cと一緒にgccでコンパイルします。

julia3bmp.c
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
#include <stdio.h>
#include <setjmp.h>
#include <time.h>
#include "complex.h"

#define USAGE "\
3次のジュリア集合の画像をBMPファイルに出力します\n\
>julia3bmp BMPファイル名 実部最小 虚部最小 辺の長さ 縦横のピクセル数 回数制限 収束半径 2次係数実部 2次係数虚部 1次係数実部 1次係数虚部 定数項実部 定数項虚部\n\
省略値        なし       -2.0     -2.0      4.0          256          256      0.01       0.0          0.0         0.0         0.0       -1.0        0.0\n\
収束半径=収束したと判断する|f(z)|の上限\n\
f(z)=z^3+A*z^2+B*z+C  f'(z)=3*z^2+2*A*z+B\n\
"

/* 描画範囲のデフォルト */
#define RE_MIN (-2.0)
#define IM_MIN (-2.0)
#define SC_LEN (4.0)

/* 縦横のピクセル数のデフォルト */
#define SC_CNT 256

/* 回数制限のデフォルト */
#define MX_CNT 256

/* 収束半径のデフォルト */
#define ROT_F_LIMIT (0.01)

/* 係数と定数項のデフォルト */
#define Re_A (0.0)
#define Im_A (0.0)
#define Re_B (0.0)
#define Im_B (0.0)
#define Re_C (-1.0)
#define Im_C (0.0)

char *bmpName;  /* BMPファイル名の先頭 */

double reMin = RE_MIN;  /* 実部最小 */
double imMin= IM_MIN;  /* 虚部最小 */
double scLen = SC_LEN;  /* 辺の長さ */
int scCnt = SC_CNT;  /* 縦横のピクセル数(1〜32767) */
int mxCnt = MX_CNT;  /* 回数制限(1〜1000000) */

double rot_f_limit = ROT_F_LIMIT;  /* 収束半径 */
double rot_f_limit2;  /* 収束半径の2乗 */

complex A, B, C;  /* 係数と定数項 */

int *cntTable;  /* 結果を記録する領域へのポインタ */
int *paTable;  /* パレットテーブルへのポインタ */

/* パレットテーブルを作る */
int *make_palet_table(int mxCnt);

/* 24ビットのBMPファイルを出力する */
int output_bmp(char *bmpName, int *paTable, int *cntTable, int width, int height);

/* 3次のジュリア集合の素となる3次関数の定義 */
/* f(z)=z^3+A*z^2+B*z+C, df_f(z)=3*z^2+2*A*z+B */
/* マクロ版と関数版のどちらか速い方を選択する */
#if 0
#define f(z) cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), A)), cmul((z), B)), C)
#define df_f(z) cadd(cadd(cremul(cpow2(z), 3.0), cremul(cmul((z), A), 2.0)), B)
#else
complex f(complex z)
{
  return cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), A)), cmul(z, B)), C);
}
complex df_f(complex z)
{
  return cadd(cadd(cremul(cpow2(z), 3.0), cremul(cmul(z, A), 2.0)), B);
}
#endif

/* 初期値cでf(z)に対するNewton法の漸化式の反復回数を返す */
int iteration(complex c, int limit_cnt)
{
  complex z;
  int cnt;

  if (setjmp(cjmpenv)) {
    return limit_cnt;  /* 計算できなかった */
  }

  /* 漸化式の反復 */
  cnt = limit_cnt;
  z = c;
  while (cnt--) {
    if (cabs2(f(z)) <= rot_f_limit2) {
      break;
    }
    z = csub(z, cdiv(f(z), df_f(z)));  /* f(z)に対するNewton法の漸化式 rec_f(z)=z-f(z)/df_f(z) */
  }
  return limit_cnt - (cnt + 1);
}

/* 描画ルーチンの本体 */
void draw_body()
{
  complex c0 = tocomplex(reMin, imMin), c;
  double step = scLen / (double)scCnt;
  int *p = cntTable;
  int reCnt, imCnt;

  /* 虚軸方向のループ */
  imCnt = scCnt;
  while (imCnt--) {
    c = c0;

    /* 実軸方向のループ */
    reCnt = scCnt;
    while (reCnt--) {
      *p++ = iteration(c, mxCnt);  /* 漸化式の反復回数を数えて記録する */

      Re(c) += step;  /* 実軸方向に進む */
    }

    Im(c0) += step;  /* 虚軸方向に進む */
  }
}

int main(int argc, char *argv[])
{
  clock_t tm0, tm1;

  /* 係数と定数項のデフォルトを設定 */
  A = tocomplex(Re_A, Im_A);
  B = tocomplex(Re_B, Im_B);
  C = tocomplex(Re_C, Im_C);

  /* コマンドラインのパラメータを取得 */
  switch (argc) {
  case 14:
    sscanf(argv[13], "%lf", &Im(C));
    sscanf(argv[12], "%lf", &Re(C));
  case 12:
    sscanf(argv[11], "%lf", &Im(B));
    sscanf(argv[10], "%lf", &Re(B));
  case 10:
    sscanf(argv[9], "%lf", &Im(A));
    sscanf(argv[8], "%lf", &Re(A));
  case 8:
    sscanf(argv[7], "%lf", &rot_f_limit);
  case 7:
    sscanf(argv[6], "%d", &mxCnt);
  case 6:
    sscanf(argv[5], "%d", &scCnt);
  case 5:
    sscanf(argv[4], "%lf", &scLen);
    sscanf(argv[3], "%lf", &imMin);
    sscanf(argv[2], "%lf", &reMin);
  case 2:
    bmpName = argv[1];
    if (bmpName[0] != '-') {
      break;
    }
  default:
    fprintf(stderr, USAGE);
    return 1;
  }

  /* 収束半径を2乗しておく */
  rot_f_limit2 = rot_f_limit * rot_f_limit;

  /* 結果を記録する領域を用意する */
  if ((cntTable = (int *)malloc(sizeof(int) * scCnt * scCnt)) == NULL) {
    fprintf(stderr, "メモリが不足しています\n");
    return 1;
  }

  /* 計測開始 */
  {
    clock_t t = clock();
    while ((tm0 = clock()) == t);
  }

  /* 描画ルーチンの本体 */
  draw_body();

  /* 計測終了 */
  tm1 = clock();

  /* 所要時間を表示 */
  printf("所要時間は %.8g 秒です.\n", ((double)(tm1 - tm0)) / CLK_TCK);

  /* パレットテーブルを作る */
  if ((paTable = make_palet_table(mxCnt)) == NULL) {
    return 1;
  }

  /* BMPファイルを出力する */
  if (output_bmp(bmpName, paTable, cntTable, scCnt, scCnt)) {
    return 1;
  }

  return 0;
}


julia3cbmp.c

説明
 julia3cbmp.cは、3次関数のジュリア集合を関数の定数項について描画したBMPデータを生成するプログラムです。complex.hが必要です。output.cと一緒にgccでコンパイルします。

julia3cbmp.c
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
#include <stdio.h>
#include <setjmp.h>
#include <time.h>
#include "complex.h"

#define USAGE "\
3次のジュリア集合の定数項の画像をBMPファイルに出力します\n\
>julia3cbmp BMPファイル名 実部最小 虚部最小 辺の長さ 縦横のピクセル数 回数制限 収束半径 2次係数実部 2次係数虚部 1次係数実部 1次係数虚部 初期値実部 初期値虚部\n\
省略値          なし       -2.0     -2.0      4.0          256          256      0.01       0.0         0.0        -1.0         0.0        0.0        0.0\n\
収束半径=収束したと判断する|f(z)|の上限\n\
f(z)=z^3+A*z^2+B*z+C  f'(z)=3*z^2+2*A*z+B\n\
"

/* 描画範囲のデフォルト */
#define RE_MIN (-2.0)
#define IM_MIN (-2.0)
#define SC_LEN (4.0)

/* 縦横のピクセル数のデフォルト */
#define SC_CNT 256

/* 回数制限のデフォルト */
#define MX_CNT 256

/* 収束半径のデフォルト */
#define ROT_F_LIMIT (0.01)

/* 係数と初期値のデフォルト */
#define Re_A (0.0)
#define Im_A (0.0)
#define Re_B (-1.0)
#define Im_B (0.0)
#define Re_Z (0.0)
#define Im_Z (0.0)

char *bmpName;  /* BMPファイル名の先頭 */

double reMin = RE_MIN;  /* 実部最小 */
double imMin= IM_MIN;  /* 虚部最小 */
double scLen = SC_LEN;  /* 辺の長さ */
int scCnt = SC_CNT;  /* 縦横のピクセル数(1〜32767) */
int mxCnt = MX_CNT;  /* 回数制限(1〜1000000) */

double rot_f_limit = ROT_F_LIMIT;  /* 収束半径 */
double rot_f_limit2;  /* 収束半径の2乗 */

complex A, B, C, Z;  /* 係数と定数項と初期値 */

int *cntTable;  /* 結果を記録する領域へのポインタ */
int *paTable;  /* パレットテーブルへのポインタ */

/* パレットテーブルを作る */
int *make_palet_table(int mxCnt);

/* 24ビットのBMPファイルを出力する */
int output_bmp(char *bmpName, int *paTable, int *cntTable, int width, int height);

/* 3次のジュリア集合の素となる3次関数の定義 */
/* f(z)=z^3+A*z^2+B*z+C, df_f(z)=3*z^2+2*A*z+B */
/* マクロ版と関数版のどちらか速い方を選択する */
#if 0
#define f(z) cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), A)), cmul((z), B)), C)
#define df_f(z) cadd(cadd(cremul(cpow2(z), 3.0), cremul(cmul((z), A), 2.0)), B)
#else
complex f(complex z)
{
  return cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), A)), cmul(z, B)), C);
}
complex df_f(complex z)
{
  return cadd(cadd(cremul(cpow2(z), 3.0), cremul(cmul(z, A), 2.0)), B);
}
#endif

/* 初期値cでf(z)に対するNewton法の漸化式の反復回数を返す */
int iteration(complex c, int limit_cnt)
{
  complex z;
  int cnt;

  if (setjmp(cjmpenv)) {
    return limit_cnt;  /* 計算できなかった */
  }

  /* 漸化式の反復 */
  cnt = limit_cnt;
  z = Z;
  C = c;
  while (cnt--) {
    if (cabs2(f(z)) <= rot_f_limit2) {
      break;
    }
    z = csub(z, cdiv(f(z), df_f(z)));  /* f(z)に対するNewton法の漸化式 rec_f(z)=z-f(z)/df_f(z) */
  }
  return limit_cnt - (cnt + 1);
}

/* 描画ルーチンの本体 */
void draw_body()
{
  complex c0 = tocomplex(reMin, imMin), c;
  double step = scLen / (double)scCnt;
  int *p = cntTable;
  int reCnt, imCnt;

  /* 虚軸方向のループ */
  imCnt = scCnt;
  while (imCnt--) {
    c = c0;

    /* 実軸方向のループ */
    reCnt = scCnt;
    while (reCnt--) {
      *p++ = iteration(c, mxCnt);  /* 漸化式の反復回数を数えて記録する */

      Re(c) += step;  /* 実軸方向に進む */
    }

    Im(c0) += step;  /* 虚軸方向に進む */
  }
}

int main(int argc, char *argv[])
{
  clock_t tm0, tm1;

  /* 係数と初期値のデフォルトを設定 */
  A = tocomplex(Re_A, Im_A);
  B = tocomplex(Re_B, Im_B);
  Z = tocomplex(Re_Z, Im_Z);

  /* コマンドラインのパラメータを取得 */
  switch (argc) {
  case 14:
    sscanf(argv[13], "%lf", &Im(Z));
    sscanf(argv[12], "%lf", &Re(Z));
  case 12:
    sscanf(argv[11], "%lf", &Im(B));
    sscanf(argv[10], "%lf", &Re(B));
  case 10:
    sscanf(argv[9], "%lf", &Im(A));
    sscanf(argv[8], "%lf", &Re(A));
  case 8:
    sscanf(argv[7], "%lf", &rot_f_limit);
  case 7:
    sscanf(argv[6], "%d", &mxCnt);
  case 6:
    sscanf(argv[5], "%d", &scCnt);
  case 5:
    sscanf(argv[4], "%lf", &scLen);
    sscanf(argv[3], "%lf", &imMin);
    sscanf(argv[2], "%lf", &reMin);
  case 2:
    bmpName = argv[1];
    if (bmpName[0] != '-') {
      break;
    }
  default:
    fprintf(stderr, USAGE);
    return 1;
  }

  /* 収束半径を2乗しておく */
  rot_f_limit2 = rot_f_limit * rot_f_limit;

  /* 結果を記録する領域を用意する */
  if ((cntTable = (int *)malloc(sizeof(int) * scCnt * scCnt)) == NULL) {
    fprintf(stderr, "メモリが不足しています\n");
    return 1;
  }

  /* 計測開始 */
  {
    clock_t t = clock();
    while ((tm0 = clock()) == t);
  }

  /* 描画ルーチンの本体 */
  draw_body();

  /* 計測終了 */
  tm1 = clock();

  /* 所要時間を表示 */
  printf("所要時間は %.8g 秒です.\n", ((double)(tm1 - tm0)) / CLK_TCK);

  /* パレットテーブルを作る */
  if ((paTable = make_palet_table(mxCnt)) == NULL) {
    return 1;
  }

  /* BMPファイルを出力する */
  if (output_bmp(bmpName, paTable, cntTable, scCnt, scCnt)) {
    return 1;
  }

  return 0;
}


j3cxbmp.c

説明
 j3cxbmp.cは、3次関数のジュリア集合のBMPデータを生成するプログラムを拡張して、f(z)とf'(z)を別々に指定できるようにしたものです。「ジュリア集合でないもの」を描画することができます。complex.hが必要です。output.cと一緒にgccでコンパイルします。

j3cxbmp.c
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
#include <stdio.h>
#include <setjmp.h>
#include <time.h>
#include "complex.h"

#define USAGE "\
3次関数のジュリア集合の定数項の画像をBMPファイルに出力します(拡張付き)\n\
>j3cxbmp FN SR SI SL PIX CNT LIM F2R F2I F1R F1I ZR ZI D2R D2I D1R D1I DCR DCI\n\
  BMPファイル  FN=BMPファイル名\n\
  定数項の範囲  SR=実部最小[-2.0]  SI=虚部最小[-2.0]  SL=辺の長さ[4.0]\n\
  ピクセル数  PIX=縦横のピクセル数[256]\n\
  反復回数  CNT=Newton法の反復回数の上限[256]\n\
  収束判定  LIM=収束半径(収束したと判断する|f(z)|の上限)[0.01]\n\
  f(z)の2次の係数  F2R=実部[0.0]  F2I=虚部[0.0]\n\
  f(z)の1次の係数  F1R=実部[-1.0]  F1I=虚部[0.0]\n\
  Newton法の初期値  ZR=実部[0.0]  ZI=虚部[0.0]\n\
  f'(z)の2次の係数  D2R=実部[3.0]  D2I=虚部[0.0]\n\
  f'(z)の1次の係数  D1R=実部[2*F2R]  D1I=虚部[2*F2I]\n\
  f'(z)の定数項  DCR=実部[F1R]  DCI=虚部[F1I]\n\
"

/* 描画範囲のデフォルト */
#define SR (-2.0)
#define SI (-2.0)
#define SL (4.0)

/* 縦横のピクセル数のデフォルト */
#define PIX 256

/* 回数制限のデフォルト */
#define CNT 256

/* 収束半径のデフォルト */
#define LIM (0.01)

/* 係数と初期値のデフォルト */
#define F2R (0.0)
#define F2I (0.0)
#define F1R (-1.0)
#define F1I (0.0)
#define ZR (0.0)
#define ZI (0.0)

char *bmpName;  /* BMPファイル名の先頭 */

double reMin = SR;  /* 実部最小 */
double imMin= SI;  /* 虚部最小 */
double scLen = SL;  /* 辺の長さ */
int scCnt = PIX;  /* 縦横のピクセル数(1〜32767) */
int mxCnt = CNT;  /* 回数制限(1〜1000000) */

double rot_f_limit = LIM;  /* 収束半径 */
double rot_f_limit2;  /* 収束半径の2乗 */

complex F2, F1, FC, Z, D2, D1, DC;  /* 係数と定数項と初期値 */

int *cntTable;  /* 結果を記録する領域へのポインタ */
int *paTable;  /* パレットテーブルへのポインタ */

/* パレットテーブルを作る */
int *make_palet_table(int mxCnt);

/* 24ビットのBMPファイルを出力する */
int output_bmp(char *bmpName, int *paTable, int *cntTable, int width, int height);

/* 3次のジュリア集合の素となる3次関数の定義 */
/* f(z)=z^3+A*z^2+B*z+C, df_f(z)=3*z^2+2*A*z+B */
/* マクロ版と関数版のどちらか速い方を選択する */
#if 0
#define f(z) cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), F2)), cmul((z), F1)), FC)
#define df_f(z) cadd(cadd(cmul(cpow2(z), D2), cmul((z), D1)), DC)
#else
complex f(complex z)
{
  return cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), F2)), cmul(z, F1)), FC);
}
complex df_f(complex z)
{
  return cadd(cadd(cmul(cpow2(z), D2), cmul(z, D1)), DC);
}
#endif

/* 初期値cでf(z)に対するNewton法の漸化式の反復回数を返す */
int iteration(complex c, int limit_cnt)
{
  complex z;
  int cnt;

  if (setjmp(cjmpenv)) {
    return limit_cnt;  /* 計算できなかった */
  }

  /* 漸化式の反復 */
  cnt = limit_cnt;
  z = Z;
  FC = c;
  while (cnt--) {
    if (cabs2(f(z)) <= rot_f_limit2) {
      break;
    }
    z = csub(z, cdiv(f(z), df_f(z)));  /* f(z)に対するNewton法の漸化式 rec_f(z)=z-f(z)/df_f(z) */
  }
  return limit_cnt - (cnt + 1);
}

/* 描画ルーチンの本体 */
void draw_body()
{
  complex c0 = tocomplex(reMin, imMin), c;
  double step = scLen / (double)scCnt;
  int *p = cntTable;
  int reCnt, imCnt;

  /* 虚軸方向のループ */
  imCnt = scCnt;
  while (imCnt--) {
    c = c0;

    /* 実軸方向のループ */
    reCnt = scCnt;
    while (reCnt--) {
      *p++ = iteration(c, mxCnt);  /* 漸化式の反復回数を数えて記録する */

      Re(c) += step;  /* 実軸方向に進む */
    }

    Im(c0) += step;  /* 虚軸方向に進む */
  }
}

int main(int argc, char *argv[])
{
  clock_t tm0, tm1;

  /* 係数と初期値のデフォルトを設定 */
  F2 = tocomplex(F2R, F2I);
  F1 = tocomplex(F1R, F1I);
  Z = tocomplex(ZR, ZI);

  /* コマンドラインのパラメータを取得 */
  switch (argc) {
  case 20:
    sscanf(argv[19], "%lf", &Im(DC));
    sscanf(argv[18], "%lf", &Re(DC));
    sscanf(argv[17], "%lf", &Im(D1));
    sscanf(argv[16], "%lf", &Re(D1));
    sscanf(argv[15], "%lf", &Im(D2));
    sscanf(argv[14], "%lf", &Re(D2));
  case 14:
    sscanf(argv[13], "%lf", &Im(Z));
    sscanf(argv[12], "%lf", &Re(Z));
  case 12:
    sscanf(argv[11], "%lf", &Im(F1));
    sscanf(argv[10], "%lf", &Re(F1));
  case 10:
    sscanf(argv[9], "%lf", &Im(F2));
    sscanf(argv[8], "%lf", &Re(F2));
  case 8:
    sscanf(argv[7], "%lf", &rot_f_limit);
  case 7:
    sscanf(argv[6], "%d", &mxCnt);
  case 6:
    sscanf(argv[5], "%d", &scCnt);
  case 5:
    sscanf(argv[4], "%lf", &scLen);
    sscanf(argv[3], "%lf", &imMin);
    sscanf(argv[2], "%lf", &reMin);
  case 2:
    bmpName = argv[1];
    if (bmpName[0] != '-') {
      break;
    }
  default:
    fprintf(stderr, USAGE);
    return 1;
  }

  /* 収束半径を2乗しておく */
  rot_f_limit2 = rot_f_limit * rot_f_limit;

  /* f'(z)の係数を決定する */
  if (argc <= 14) {
    D2 = tocomplex(3.0, 0.0);
    D1 = cremul(F2, 2.0);
    DC = F1;
  }

  /* 結果を記録する領域を用意する */
  if ((cntTable = (int *)malloc(sizeof(int) * scCnt * scCnt)) == NULL) {
    fprintf(stderr, "メモリが不足しています\n");
    return 1;
  }

  /* 計測開始 */
  {
    clock_t t = clock();
    while ((tm0 = clock()) == t);
  }

  /* 描画ルーチンの本体 */
  draw_body();

  /* 計測終了 */
  tm1 = clock();

  /* 所要時間を表示 */
  printf("所要時間は %.8g 秒です.\n", ((double)(tm1 - tm0)) / CLK_TCK);

  /* パレットテーブルを作る */
  if ((paTable = make_palet_table(mxCnt)) == NULL) {
    return 1;
  }

  /* BMPファイルを出力する */
  if (output_bmp(bmpName, paTable, cntTable, scCnt, scCnt)) {
    return 1;
  }

  return 0;
}


complex.h

説明
 complex.hは、複素数演算マクロ定義ファイルです。

complex.h
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
/* 複素数演算マクロ集 */
/* 2000.05.21 by M.Kamada */

#include <math.h>
#include <setjmp.h>

/* 演算を続行できなくなったときのためにsetjmp(cjmpenv)で環境を準備しておくこと */
jmp_buf cjmpenv;

/* setjmp(cjmpenv)の返却値 */
#define C_DIVIDE_BY_ZERO 1  /* 0で除算しようとした */
#define C_ARG_ZERO 2  /* 0の角度を求めようとした */
#define C_LOG_ZERO 3  /* log(0)を求めようとした */

#ifndef __GNUC__
#error "Sorry, GCC only"
#endif

#if __GNUC__==2
#define complex __complex__
#define Re(z) (__real__ (z))
#define Im(z) (__imag__ (z))
#else
typedef struct {
  double re;
  double im;
} complex;
#define Re(z) ((z).re)
#define Im(z) ((z).im)
#endif

#define tocomplex(r,i) \
({ \
  complex _z_; \
  Re(_z_) = (r); \
  Im(_z_) = (i); \
  _z_; \
})

#if __GNUC__==2
#define C_1I (1.0i)
#define C_2I (2.0i)
#define C_3I (3.0i)
#define C_4I (4.0i)
#define C_5I (5.0i)
#define C_6I (6.0i)
#else
#define C_1I tocomplex(0.0, 1.0)
#define C_2I tocomplex(0.0, 2.0)
#define C_3I tocomplex(0.0, 3.0)
#define C_4I tocomplex(0.0, 4.0)
#define C_5I tocomplex(0.0, 5.0)
#define C_6I tocomplex(0.0, 6.0)
#endif

#define C_2PI (6.28318530717958647692528676656)  /* 2*pi */
#define C_PI (3.14159265358979323846264338328)  /* pi */
#define C_PI_2 (1.57079632679489661923132169164)  /* pi/2 */
#define C_PI_4 (0.78539816339744830961566084582)  /* pi/4 */
#define C_1_PI (0.318309886183790671537767526745)  /* 1/pi */
#define C_2_PI (0.63661977236758134307553505349)  /* 2/pi */
#define C_2_SQRTPI (1.12837916709551257389615890312)  /* 2/sqrt(pi) */
#define C_E (2.71828182845904523536028747135)  /* e */
#define C_LN10 (2.30258509299404568401799145468)  /* ln(10) */
#define C_LN2 (0.693147180559945309417232121458)  /* ln(2) */
#define C_LOG10E (0.434294481903251827651128918917)  /* log10(e) */
#define C_LOG2E (1.442695040888963407359924681)  /* log2(e) */
#define C_SQRT1_2 (0.707106781186547524400844362105)  /* 1/sqrt(2) */
#define C_SQRT2 (1.41421356237309504880168872421)  /* sqrt(2) */
#define C_SQRT3 (1.73205080756887729352744634151)  /* sqrt(3) */
#define C_SQRT5 (2.23606797749978969640917366873)  /* sqrt(5) */
#define C_SQRT6 (2.44948974278317809819728407471)  /* sqrt(6) */
#define C_SQRT7 (2.64575131106459059050161575364)  /* sqrt(7) */
#define C_SQRT8 (2.82842712474619009760337744842)  /* sqrt(8) */
#define C_SQRT10 (3.16227766016837933199889354443)  /* sqrt(10) */

/* xの実数部にsを加える */
#if __GNUC__==2
#define creadd(x,s) ((x)+(s))
#else
#define creadd(x,s) \
({ \
  complex _z_ = (x); \
  Re(_z_) += (s); \
  _z_; \
})
#endif

/* xの虚数部にsを加える */
#define cimadd(x,s) \
({ \
  complex _z_ = (x); \
  Im(_z_) += (s); \
  _z_; \
})

/* xの実数部からsを引く */
#if __GNUC__==2
#define cresub(x,s) ((x)-(s))
#else
#define cresub(x,s) \
({ \
  complex _z_ = (x); \
  Re(_z_) -= (s); \
  _z_; \
})
#endif

/* xの虚数部からsを引く */
#define cimsub(x,s) \
({ \
  complex _z_ = (x); \
  Im(_z_) -= (s); \
  _z_; \
})

/* xの実数s倍 */
#if __GNUC__==2
#define cremul(x,s) ((x)*(s))
#else
#define cremul(x,s) \
({ \
  complex _z_ = (x); \
  double _s_ = (s); \
  Re(_z_) *= _s_; \
  Im(_z_) *= _s_; \
  _z_; \
})
#endif

/* xの虚数s*i倍 */
#define cimmul(x,s) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  double _s_ = (s); \
  Re(_z_) = Im(_x_) * (-_s_); \
  Im(_z_) = Re(_x_) * _s_; \
  _z_; \
})

/* xのスカラーs倍 */
#if __GNUC__==2
#define cscal(x,s) ((x)*(s))
#else
#define cscal(x,s) \
({ \
  complex _z_ = (x); \
  double _s_ = (s); \
  Re(_z_) *= _s_; \
  Im(_z_) *= _s_; \
  _z_; \
})
#endif

/* xの共役複素数 */
#if __GNUC__==2
#define conj(x) (~((__complex__)(x)))
#else
#define conj(x) \
({ \
  complex _z_ = (x); \
  Im(_z_) = -Im(_z_); \
  _z_; \
})
#endif

/* xの絶対値の2乗 */
#define cabs2(x) \
({ \
  complex _x_ = (x); \
  Re(_x_) * Re(_x_) + Im(_x_) * Im(_x_); \
})

/* xの絶対値 */
#define cabs(x) \
({ \
  complex _x_ = (x); \
  sqrt(Re(_x_) * Re(_x_) + Im(_x_) * Im(_x_)); \
})

/* xの角度 */
#define carg(x) \
({ \
  complex _x_ = (x); \
  if (Re(_x_) == 0.0 && Im(_x_) == 0.0) { \
    longjmp(cjmpenv, C_ARG_ZERO); \
  } \
  atan2(Im(_x_), Re(_x_)); \
})

/* x+y */
#if __GNUC__==2
#define cadd(x,y) ((x)+(y))
#else
#define cadd(x,y) \
({ \
  complex _z_ = (x); \
  complex _y_ = (y); \
  Re(_z_) += Re(_y_); \
  Im(_z_) += Im(_y_); \
  _z_; \
})
#endif

/* x-y */
#if __GNUC__==2
#define csub(x,y) ((x)-(y))
#else
#define csub(x,y) \
({ \
  complex _z_ = (x); \
  complex _y_ = (y); \
  Re(_z_) -= Re(_y_); \
  Im(_z_) -= Im(_y_); \
  _z_; \
})
#endif

/* x*y */
#if __GNUC__==2
#define cmul(x,y) ((x)*(y))
#else
#define cmul(x,y) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  complex _y_ = (y); \
  Re(_z_) = Re(_x_) * Re(_y_) - Im(_x_) * Im(_y_); \
  Im(_z_) = Re(_x_) * Im(_y_) + Im(_x_) * Re(_y_); \
  _z_; \
})
#endif

/* x/y */
#define cdiv(x,y) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  complex _y_ = (y); \
  double _r2_; \
  if ((_r2_ = cabs2(_y_)) == 0.0) { \
    longjmp(cjmpenv, C_DIVIDE_BY_ZERO); \
  } \
  Re(_z_) = (Re(_x_) * Re(_y_) + Im(_x_) * Im(_y_)) / _r2_; \
  Im(_z_) = (Im(_x_) * Re(_y_) - Re(_x_) * Im(_y_)) / _r2_; \
  _z_; \
})

/* xの2乗 (a+bi)^2=aa+2abi-bb */
#define cpow2(x) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  Re(_z_) = Re(_x_) * Re(_x_) - Im(_x_) * Im(_x_); \
  Im(_z_) = (Re(_x_) * Im(_x_)) * 2.0; \
  _z_; \
})

/* xの3乗 (a+bi)^3=aaa+3aabi-3abb-bbbi */
#define cpow3(x) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  double _aa_, _bb_; \
  _aa_ = Re(_x_) * Re(_x_); \
  _bb_ = Im(_x_) * Im(_x_); \
  Re(_z_) = Re(_x_) * (_aa_ - 3.0 * _bb_); /* aaa-3abb */ \
  Im(_z_) = Im(_x_) * (3.0 * _aa_ - _bb_); /* 3aab-bbb */ \
  _z_; \
})

/* xの4乗 (a+bi)^4=aaaa+4aaabi-6aabb-4abbbi+bbbb */
#define cpow4(x) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  double _aa_, _bb_, _ab_; \
  _aa_ = Re(_x_) * Re(_x_); \
  _bb_ = Im(_x_) * Im(_x_); \
  _ab_ = Re(_x_) * Im(_x_); \
  Re(_z_) = _aa_ * (_aa_ - 6.0 * _bb_) + _bb_ * _bb_; /* aaaa-6aabb+bbbb */ \
  Im(_z_) = 4.0 * _ab_ * (_aa_ - _bb_); /* 4aaab-4abbb */ \
  _z_; \
})

/* xの5乗 (a+bi)^5=aaaaa+5aaaabi-10aaabb-10aabbbi+5abbbb+bbbbbi */
#define cpow5(x) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  double _aa_, _bb_, _ab_, _aaa_, _bbb_; \
  _aa_ = Re(_x_) * Re(_x_); \
  _bb_ = Im(_x_) * Im(_x_); \
  _ab_ = Re(_x_) * Im(_x_); \
  _aaa_ = _aa_ * Re(_x_); \
  _bbb_ = _bb_ * Im(_x_); \
  Re(_z_) = _aaa_ * (_aa_ - 10.0 * _bb_) + 5.0 * _ab_ * _bbb_; /* aaaaa-10aaabb+5abbbb */ \
  Im(_z_) = 5.0 * _aaa_ * _ab_ + _bbb_ * (_bb_ - 10.0 * _aa_); /* 5aaaab+bbbbb-10aabbb */ \
  _z_; \
})

/* xの6乗 (a+bi)^6=aaaaaa+6aaaaabi-15aaaabb-20aaabbbi+15aabbbb+6abbbbbi-bbbbbb */
#define cpow6(x) \
({ \
  complex _z_; \
  complex _x_ = (x); \
  double _aa_, _bb_, _ab_, _aaaa_, _bbbb_; \
  _aa_ = Re(_x_) * Re(_x_); \
  _bb_ = Im(_x_) * Im(_x_); \
  _ab_ = Re(_x_) * Im(_x_); \
  _aaaa_ = _aa_ * _aa_; \
  _bbbb_ = _bb_ * _bb_; \
  Re(_z_) = _aa_ * (_aa_ * (_aa_ - 15.0 * _bb_) + 15.0 * _bbbb_) - _bb_ * _bbbb_; /* aaaaaa-15aaaabb+15aabbbb-bbbbbb */ \
  Im(_z_) = _ab_ * (6.0 * (_aaaa_ + _bbbb_) - 20.0 * _aa_ * _bb_); /* 6aaaaab+6abbbbb-20aaabbb */ \
  _z_; \
})

/* xのe乗 */
#define cpow(x,e) \
({ \
  complex _z_ = (x); \
  double _e_ = (e); \
  double _r_, _t_; \
  _r_ = pow(cabs(_z_), _e_); \
  _t_ = carg(_z_) * _e_; \
  Re(_z_) = _r_ * cos(_t_); \
  Im(_z_) = _r_ * sin(_t_); \
  _z_; \
})

/* xのm番目のn乗根 */
#define crot(x,n,m) \
({ \
  complex _z_ = (x); \
  int _n_ = (n); \
  double _r_, _t_; \
  _r_ = pow(cabs(_z_), 1.0 / (double)_n_); \
  _t_ = carg(_z_) + C_2PI * (double)(m) / (double)_n_; \
  Re(_z_) = _r_ * cos(_t_); \
  Im(_z_) = _r_ * sin(_t_); \
  _z_; \
})

/* 距離の2乗 */
#define cdist2(x,y) \
({ \
  complex _x_ = (x); \
  complex _y_ = (y); \
  double _r_ = Re(_x_) - Re(_y_); \
  double _i_ = Im(_x_) - Im(_y_); \
  _r_ * _r_ + _i_ * _i_; \
})

/* 距離 */
#define cdist(x,y) \
({ \
  complex _x_ = (x); \
  complex _y_ = (y); \
  double _r_ = Re(_x_) - Re(_y_); \
  double _i_ = Im(_x_) - Im(_y_); \
  sqrt(_r_ * _r_ + _i_ * _i_); \
})

/* 複素数zの複素数c乗 */
/* z^(a+i*b)
     = (abs(z)*e^(i*arg(z)))^(a+i*b)
     = abs(z)^(a+i*b)*e^(i*(a+i*b)*arg(z))
     = abs(z)^a*(e^(i*b*log(abs(z))))*e^(i*a*arg(z))*e^(-b*arg(z))
     = abs(z)^a*e^(-b*arg(z))*e^(i*(a*arg(z)+b*log(abs(z))))
     = abs(z)^a*e^(-b*arg(z))*(cos(a*arg(z)+b*log(abs(z)))+i*sin(a*arg(z)+b*log(abs(z)))) */
#define cpowc(z,c) \
({ \
  complex _z_ = (z); \
  complex _c_ = (c); \
  double _a_, _b_; \
  double _abs_z_, _arg_z_; \
  double _r_, _t_; \
  _abs_z_ = cabs(_z_); \
  if (_abs_z_ == 0.0) { \
    if (Re(_c_) == 0.0) { \
      longjmp(cjmpenv, C_LOG_ZERO); \
    } \
    _z_ = 0.0; \
  } else { \
    _a_ = Re(_c_); \
    _b_ = Im(_c_); \
    _arg_z_ = carg(_z_); \
    _r_ = pow(_abs_z_, _a_) * pow(C_E, -_b_*_arg_z_); \
    _t_ = _a_*_arg_z_ + _b_*log(_abs_z_); \
    Re(_z_) = _r_ * cos(_t_); \
    Im(_z_) = _r_ * sin(_t_); \
  } \
  _z_; \
})


output.c

説明
 output.cは、マンデルブロ集合やジュリア集合の演算結果に彩色を施してBMPファイルを生成するプログラムです。

output.c
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#include <stdio.h>
#include <math.h>

#define C_PI (3.14159265358979323846264338328)  /* pi */

#define RGB(r,g,b) (((r) << 16) | ((g) << 8) | (b))

/*
  int hsvtorgb(int h, int s, int v);

    h = 0〜1535  色相
        0〜255   赤〜黄色
      256〜511   黄色〜緑
      512〜767   緑〜シアン
      768〜1023  シアン〜青
     1024〜1279  青〜マゼンダ
     1280〜1535  マゼンダ〜赤

    s = 0〜255  彩度(白レベル)
        0〜255  白〜原色

    v = 0〜255  明度(黒レベル)
        0〜255  黒〜原色

  hsvtorgb() = (r << 16) | (g << 8) | b
    r = 0〜255  赤の成分
    g = 0〜255  緑の成分
    b = 0〜255  青の成分
*/

int hsvtorgb(int h, int s, int v)
{
  int f, f1, s1, v1, w, d;

  if (s == 0) {
    return RGB(v, v, v);
  }

  f = h & 255;
  f1 = (f << 1) | 1;
  s1 = (s << 1) | 1;
  v1 = (v << 1) | 1;
  w = v - ((s1 * v1 + (1 << 9)) >> 10);
  d = (f1 * s1 * v1 + (1 << 18)) >> 19;
  switch (h >> 8) {
  case 0: return RGB(v, w+d, w);
  case 1: return RGB(v-d, v, w);
  case 2: return RGB(w, v, w+d);
  case 3: return RGB(w, v-d, v);
  case 4: return RGB(w+d, w, v);
  case 5: return RGB(v, w, v-d);
  }
}

/* パレットテーブルを作る */
/* 返却値: パレットテーブルの先頭,NULL=メモリが不足している */
int *make_palet_table(int mxCnt)
{
  int *paTable, *p;
  int n;

  if ((paTable = (int *)malloc((mxCnt + 1) * sizeof(int))) == NULL) {
    fprintf(stderr, "メモリが不足しています\n");
    return NULL;
  }

  p = paTable;

  for (n = 0; n < mxCnt; n++) {
    double t = n / (mxCnt / 15.0);
    int u, v;

    u = (int)(cos((t) * C_PI) * 112.0 + 212.0);
    u = (u <= 255) ? u : 255;
    v = (int)(cos((t + 1.0) * C_PI) * 72.0 + 188.0);
    v = (v <= 255) ? v : 255;
    *p++ = hsvtorgb((1536 * n) / mxCnt, u, v);
  }
  *p = 0;

  return paTable;
}

/* BMPファイルのヘッダ */
char bmpHeader[54] = {
  'B', 'M',          /* [ 0] ファイルタイプ */
  0, 0, 0, 0,        /* [ 2] ファイルサイズ */
  0, 0, 0, 0,        /* [ 6] 予約 */
  0x36, 0, 0, 0,     /* [10] ビットマップデータのシーク位置 */
  0x28, 0, 0, 0,     /* [14] ここから始まるヘッダの長さ */
  0, 0, 0, 0,        /* [18] ビットマップの幅 */
  0, 0, 0, 0,        /* [22] ビットマップの高さ */
  0x01, 0,           /* [26] プレーン数 */
  0x18, 0,           /* [28] 1ピクセルあたりのビット数(1,4,8,24) */
  0, 0, 0, 0,        /* [30] 圧縮タイプ(0=非圧縮) */
  0, 0, 0, 0,        /* [34] ビットマップデータの長さ */
  0x20, 0x2e, 0, 0,  /* [38] 水平解像度(px/m) */
  0x20, 0x2e, 0, 0,  /* [42] 垂直解像度(px/m) */
  0, 0, 0, 0,        /* [46] カラーインデックス数 */
  0, 0, 0, 0         /* [50] 重要なカラーインデックス数 */
};                   /* [54] */

/* 32ビットの数値をlittle-endianで書き込む */
void wl4(char *p, int value)
{
  *p++ =  value        & 0xff;
  *p++ = (value >>  8) & 0xff;
  *p++ = (value >> 16) & 0xff;
  *p++ = (value >> 24) & 0xff;
}

/* 24ビットのBMPファイルを出力する */
/* 返却値:0=正常終了,1=異常終了 */
int output_bmp(char *bmpName, int *paTable, int *cntTable, int width, int height)
{
  FILE *fp;
  int *p;
  int n;

  /* BMPファイルのヘッダを作る */
  wl4(&bmpHeader[ 2], sizeof(bmpHeader) + 3 * width * height);
  wl4(&bmpHeader[34],                     3 * width * height);
  wl4(&bmpHeader[18], width);
  wl4(&bmpHeader[22], height);

  /* BMPファイルを作る */
  if ((fp = fopen(bmpName, "wb")) == NULL) {
    fprintf(stderr, "ファイルが作れません … %s\n", bmpName);
    return 1;
  }

  /* BMPファイルのヘッダを出力する */
  for (n = 0; n < sizeof(bmpHeader); n++) {
    fputc(bmpHeader[n], fp);
  }

  /* ビットマップデータを出力する */
  p = cntTable;
  n = width * height;
  while (n--) {
    int rgb = paTable[*p++];
    fputc(rgb & 255, fp);          /* B */
    fputc((rgb >> 8) & 255, fp);   /* G */
    fputc((rgb >> 16) & 255, fp);  /* R */
  }

  fclose(fp);

  return 0;
}

counter Since December 18, 2000
フラクタル | サイトマップ | ホーム
E-mail: m_kamada@nifty.comMirrorCopyright (C) 1999-2001 M.Kamada All Rights Reserved.
[PR] | インプラント債務整理 相談債務整理インプラント転職サイトSEOアクセス解析ハウスメーカーレンタルオフィスSEO対策消費者金融不動産担保ローン時計車 買取ハワイハワイ挙式アスクル転職生命保険テンプレート沖縄旅行動画FX免許合宿二輪引越し消費者金融税理士ゴルフ会員権留学レーシックマッサージ貸し店舗FX投資信託くりっく365アフィリエイト育毛剤FXホームページ制作デイトレードFXホノルルマラソンベスト ハワイ ホテル レーツバリ島ハワイウエディングHawaii hotelsHawaii Activitiesbhhr
【運営会社「パラダイムシフト」サービス】 ハワイ現地オプショナルツアーリラックマ.ハワイホテルガイド) - ビジネスクラス航空券 - 格安航空券(1) - 格安航空券(2) - 海外ホテル - 韓国旅行
無料ホームページ作成 - レンタルサーバー - 携帯ホームページ - ブログ - ホテル 予約 - 格安航空券 - 長期滞在 - タイムシェア - ヴィラ - ハワイ コンドミニアム - バリ島 ホテル
[PR] 包茎でお困りの方必見!厳選したクリニックをご紹介!