代数方程式ライブラリ

［ 簡単な説明 ］

代数方程式解法ルーチンのライブラリです。

代数方程式とは、

a0・xn ＋ a1・xn−1 ＋ a2・xn−2 ＋ ... ＋ an ＝ 0

の形をした方程式のことで、a0〜an を係数、n を代数方程式の次数といいます。
ｎ次の代数方程式は重根まで含めると必ずｎ個の根を持ちます。（ただし、必ずしも実数とは限りません。）

高次の代数方程式を解くルーチン（特に、複数個の根を計算するルーチン）は、１個毎に根を求め、その根を用いて次数が１つ低下した代数方程式を再形成し、次にその根を求めるというアルゴリズムになっています。従って、一般的に後に求める根ほど誤差が大きくなりますので、高次の代数方程式を解く場合には必ず誤差の評価を行って下さい。

なお、求められた根を使って、次数が１次低下した代数方程式の係数を求めるのは、高次の係数から順に計算する方法と、低次の係数から計算する方法がありますが、絶対値が小さい順に根を求めて行く場合（既知の根より次に求める根の方が絶対値が大きい場合）には、高次の係数から計算する方が累積誤差が小さくなります。逆に、絶対値が大きい順に根を求める場合には、低次の係数から計算する方が誤差が小さくなります。（証明略）

また、各関数は、収束判定値と繰返し回数許容値を要求します。各関数は繰返し計算により根の精度を高めていくアルゴリズムをとっているので、繰返し処理を停止させる条件を収束判定値により決めます。
一般的には収束判定値が小さいほど、根の精度は向上しますが、コンピュータの演算は有限桁なので精度向上にも限界があります。そのため、収束条件を満足しないで繰り返す場合の最大値を第２パラメータとして与え、プログラムが無限ループに陥るのを回避する必要があります。

下記の No.2〜5 は実数係数、No.1,6,7 は複素係数の代数方程式を解きます。
特に、No.5 に関しては、実数の範囲に限り、代数方程式以外の方程式も解けます。

なお、複素数の要素は、複素数ライブラリと同じ、Complex 構造体で表しています。
No. 関数名 関数の機能 呼び出し例 説明
1 qurt( ) ２次
代数方程式		 qurt(a, b, c, x); a・x2 ＋ b・x ＋ c ＝ 0 の根を配列 x に
入れて返します。
a, b, c, x はすべて Complex 型です。
2 carda( ) ３次
代数方程式		 カルダノ法 carda(a, x); a0・x3 ＋ a1・x2 ＋ a2・x ＋ a3 ＝ 0 の根を
配列 x に入れて返します。
a は double、x は Complex 型です。
3 newton( ) 高次
代数方程式		 ニュートン法 x = newton(a,
n, eps, iter);		 係数を配列 a、次数を n、収束判定値を eps、
繰返し許容数を iter に入れて呼び出すと、
実根１つを戻り値として返します。
a, eps は double 型です。
4 bairs( ) ベアストウ法 bairs(a, n, eps, iter, x); 全ての根を Complex 配列 x に入れて返します。
それ以外のパラメータは上と同じです。
5 regula( ) 超越
代数方程式		 レギュラ・
ファルシ法		 x = regula(xs,
xe, h, eps, iter);		 関数 _f(x) ＝ 0 となる実根を１つ求めます。
xs, xe は探索区間の下限と上限，
h はきざみ幅、eps は収束判定値、
iter は繰返し回数許容値です。
6 cnewton( ) 複素係数
高次
代数方程式		 ニュートン法 cnewton(a, x,
n, eps, iter);		 a は係数配列、n は次数、eps は収束判定値、
iter は繰返し回数許容値、x は根を入れる
配列で全ての根を求めます。
a, x は Complex 型です。
7 dka( ) ＤＫＡ法 dka(a, x, n, eps, iter);
プログラム・ソース（"poly.c"） 　　　　　　　　　　top （トップに戻る） 
/*		poly.c		*/
#include <stdio.h>
#include <stdlib.h>
#include "sslib.h"

