Projekty z programování na FIT VUT: Projekt 2 — iterační výpočty

V předmětu Základy programování vyučovaném v prvním ročníku FITu. Druhý projekt se zaměřuje na základy iteračních výpočtů — tedy kdy se každým krokem zvyšuje přesnost daného čísla. Úkolem bylo znovu-objevit ameriku: výpočet tanh(x) a loga(x) — tangens hyperbolický a logaritmus o libovolném základu — a vážený aritmetický průměr wam a vážený kvadratický průměr — wqm. Součástí projektu byla nutnost vypracovat dokumentace — na konci tohoto článku je přiložené i PDF s dokumentací.

Dobrá rada ze začátku zdarma: NEOPISUJTE!!! Studujete-li FIT na VUT v Brně a řešíte projekt, neopisujte můj kód. Ačkoli je program licencován GPLv3+, tak nedoporučuji opisovat. Prvně — já již projekt odevzdal a pokud vím, tak ten FIŤácký soft pro odhalování plagiátů má ten „můj“ vzorek již uložen. A dále — pokud sami nebudete programovat, nenaučíte se to.

Zadání projektu je poměrně jednoduché — načti ze standardního vstupu hodnotu a vypiš odpovídající hodnotu logaritmu/tangentu. Opět zde nebudu uvádět detaily main() a ani načítání a kontrola parametrů. Mimo jiné jsem naprogramoval matematické funkce pow(), factorial() a exp().

izppow(double x, double y)

Jednoduchá funkce, která počítá xy. Hodnota x je libovolné desetinné číslo; y je libovolné celé číslo.

double izppow(double x, int y)
{
	double vysledek;

	if ( (x == 0.) && (y == 0.) )
	{
		vysledek = (NAN);
	}
	else
	{
		if (y > 0.)
		{
			vysledek = x * izppow(x,y-1.);
		}
		else if (y < 0.)
		{
			vysledek = 1 / izppow(x, -y);
		}
		else
		{
			vysledek = 1;
		}
	}
	return vysledek;
}

izpfactorial(int n)

Faktorial, neboli n! = n(n-1)(n-2)…(n-n+1) je matematická rekursivní funkce.

unsigned long izpfactoriall (int n)
{
	unsigned long vysledek;

	n = (unsigned long) n;

	if (n < 0)
		return -1;

	if ( n == 0)
	{
		return 1;
	}
	else
	{
		vysledek = n * izpfactoriall(n-1);
	}
	return vysledek;
}

izpexp(double x, int sigdig)

Funkce izpexp() počítá hodnotu ex na zadaný počet desetinných míst (sigdig). Funkce využívá dvou různých algoritmů. První vychází z taylorovského rozvoje funkce kolem nuly:

exp(x) = SUMn->inf ( xn / n! )

...
do{
	y = yold + izppow(x,n)/izpfactoriall( n );
	n++;
	delta = fabs( rad * ( y - yold ) );
	if (n > MAXFACTCNT)
		return -1.;

	yold = y;
} while ( delta > 1. );

	return y;

Druhý využívá limitní definice exp():

limn->inf = ( 1 + x/n )n

Tedy pěkný, ale přesto není nejrychlejší.

...
double n = 10.;
do
{
	y = izppow((1+x/n), (int) n);
	n = n + 10.;
	delta = fabs( rad * ( y - yold ) );
	yold = y;
} while ( delta > 1. );
return yold;

Tento algoritmus je ovšem stabilní — nepotřebuje výpočet faktoriálu — a je proto použit v intervalu (-inf; -1)  a (+1; +inf).

Funkce tanh

Konečně se dostáváme k funkci, kterou jsme měli za úkol naprogramovat. Tanh je možné definovat buďto rozvojem a nebo

tanh(x) = (e2x – 1) / (e2x + 1)

Funkce pro výpočet tanh(x) jen volá funkci e2x a zvyšuje případně požadovanou přesnost.

Funkce logax

Pro výpočet logaritmu je možné využít několik přístupů, já jsem ovšem zvolil ten, který je pro počítače nejpřirozenější — výpočet log2(x) — tedy binární logaritmus. Počítá velice rychle a princip výpočtu je například ve knížce od Knutha — Umění programování 1.

