その他関数ライブラリ



[ 簡単な説明 ]

特に分類できない小さなルーチンを集めました。

関数名 関数の機能 呼び出し例
normal( ) 乱数 正規乱数 x = normal(av, st);
rnd( ) 一様乱数 x = rnd( );
poison( ) ポアソン乱数 k = poison(av);
bino( ) 2項係数 k = bino(m, n);
filt( ) 逆ラプラス変換 filt(ts, te, f, nt, n, m, a1, b1, a, k1, k2, p);
kaijo1( ) 階乗 n!の計算 k = kaijo1(n);
kaijo2( ) 10log_kai2(n) による計算
(近似値)
k = kaijo2(n);
log_kai1( ) 階乗の常用対数 Σlog10( i ) による計算 a = kog_kai1(n);
log_kai2( ) 近似式使用 a = log_kai2(n);
minmax_i( ) 最大値、最小値 int 型配列 minmax_i(a, n, &min, &max);
minmax_d( ) double 型配列 minmax_d(a, n, &min, &max);
minmax_c( ) unsigned char 型配列 minmax_c(a, n, &min, &max);

プログラム・ソース("etc.c")           top (トップに戻る)
/*		etc.c		*/
#include <stdio.h>
#include <stdlib.h>
#include <values.h>
#include "sslib.h"

#define		uchar		unsigned char

int strcomp(uchar *a, uchar *b);

double normal(double av, double st)
/* av:平均値 , st:標準偏差 */
{
	float r1, r2, z1, z2;

	r1 = rnd();
	r2 = rnd();
	z1 = sqrt(-2. * log(r1));
	z2 = sin(6.283185307179586 * r2);
	return st * z1 * z2 + av;
}

double rnd(void)	/*  0  <=  rnd()  <  1.0  */
{
	return (double)rand() / RAND_MAX;
}

int poison(double av)
/* av:平均値	*/
{
	double c, x, w;
	int rn;

	if(av <= 0.)
	{
		fprintf(stderr, "Error : av <= 0.  in poison()\n");
		return 0.;
	}
	if(av >= 10.)
	{
		fprintf(stderr, "Warning : av >= 10.  in poison(). Use normal Random number!\n");
		return 0.;
	}
	c = exp(-av);
	x = 1.;
	rn = -1;
	while(x > c)
	{
		w = rnd();
		x *= w;
		rn++;
	}
	return rn;
}

double bino(int m, int n)
{
	double a[3], w1, w2;
	int k[3], i, j, l;

	if(m < n || n < 0)
	{
		fprintf(stderr, "Error : Illegal parameter  in bino()\n");
		return 0.;
	}
	if(m == n || n == 0)	return 1.;
	k[0] = m;
	k[1] = n;
	k[2] = m - n;
	for(i = 0; i < 3; i++)
	{
		w1 = 0.;
		l = k[i];
		for(j = 2; j <= l; j++)
		{
			w2 = l + 2 - j;
			w1 += log(w2);
		}
		a[i] = w1;
	}
	return exp(a[0] - a[1] - a[2]);
}

/*					逆ラプラス変換					*/

void filt(double ts, double te, double *f, int nt, int n, int m, double *a1, double *b1, double a, double k1, double k2, int p)
{
	double *am, *ap, expa, *ff, pit, *q, *r, sg, *t, td, w, x, y;
	double sr, si, z1r, z1i, z2r, z2i;
	int i, i1, j, k, l;

	if(ts >= te || nt < 1 || n < 0 || m < 0)
	{
		fprintf(stderr, "Error : Illegal parameter  in filt()\n");
		return;
	}
	t = (double *)malloc(nt * sizeof(double));
	if(t == NULL)
	{
		fprintf(stderr, "Error : Out of memory  in filt()\n");
		return;
	}
	ap = (double *)malloc((p + 1) * sizeof(double));
	if(ap == NULL)
	{
		fprintf(stderr, "Error : Out of memory  in filt()\n");
		free((char *)t);
		return;
	}

/* calc Euler coeffs */
	x = 1.;
	r = ap + p;
	*r-- = x;
	for(y = (double)p, w = 1.; y > 0; y--, w++, r--)
		*r = *(r + 1) + (x *= (y / w));

	expa = exp(a) / *ap;

	td = (te - ts) / (double)nt;
	for(i = 0, q = t; i < nt; i++)	*q++ = ts + td * (double)i;
	q = t;
	ff = f;
	if(*q == 0.)
	{
		*ff++ = 0.;
		i1 = 1;
		q++;
	}
	else	i1 = 0;
	for(i = i1; i < nt; i++)
	{
		k = k1 + *q * k2;
		pit = M_PI / *q;
		sr = a / *q;
		si = pit / 2.;
		sg = -1.;
		w = 0.;

		for(j = 1; j <= k; j++)
		{
			am = a1 + n - 1;
			z1r = *am--;
			z1i = 0.;
			for(l = n; l > 1; l--)
			{
				x   = sr * z1r - si * z1i;
				z1i = sr * z1i + si * z1r;
				z1r = x + *am--;
			}
			am = b1 + m - 1;
			z2r = *am--;
			z2i = 0.;
			for(l = m; l > 1; l--)
			{
				x   = sr * z2r - si * z2i;
				z2i = sr * z2i + si * z2r;
				z2r = x + *am--;
			}
			w += sg * (z1i * z2r - z1r * z2i) / (z2r * z2r + z2i * z2i);
			sg = - sg;
			si += pit;
		}

/* Euler translation */
		r = ap;
		w *= *r++;
		for(j = 1; j <= p; j++)
		{
			am = a1 + n - 1;
			z1r = *am--;
			z1i = 0.;
			for(l = n; l > 1; l--)
			{
				x   = sr * z1r - si * z1i;
				z1i = sr * z1i + si * z1r;
				z1r = x + *am--;
			}
			am = b1 + m - 1;
			z2r = *am--;
			z2i = 0.;
			for(l = m; l > 1; l--)
			{
				x   = sr * z2r - si * z2i;
				z2i = sr * z2i + si * z2r;
				z2r = x + *am--;
			}
			w += sg * (z1i * z2r - z1r * z2i) / (z2r * z2r + z2i * z2i) * *r++;
			sg = - sg;
			si += pit;
		}
		*ff++ = w * expa / *q++;
	}
	free((char *)t);
	free((char *)ap);
	return;
}