void qurt(Complex a, Complex b, Complex c, Complex x[])
{
	Complex d, d1, d2, w;

	if(a.r == 0. && a.i == 0.)
	{
		fprintf(stderr, "Warning : a = 0  in qurt()\n");
		if(b.r == 0. && b.i == 0.)
		{
			fprintf(stderr, "Error : a = b = 0  in qurt()\n");
			x[0] = x[1] = tocomplex(0., 0.);
			return;
		}
		x[0] = cdiv(c, b);
		x[1] = tocomplex(0., 0.);
		return;
	}
	d = cmul1(a, -4.);
	d = cmul(d, c);
	w = cmul(b, b);
	d = cadd(d, w);
	if(d.r == 0. && d.i == 0.)
	{
		w = tocomplex(-b.r, -b.i);
		d = cmul1(a, 2.);
		x[0] = x[1] = cdiv(w, d);
		return;
	}
	d = csqrt(d);
	d1 = d2 = tocomplex(-b.r, -b.i);
	d1 = cadd(d1, d);
	d2 = csub(d2, d);
	if(cabslt(d2) > cabslt(d1))
	{
		w = d1;
		d1 = d2;
		d2 = w;
	}
	w = cmul1(a, 2.);
	x[0] = cdiv(d1, w);
	w = cmul1(c, 2.);
	x[1] = cdiv(w, d1);
	return;
}

void carda(double a[], Complex x[])
{
	int i;
	double alfa, beta, b1, b2, b3, d, p, q, theta, u3, v3, w;

	if(a[0] == 0.)
	{
		fprintf(stderr, "Error : a[0] = 0  in carda()\n");
		return;
	}
	b1 = *(a + 1) / *a / 3.;
	b2 = *(a + 2)/ *a;
	b3 = *(a + 3)/ *a;
	p = b2 / 3. - b1 * b1;
	q = -(b1 * (2. * b1 * b1 - b2) + b3);
	d = q * q + 4. * p * p * p;
	w = sqrt(fabs(d));
	if(d >= 0.)
	{
		u3 = (q + w) / 2.;
		v3 = (q - w) / 2.;
		alfa = cbrt(u3);
		beta = cbrt(v3);
		(*x).r = alfa + beta - b1;
		(*x).i = 0.;
		(*(x + 1)).r = - (alfa + beta) / 2. - b1;
		(*(x + 1)).i = 0.86602540378443864 * (alfa - beta);
		*(x + 2) = conj(*(x + 1));
		return;
	}
	theta = M_PI_2;
	if(q != 0.)	theta = atan(w / q);
	w = 2. * sqrt(-p);
	(*x).r = w * cos(theta / 3.) - b1;
	(*(x + 1)).r = - (w * cos((M_PI - theta) / 3.)) - b1;
	(*(x + 2)).r = - (w * cos((M_PI + theta) / 3.)) - b1;
	(*x).i = (*(x + 1)).i = (*(x + 2)).i = 0.;
	return;
}

double newton(double a[], int n, double eps, int iter)
/* １つの実根だけ求める	*/
{
	int i, ic, m;
	double a1, d, dx, p, q, r, w, wx;

	a1 = fabs(*a);
	if(a1 == 0. || n < 2 || eps <= 0. || iter < 1)
	{
		fprintf(stderr, "Error : Illegal parameter in newton()\n");
		return 0;
	}
	for(i = 1, wx = *a; i <= n; i++)
	{
		w = *(a + i);
		wx += w;
		if (a1 < fabs(w))	a1 = fabs(w);
	}
	wx = - wx / (double)(n + 1);
	for(ic = 1; ic <= iter; ic++)
	{
		p = *a;
		for(i = 1; i <= n; i++)	p = p * wx + *(a + i);
		m = n;
		q = *a * (double)m;
		for(i = 1; i < n; i++)
		{
			m--;
			q = q * wx + (double)m * *(a + i);
		}
		if(q != 0.)	dx = - p / q;
		else
		{
			m = n;
			r = *a * (double)(m * (n - 1));
			for(i = 1; i < n - 1; i++)
			{
				m--;
				r = r * wx + (double)(m * (m - 1)) * *(a + i);
			}
			d = p * (p - 2. * r);
			if(d <= 0.)	dx = - q / p;
			else
			{
				w = sqrt(d);
				if(q >= 0.)	dx = 2. * p / (- q - w);
				else		dx = 2. * p / (- q + w);
			}
		}
		wx += dx;
		if(fabs(dx) / a1 <= eps)	return wx;
	}
	fprintf(stderr, "Error : No convergence in newton()\n");
	return wx;
}

