:: алгоритмы  и методы :: :: олимпиадные задачи :: :: связь :: :: о сайте ::
Путь: Математика » Разное » Комплексная арифметика
На правах рекламы
Для вас со скидками интернет магазин детских кроватей машин по привлекательной цене.
  Комплексная арифметика



Кантор Илья по материалам
Numerical Recipes

Некоторые языки программирования обладают встроенной реализацией комплексных чисел, однако в наиболее распространенных типа Си, Бэйсик, и т.п. она напрочь отсутствует. Здесь мы рассмотрим, как качественно создать тип COMPLEX самому.

Cо сложением, как правило, вопросов не возникает. Стандартный путь

(a + ib)(c + id) = (ac - bd) + i(bc + ad)

вполне годится.

Не все так просто с модулем числа: запись |a + ib| = SQRT(a2 + b2) вызовет переполнение, если a или b близки к корню из самого большого возможного числа.

Правильный способ такой:

|a + ib| = |a|*SQRT(1+(b/a)2), если |a| >= |b|
|a + ib| = |b|*SQRT(1+(a/b)2), если |a| <   |b|

Деление:

division

Квадратный корень: сначала вычислить

root

Затем

Исходники даны в простейшем виде, в них отсутствует ряд проверок, например, на деление на ноль..


#include <math.h>

typedef struct FCOMPLEX {float r,i;} fcomplex;

fcomplex Cadd(fcomplex a, fcomplex b)
{
	fcomplex c;
	c.r=a.r+b.r;
	c.i=a.i+b.i;
	return c;
}

fcomplex Csub(fcomplex a, fcomplex b)
{
	fcomplex c;
	c.r=a.r-b.r;
	c.i=a.i-b.i;
	return c;
}
fcomplex Cmul(fcomplex a, fcomplex b)
{
	fcomplex c;
	c.r=a.r*b.r-a.i*b.i;
	c.i=a.i*b.r+a.r*b.i;
	return c;
}

fcomplex Complex(float re, float im)
{
	fcomplex c;
	c.r=re;
	c.i=im;
	return c;
}

fcomplex Conjg(fcomplex z)
{
	fcomplex c;
	c.r=z.r;
	c.i = -z.i;
	return c;
}

fcomplex Cdiv(fcomplex a, fcomplex b)
{
	fcomplex c;
	float r,den;
	if (fabs(b.r) >= fabs(b.i)) {
		r=b.i/b.r;
		den=b.r+r*b.i;
		c.r=(a.r+r*a.i)/den;
		c.i=(a.i-r*a.r)/den;
	} else {
		r=b.r/b.i;
		den=b.i+r*b.r;
		c.r=(a.r*r+a.i)/den;
		c.i=(a.i*r-a.r)/den;
	}
	return c;
}

float Cabs(fcomplex z)
{
	float x,y,ans,temp;
	x=fabs(z.r);
	y=fabs(z.i);
	if (x == 0.0)
		ans=y;
	else if (y == 0.0)
		ans=x;
	else if (x > y) {
		temp=y/x;
		ans=x*sqrt(1.0+temp*temp);
	} else {
		temp=x/y;
		ans=y*sqrt(1.0+temp*temp);
	}
	return ans;
}

fcomplex Csqrt(fcomplex z)
{
	fcomplex c;
	float x,y,w,r;
	if ((z.r == 0.0) && (z.i == 0.0)) {
		c.r=0.0;
		c.i=0.0;
		return c;
	} else {
		x=fabs(z.r);
		y=fabs(z.i);
		if (x >= y) {
			r=y/x;
			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
		} else {
			r=x/y;
			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
		}
		if (z.r >= 0.0) {
			c.r=w;
			c.i=z.i/(2.0*w);
		} else {
			c.i=(z.i >= 0) ? w : -w;
			c.r=z.i/(2.0*c.i);
		}
		return c;
	}
}

fcomplex RCmul(float x, fcomplex a)
{
	fcomplex c;
	c.r=x*a.r;
	c.i=x*a.i;
	return c;
}