Ｃ言語における複素数の利用方法

C89 の複素数

C89（ANSI C）等の古い標準Ｃでは，複素数用のデータ型や関数を各プログラマが独自に実装していた．

複素数型の定義例：

#include <math.h>

// 複素数のデータ型
typedef struct { double	re, im; } Complex;

// 直交座標形式の数値設定関数
//	数学：z = x + i y	コード：z = cValue(x, y);
Complex cValue(double x, double y)
{
	Complex	z;
	z.re = x;
	z.im = y;
	return (z);
}

// 極座標形式の数値設定関数
//	数学：z = r exp(j a)	コード：z = cPolar(r, a);
Complex cPolar(double r, double a)
{
	return (cValue(r*cos(a), r*sin(a)));
}

// 加算関数
//	数学：c = a + b		コード：c = cAdd(a, b);
Complex cAdd(Complex a, Complex b)
{
	return (cValue(a.re + b.re, a.im + b.im));
}
...

C99 の複素数

組込機能

C99 標準（や GCC）では，実数の拡張版として，複素数を取り扱えるようになった．もちろん，複素数同士の四則演算や変数への代入も可能である．

使用例：

名称	コード書式	コード例	補足
複素数型	`double _Complex`	`double _Complex z;`	複素変数 `z` の宣言
虚数接尾辞	`i`	`z = 1.0 + 2.0i;`	変数 `z` への複素数 $ 1 + i 2 $ の代入

型名にアンダーバーが付けられているのは，古い独自実装との互換性確保（名前衝突回避）のための配慮でしょう．

これで，複素数の独自実装は不要となった．ただし，さらに新しい C11 標準では，必須機能から外れてしまったので，独自実装が再度必要な環境もあり得る．互換性や規模感などの問題で揉めたのでしょう．

標準ライブラリ

複素数用の数学関数等も標準ライブラリとして提供されるようになった．

ヘッダファイル：complex.h

利用方法：

#include <complex.h>

複素数のデータ型：double complex
アンダーバーが不要になった．しかし，double は必要なので，まだ長たらしい...
虚数単位のマクロ定数：I

互換性問題が頻出しそうな短さ... 工学分野では記号 $ i $ や $ I $ は電流とか整数を表すものだし．

使用例：

名称	数学的表記例	コード例
直交座標形式	$ \dot{z} = x + i y $	`z = x + I*y;`
実部	$ x = \RE{\dot{z}} $	`x = creal(z);`
虚部	$ y = \IM{\dot{z}} $	`y = cimag(z);`
極座標形式	$ \dot{z} = r \exp{j a} $	`z = rcexp(Ia);`
絶対値	$ r = \N{\dot{z}} $	`r = cabs(z);`
偏角	$ a = \arg{\dot{z}} $	`a = carg(z);`
複素共役	$ \dot{z}_2 = \overline{\dot{z}_1} $	`z2 = conj(z1);`

ライブラリファイル：libm.so
コンパイル方法：
```
$ cc  ソース.c -std=c99  -o プログラム  -lm
```
オンラインマニュアル：complex
```
$ man 7 complex
```

参考文献

C, C++ で複素数／桂田先生＠明治大学??

C99 の魅惑の機能

C99 では，複素数の他にも便利（かつ危険）な機能が提供されている．

任意位置の変数宣言：ブロックの先頭以外でも変数を宣言可能．
C89 では，変数宣言の記述位置はブロックの先頭部分に限定されていた．どこにでも書き散らかしてよくなった？わけではない．

どの変数をどの範囲内で利用するのか把握し易くするため，宣言位置は可能な限り集結すべきでしょう．
可変長配列：配列の要素数を実行時にも指定可能．
C89 では，配列サイズはコンパイル時に指定される固定長なので，可変長配列にはメモリ確保関数 malloc() とかを使う必要があった． C99 では，それが不要になった？わけではない．

試作・実験レベルでの利用にとどめましょう．この可変長配列は，お手軽ではあるが，メモリ不足になると実行中のプログラムが突然死するので危険．エラーチェックできないので，安全性が必要な場合には使えませんよ．という訳で，その後の C11 標準では必須機能から外されたのでしょう．

コード例：

...
int main(void)
{
	printf("要素数 > ");
	int	n;			// 任意位置
	scanf("%d", &n);
	int	x[n];		// 任意位置，動的配列
	for (int i = 0; i < n; i++) {	// 任意位置
		x[i] = i;
	}
	...
}

C99 複素数のカスタマイズ

工学分野では，虚数単位の記号として $ i $ とか I ではなく $ j $ とか J を使いたい．また，複素数利用時には，しばしば，円周率πも使いたい．そして，型名の double complex は長すぎる...

カスタムヘッダファイルの定義例：cx.h

#ifndef	CX_H
#define	CX_H

#ifdef __cplusplus
extern "C" {
#endif

// 円周率マクロ定数 M_PI を利用可能に
#define	_USE_MATH_DEFINES	// for M_PI on Windows
#define	__USE_XOPEN		// for M_PI on GCC/Linux
#include <math.h>

// よく利用する複素数計算を簡便に
#include <complex.h>
typedef	double complex	Cx;	// 複素数型：Cx
#define	J	1.0i		// 虚数単位jのマクロ定数：J
#define	J2Pi	2.0i*M_PI	// j 2πのマクロ定数：J2Pi
#define	Pi2	2.0*M_PI	// 2πのマクロ定数：Pi2
#define	ExpJ(rad)	cexp(J*(rad))	// exp(j a)の引数付きマクロ：ExpJ(a)
#define	Rot(N, kn)	cexp(-J2Pi*(kn)/(N))	// 回転子W_N^{kn}の引数付きマクロ：Rot(N, kn)

#ifdef __cplusplus
}
#endif