double kaijo1(int n)
{
	int i;
	double k;

	if(n < 0)
	{
		fprintf(stderr, "Error : n < 0  in kaijo1()\n");
		return 0.;
	}
	if(n > 170)
	{
		fprintf(stderr, "Error : Too large n (n > 170)  in kaijo1()\n");
		return MAXDOUBLE;
	}
	k = 1.;
	if(n <= 1)	return k;
	for(i = 2; i <= n; i++)	k *= (double)i;
	return k;
}

double kaijo2(int n)
{
	if(n < 0)
	{
		fprintf(stderr, "Error : n < 0  in kaijo2()\n");
		return 0.;
	}
	if(n > 170)
	{
		fprintf(stderr, "Error : Too large n (n > 170)  in kaijo2()\n");
		return MAXDOUBLE;
	}
	return pow(10., log_kai2(n));
}

double log_kai1(int n)
{
	int i;
	double k;

	if(n < 0)
	{
		fprintf(stderr, "Error : n < 0  in log_kai1()\n");
		return 0.;
	}
	k = 0.;
	if(n <= 1)	return k;
	for(i = 2; i <= n; i++)	k += log10((double)i);
	return k;
}

double log_kai2(int n)
{
	static double lp = 0.39908993417905751;
	double le, lf, lnn, nn;

	if(n < 0)
	{
		fprintf(stderr, "Error : n < 0  in log_kai2()\n");
		return 0.;
	}
	nn = (double)n;
	le = nn * M_LOG10E;
	lf = log10((((-571. / 2488320. / nn - 139. / 51840.) / nn + 1. / 288.) / nn
				 + 1. / 12.) / nn + 1.);
	lnn = (nn + 0.5) * log10(nn);
	return lnn - le + lp + lf;
}

void minmax_i(int a[], int n, int *min, int *max)
{
	int i;

	if(n == 1)
	{
		*min = *max = a[0];
		return;
	}
	if(a[0] == a[1])	*min = *max = a[0];
	else if(a[0] < a[1])
	{
		*min = a[0];
		*max = a[1];
	}
	else
	{
		*min = a[1];
		*max = a[0];
	}
	for(i = n - 1; i > 2; i -= 2)
	{
		if(a[i] >= a[i - 1])
		{
			if(*max < a[i])				*max = a[i];
			else if(*min > a[i - 1])	*min = a[i - 1];
		}
		else
		{
			if(*max < a[i - 1])		*max = a[i - 1];
			else if(*min > a[i])	*min = a[i];
		}
	}
	if(n % 2)
		if(*max < a[2])			*max = a[2];
		else if(*min > a[2])	*min = a[2];
}

void minmax_d(double a[], int n, double *min, double *max)
{
	int i;

	if(n == 1)
	{
		*min = *max = a[0];
		return;
	}
	if(a[0] == a[1])	*min = *max = a[0];
	else if(a[0] < a[1])
	{
		*min = a[0];
		*max = a[1];
	}
	else
	{
		*min = a[1];
		*max = a[0];
	}
	for(i = n - 1; i > 2; i -= 2)
	{
		if(a[i] >= a[i - 1])
		{
			if(*max < a[i])				*max = a[i];
			else if(*min > a[i - 1])	*min = a[i - 1];
		}
		else
		{
			if(*max < a[i - 1])		*max = a[i - 1];
			else if(*min > a[i])	*min = a[i];
		}
	}
	if(n % 2)
		if(*max < a[2])			*max = a[2];
		else if(*min > a[2])	*min = a[2];
}

void minmax_c(uchar *a[], int n, uchar **min, uchar **max)
{
	int i, c;

	if(n == 1)
	{
		*min = *max = a[0];
		return;
	}
	c = strcomp(a[0], a[1]);
	if(c == 0)	*min = *max = a[0];
	else if(c < 0)
	{
		*min = a[0];
		*max = a[1];
	}
	else
	{
		*min = a[1];
		*max = a[0];
	}
	for(i = n - 1; i > 2; i -= 2)
	{
		if(strcomp(a[i], a[i - 1]) >= 0)
		{
			if(strcomp(*max, a[i]))				*max = a[i];
			else if(strcomp(*min, a[i - 1]))	*min = a[i - 1];
		}
		else
		{
			if(strcomp(*max, a[i - 1]))		*max = a[i - 1];
			else if(strcomp(*min, a[i]))	*min = a[i];
		}
	}
	if(n % 2)
		if(strcomp(*max, a[2]))			*max = a[2];
		else if(strcomp(*min, a[2]))	*min = a[2];
}

int strcomp(uchar *a, uchar *b)
{
	uchar *p, *q;

	p = a;
	q = b;
	while(*p)
	{
		if(*p > *q)	return 1;
		if(*p++ < *q++)	return -1;
	}
	if(*q)	return -1;
	return 0;
}