void bairs(double a[], int n, double eps, int iter, Complex x[])
{
	int flag, i, ic, j, k, nm, nm1, nn;
	double amax, bk, b1, b2, ck, c1, c2, c3, d, p, q, w, xx;

	for(i = 0; i < n; i++)	*(x + i) = tocomplex(0., 0.);
	if(*a == 0. || n < 2 || eps <= 0. || iter < 1)
	{
		fprintf(stderr, "Error : Illegal parameter  in bairs()\n");
		return;
	}
	nm = n;
	nm1 = n - 1;
	if((n % 2) != 0)
	{
		xx = newton(a, n, eps, iter);
		*(x + nm1) = tocomplex(xx, 0.);
		w = 0.;
		for(i = 0; i <= nm; i++)
		{
			w = w * xx + *(a + i);
			*(a + i) = w;
		}
		nm = nm1;
		nm1--;
	}
	amax = fabs(*a);
	for(i = 1; i <= nm; i++)
		if (amax < fabs(*(a + i)))	amax = fabs(*(a + i));
	i = nm1;
	while(i >= 3)
	{
		p = q = 1.;
		flag = 0;
		for(ic = 0; ic < iter; ic++)
		{
			b1 = b2 = c1 = c2 = 0.;
			nn = i + 1;
			for(k = 0; k <= nn; k++)
			{
				bk = *(a + k) - (p * b1 + q * b2);
				ck = bk       - (p * c1 + q * c2);
				if(k == (nn - 3))	c3 = ck;
				if(k != nn)
				{
					b2 = b1;
					b1 = bk;
					c2 = c1;
					c1 = ck;
				}
			}
			c1 = - (p * c2 + q * c3);
			d = c2 * c2 - c1 * c3;
			if(d != 0.)
			{
				p += (b1 * c2 - bk * c3) / d;
				q += (bk * c2 - b1 * c1) / d;
			}
			else
			{
				p++;
				q++;
			}
			if((fabs(b1) + fabs(bk + p * b1)) / amax <= eps)
			{
				flag = 1;
				break;
			}
		}
		if(!flag)
		{
			fprintf(stderr, "Error : No convergence in bairs()\n");
			return;
		}
		b1 = b2 = 0.;
		j = i - 1;
		for(k = 0; k <= j; k++)
		{
			bk = *(a + k) - (p * b1 + q * b2);
			*(a + k) = bk;
			b2 = b1;
			b1 = bk;
		}
		*(a + i) = p;
		*(a + i + 1) = q;
		i -= 2;
	}
	*(a + 1) = *(a + 1) / *a;
	*(a + 2) = *(a + 2) / *a;
	for(i = 1; i <= nm1; i += 2)
	{
		d = *(a + i) * *(a + i)	- 4. * *(a + i + 1);
		w = sqrt(fabs(d));
		if(d >= 0.)
		{
			*(x + i - 1) = tocomplex(- (*(a + i) - w) / 2., 0.);
			*(x + i)     = tocomplex(- (*(a + i) + w) / 2., 0.);
		}
		else
		{
			*(x + i - 1) = tocomplex(- *(a + i) / 2., w / 2.);
			*(x + i)     = conj(*(x + i - 1));
		}
	}
	return;
}

int cnewton(Complex a[], Complex r[], int n, double eps, int iter)
{
	Complex *d, w, x, p, q;
	int root, nn, i, ic;
	double ww, amax, am;

	d = (Complex *)malloc((n + 1)*sizeof(Complex));
	if(d == NULL)
	{
		fprintf(stderr, "Error : Out of Memory in cnwtn()\n");
		return(997);
	}

	if(ceq(a[0], tocomplex(0.0, 0.0)) || n < 2 || eps <= 0.0 || iter < 1)
	{
		fprintf(stderr, "Error : Illigal parameter in cnwtn()\n");
		free((char *)d);
		return(999);
	}

	root = 0;
	nn = n;

	w = a[0];
	if(cne(w, tocomplex(1.0, 0.0)))
	{
		r[0] = tocomplex(1.0, 0.0);
		for(i = 1; i <= n; i++)	r[i] = cdiv(a[i],w);
				/* normalizing coefficients */
	}
	else	for(i = 0; i <= n; i++)	r[i] = a[i];

	while(ceq(r[nn], tocomplex(0.0, 0.0)))
	{
		root++;
		nn--;
	}

	while(nn > 1)
	{
		amax = 1.0;
		for(i = 1; i <= nn; i++)
		{
			ww = r[i].r * r[i].r + r[i].i * r[i].i;
			if (amax < ww) amax = ww;
		}

		for(i = 0, am = nn; i < nn; i++, am--)
			d[i] = tocomplex(am * r[i].r, am * r[i].i);

		x = tocomplex(0.0, -1.0);

		for(ic = 0; ic < iter; ic++)
		{
			p = cfunc(r, nn, x);
			if(sqrt((p.r * p.r + p.i * p.i) / amax) <= eps)	break;
			q = cfunc(d, nn - 1, x); 
			w = cdiv(p, q);
			x = csub(x, w);
		}
		if(ic >= iter)
		{
			free((char *)d);
			if (nn == n)	return(-1);		/* no convergence , no root */
			return(root);					/* no convergence , any root */
		}
		r[nn] = x;
		nn--;
		root++;
		for(i = 1; i <= nn; i++)
		{
			w = cmul(r[i - 1], x);
			r[i] = cadd(r[i], w);
		}
	}
	w = cdiv(r[1], r[0]);
	r[1] = tocomplex(-w.r, -w.i);
	for(i = 1; i <= n; i++)	r[i - 1] = r[i];
	free((char *)d);
	return 0;
}