#endif

カスタムヘッダファイルの利用例：

ソース：osc.c

// RL直列交流回路の電圧・電流のシミュレーションと
// 	オシロスコープ（oscilloscope）的なグラフ表示
// コンパイル：$ cc osc.c -std=c99 -lm -o osc
// 実行例：
// 	$ ./osc
// 	$ ./osc 1
// 	$ ./osc > osc.txt; ./osc.plot osc.txt &

#include <stdio.h>
#include "cx.h"

#define	GMAX	20	// グラフ表示の最大振幅

void plot(double t, double v, double i, double max, int graph)
{
	printf("%6.3f\t%6.2f\t%6.2f\t", t, v, i);
	if (graph) {
		int	pv = (int)(GMAX*v/max);
		int	pi = (int)(GMAX*i/max);

		for (int k = -GMAX; k <= GMAX; k++) {
			int	c;
			if (k == pi) c = 'I';
			else if (k == pv) c = 'V';
			else if (k == 0) c = '.';
			else c = ' ';
			putchar(c);
		}
	}
	putchar('\n');
}
int main(int argc, char *argv[])
{
	int	graph = 0;	// グラフ表示フラグ
	if (argc > 1) graph = 1;	// コマンド引数あり → グラフ表示へ

	double	v0 = 10.0;	// 電圧振幅 V_0 [V]
	double	f = 50.0;	// 周波数 f [Hz]
	double	r = 50.0;	// 抵抗 R [Ω]
	double	l = 10.0;	// インダクタンス L [H]

	double	w = Pi2*f;		// 角周波数
	Cx	z = (r + J*w*l)/1000.0;	// インピーダンス [kΩ]
	double	dt = 1.0/f/20.0;	// 刻み時間
	double	tmax = 2.0/f;		// 表示時間

	printf("t[s]\tv(t)[V]\ti(t)[mA]\n");
	for (double t = 0.0; t < tmax; t += dt) {	// 時刻 t [s]
		Cx	v = v0*ExpJ(w*t);	// 電圧フェーザ v(t) [V]
		Cx	i = v/z;		// 電流フェーザ i(t) [mA]
		plot(t, creal(v), creal(i), v0, graph);
	}
	return (0);
}

グラフ表示用 gnuplot スクリプト：osc.plot

実行例：

$ chmod +x osc.plot

$ cc osc.c -std=c99  -o osc  -lm

$ ./osc 1		# 文字端末でグラフ表示
t[s]	v(t)[V]	i(t)[mA]
 0.000	 10.00	  0.05	                    I                   V
 0.001	  9.51	  1.03	                    . I                V 
 0.002	  8.09	  1.91	                    .  I            V    
 0.003	  5.88	  2.60	                    .    I     V         
 0.004	  3.09	  3.04	                    .     I              
 0.005	  0.00	  3.18	                    V     I              
 0.006	 -3.09	  3.01	              V     .     I              
 0.007	 -5.88	  2.54	         V          .    I               
 0.008	 -8.09	  1.83	    V               .  I                 
 0.009	 -9.51	  0.94	 V                  .I                   
 0.010	-10.00	 -0.05	V                   I                    
 0.011	 -9.51	 -1.03	 V                I .                    
 0.012	 -8.09	 -1.91	    V            I  .                    
 0.013	 -5.88	 -2.60	         V     I    .                    
 0.014	 -3.09	 -3.04	              I     .                    
 0.015	  0.00	 -3.18	              I     V                    
 0.016	  3.09	 -3.01	              I     .     V              
 0.017	  5.88	 -2.54	               I    .          V         
 0.018	  8.09	 -1.83	                 I  .               V    
 0.019	  9.51	 -0.94	                   I.                  V 
 0.020	 10.00	  0.05	                    I                   V
 0.021	  9.51	  1.03	                    . I                V 
 0.022	  8.09	  1.91	                    .  I            V    
 0.023	  5.88	  2.60	                    .    I     V         
 0.024	  3.09	  3.04	                    .     I              
 0.025	 -0.00	  3.18	                    V     I              
 0.026	 -3.09	  3.01	              V     .     I              
 0.027	 -5.88	  2.54	         V          .    I               
 0.028	 -8.09	  1.83	    V               .  I                 
 0.029	 -9.51	  0.94	 V                  .I                   
 0.030	-10.00	 -0.05	V                   I                    
 0.031	 -9.51	 -1.03	 V                I .                    
 0.032	 -8.09	 -1.91	    V            I  .                    
 0.033	 -5.88	 -2.60	         V     I    .                    
 0.034	 -3.09	 -3.04	              I     .                    
 0.035	  0.00	 -3.18	              I     V                    
 0.036	  3.09	 -3.01	              I     .     V              
 0.037	  5.88	 -2.54	               I    .          V         
 0.038	  8.09	 -1.83	                 I  .               V    
 0.039	  9.51	 -0.94	                   I.                  V 

$ ./osc > osc.txt; ./osc.plot osc.txt &	# gnuplot でグラフ表示

名称	数学的表記例	コード例
直交座標形式	\( \dot{z} = x + i y \)	`z = x + I*y;`
実部	\( x = \RE{\dot{z}} \)	`x = creal(z);`
虚部	\( y = \IM{\dot{z}} \)	`y = cimag(z);`
極座標形式	\( \dot{z} = r \exp{j a} \)	`z = rcexp(Ia);`
絶対値	\( r = \N{\dot{z}} \)	`r = cabs(z);`
偏角	\( a = \arg{\dot{z}} \)	`a = carg(z);`
複素共役	\( \dot{z}_2 = \overline{\dot{z}_1} \)	`z2 = conj(z1);`