Celý algoritmus je rozdělen na 2 části — výpočet celé části logaritmu a výpočet desetinné části. Výpočet celé části je poměrně jednoduchý — pokud je možné dělit dané číslo 216k, poté je celá část rovna 2k. Poté se číslo vydělí a kontroluje se 28k atd.

Celá část tedy bude

int izplgxbin(int x)
{
	if (x == 0)
		return -1;

	int bin = 0;
	if (x >= 65536)
	{
		x /=65536;
		bin += 16;
	}
	if (x >= 256)
	{
		x /= 256;
		bin +=  8;
	}
	if (x >= 16)
	{
		x /= 16;
		bin +=  4;
	}
	if (x >= 4)
	{
		x /= 4;
		bin +=  2;
	}
	if (x >= 2)
		bin +=  1;
	return bin;
}

Poměrně jednoduché, ne? Teď ta desetinná část — její princip je částečně popsán v dokumentaci: pokud je hodnota v rozmezí (1–2), tak ji násobíme samu sebou tak dlouho, dokud není větší jako dva a zmenšujeme dělíme dvěma zvanou zlomek. Poté se opět zavolá tatáž funkce s parametry x/2 (tedy opět dostaneme x do intervalu (1–2) ), parametr zlomek a nakonec nějakým způsobem předáváme hodnotu přesnosti. V mém pojetí byla přesnost počet iterací — při každé iteraci se snížila.

double izplgxfrac(double x,double zlomek,int presnost)
{
	if (presnost == 0)
		return 0;
	if (x == 1.)
		return 0.;
	presnost--;
    while (x < 2.)
	{
		zlomek /= 2.;
     		x *= x;
	}
	return zlomek + izplgxfrac(x/2,zlomek,presnost);
}

Statistické funkce

Nakonec byly požadovány statistické funkce. Vzhledem k tomu, že se jednalo o opravdu jednoduché záležitosti, nebudu zde nijak moc rozebírat. Pouze uvedu, že výpočet probíhal v jednom cyklu a pomocí struktury STATISTICS se předávala data. Poslední dvě subroutiny tedy vypadaly

struct STATISTICS
{
	double x;
	double h;
	double citatel;
	double jmenovatel;
	double vysledek;
};

int izpwam(struct STATISTICS *st)
{
	st->citatel += (st->x*st->h);
	st->jmenovatel += st->h;
	st->vysledek = st->citatel / st->jmenovatel;
	return 0;
}

int izpwqm(struct STATISTICS *st)
{
	st->citatel += (st->x*st->x*st->h);
	st->jmenovatel += st->h;
	st->vysledek =  sqrt( (st->citatel)/(st->jmenovatel) );
	return 0;
}

Při načítání dat bylo samozřejmě ověřováno, jestli je váha větší (nebo rovna) nule, jestli je x a h opravdu číslo a podobně. Ovšem nebudu zde vypisovat seznam nejrůznějších podmínek, stejně jako jsem nevypisoval podmínky pro logaritmus (jaký v v nule, v nekonečnu a to v závislosti na základu a) a nebo zkratku v tanh(), kde funkce od jisté hodnoty vracela hodnotu asymptotickou — 1.

Celý program se zdrojovými kódy spolu s dokumentací je možné stáhnout, rozbalit a přeložit. Oproti minulé verzi jsem komentáři mnohem více šetřil, protože minule mi bylo vyčteno, že kód musí být samovysvětlující (souhlasím) a že takřka nejsou potřeba (nesouhlasím). A nebo si přečtěte samotnou dokumentaci (jednalo se o projekt patřící do jiného předmětu — Softwarové inženýrství). Pokud vás zajímají detaily k dokumentaci, brzy vydám článek rozebírající i toto. Takže to chce klid, nohy v teple.

V případě, že byste rádi použili některou část, dodržte prosím licenci GPL verze 3 a nebo pozdější.

Napsat komentář

Vaše emailová adresa nebude zveřejněna. Vyžadované informace jsou označeny *