int dka(Complex a[], Complex r[], int n, double eps, int iter)
{
	Complex b, *c, f1, f2, *p, *q, *ratio, ri, z, zmax;
	double cf, di, dn, p0, p1, rmax;
	int flag, i, j, k;

	if(n < 1 || eps <= 0. || iter < 1)
	{
		fprintf(stderr, "Error : Illegal parameter  in dka()\n");
		return 999;
	}
	ratio = (Complex *)malloc((n + 1) * sizeof(Complex));
	if(ratio == NULL)
	{
		fprintf(stderr, "Error : Out of Memory  in dka()!!\n");
		return 997;
	}
	c = (Complex *)malloc((n + 1) * sizeof(Complex));
	if(c == NULL)
	{
		fprintf(stderr, "Error : Out of Memory  in dka()!!\n");
		free((char *)ratio);
		return 997;
	}

/* Normalize coefficients */

	z = *a;
	if(z.r != 1. || z.i != 0.)
			for(i = 0, p = a, q = c; i <= n; i++)	*q++ = cdiv(*p++, z);
	else	for(i = 0, p = a, q = c; i <= n; i++)	*q++ = *p++;

/* Get initial roots */

	dn = (double)n;
	rmax = 0.;
	b.r = (*(c + 1)).r / dn;
	b.i = (*(c + 1)).i / dn;
	for(i = 2, di = 2., p = c + 2; i <= n; i++, di++, p++)
	{
		cf = pow(cabslt(*p), 1.0 / di);
		if(rmax < cf)	rmax = cf;
	}

	p0 = 1.5 / dn;
	p1 = 2.0 * M_PI / dn;
	zmax = tocomplex(rmax, 0.);
	for(i = 1, di = 0.; i <= n; i++, di++)
	{
		z = cexp(tocomplex(0.0, p0 + di * p1));
		*(r + i) = csub(cmul(z, zmax), b);
	}

/* Main loop */

	k = 0;
	while(k++ < iter)
	{
		for(i = 1; i <= n; i++)
		{
			f1 = cfunc(c, n, (ri = *(r + i)));
			f2 = tocomplex(1.0, 0.0);
			for(j = 1; j <= n; j++)
				if (j != i)	f2 = cmul(f2, csub(ri, *(r + j)));
			*(ratio + i) = cdiv(f1, f2);
			*(r + i) = csub(ri, *(ratio + i));
		}

/* test convergence */

		flag = 1;
		for(i = 1; i <= n; i++)
		{
			if(cabslt(*(ratio + i)) > eps)
			{
				flag = 0;
				break;
			}
		}
		if(flag)
		{
			free((char *)c);
			free((char *)ratio);
			return 0;
		}
	}
	fprintf(stderr, "Error :  No Convergence in dka()\n");
	free((char *)c);
	free((char *)ratio);
	return -1;
}
プログラム・ソース（"regula.c"）　　　　　　　　　　　　　　　　　　　　　　　　　　　　top （トップに戻る） 
/*		regula.c		*/
#include <stdio.h>
#include "sslib.h"

double regula(double xs, double xe, double h, double eps, int iter)
{
	int ic;
	double x0, x1, x2, y0, y1, y2;

	if(xs == xe || h <= 0. || eps <= 0. || iter < 1)
	{
		fprintf(stderr, "Error : Illegal parameter in regula()\n");
		return 0.;
	}
	if(xs > xe)	swapd(&xs, &xe);
	x0 = xs;
	x1 = xe;
	y1 = _f(x1);
	do
	{
		y0 = _f(x0);
		if(y0 * y1 > 0.)	x0 += h;
		else
		{
			for(ic = 0; ic < iter; ic++)
			{
				x2 = x0 - y0 * (x1 - x0) / (y1 - y0);
				y2 = _f(x2);
				if(y0 * y2 >= 0.)
				{
					x0 = x2;
					y0 = y2;
				}
				else
				{
					x1 = x2;
					y1 = y2;
				}
				if(fabs(y2) <= eps)	return x2;
			}
			fprintf(stderr, "Error : No convergence in regula()\n");
			return x2;
		}
	} while(x0 <= xe);
	fprintf(stderr, "Error : No root in range(%e - %e) on regula()\n", xs, xe);
	return 0.;
}