Рекурсивный алгоритм решения задачи Ханойская башня

⇐ Назад

Далее ⇒

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include <string.h>

void move(int i, int j, int d)

{

int m;

for(m = 0; m < d; m++) printf(" ");

printf("перемещение(%d,%d)\n", i, j);

}

void hanoy(int i, int j, int k, int d)

{

// для отладки

//int m ;

//for (m = 0; m < d; m++) printf(" ");

//printf("hanoy(%d,%d,%d)\n", i, j, k);

if (k == 1)

move(i, j, d);

else

{

hanoy(i, 6-i-j, k-1, d+1);

hanoy(i, j, 1, d+1);

hanoy(6-i-j, j, k-1, d+1);

}

int main()

{

int n = 4;

printf("Введите количество дисков : ");

scanf("%d",&n);

printf("\nХаной %d дисками :\n", n);

hanoy(1, 3, n, 0);

getch();

// system ("pause");

}

Результат выполнения программы

Введите количество дисков : 4

Ханой с 4 дисками :

перемещение(1,2)

перемещение(1,3)

перемещение(2,3)

перемещение(1,2)

перемещение(3,1)

перемещение(3,2)

перемещение(1,2)

перемещение(1,3)

перемещение(2,3)

перемещение(2,1)

перемещение(3,1)

перемещение(2,3)

перемещение(1,2)

перемещение(1,3)

перемещение(2,3)

Пример 13. Функция решение кубического уравнения с действительными коэффициентами методом Виета-Кардано. Корни могут быть комплексными.

Кубическое уравнение записывается в виде:

x³+a*x²+b*x+c=0.

Для нахождения его корней, в случае действительных коэффициентов, вначале вычисляются:

Q = (a2-3b)/9

R = (2a3-9ab+27c)/54.

Далее, если R2 < Q3, то уравнение имеет три действительных корня, вычисляющихся по формулам (Виета):

е = acos(R/sqrt(Q3))/3,

x₁ = -2*sqrt(Q)cos(t)-a/3,

x₂ = -2*sqrt(Q)cos(t+(2*pi/3))-a/3,

x₃ = -2*sqrt(Q)cos(t-(2*pi/3))-a/3.

В том случае, когда R2 >= Q3, то действительных корней один (общий случай) или два (вырожденные случаи). Кроме действительного корня, имеется два комплексно-сопряженных. Для их нахождения вычисляются (формула Кардано):

A = -sign(R)[|R|+sqrt(R2-Q3)]1/3,

B = Q/A при A != 0 или B = 0 при A = 0.

Действительный корень будет:

x₁ = (A+B)-a/3.

Комплексно-сопряженные корни:

x_2,3 = -(A+B)/2-a/3 + i*sqrt(3)*(A-B)/2

В том случае, когда A = B, то комплексно-сопряженные корни вырождаются в действительный:

x2=-A-a/3.

Формулы Кардано и Виета требуют применения специальных функций, и в том случае, когда требуется провести большую серию вычислений корней кубического уравнения с не слишком сильно меняющимися коэффициентами, более быстрым алгоритмом является использование метода Ньютона или других итерационных методов (с нахождением начального приближения по формулам Кардано-Виета).

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include <string.h>

#define M_PI (3.141592653589793)

#define M_2PI (2.*M_PI)

int Cubic(double *x,double a,double b,double c);

Параметры:

a, b, c - коэффициенты

x - массив решения (size 3).

На выходе:

3 действительных корня -> затем x ими заполняется;

1 действительный + 2 комплексных -> x[0] - действительный, x[1] действительная часть комплексных корней, x[2] - неотрицательная мнимая часть.

Результаты:

1 - 1 действительный + 2 комплексных;

2 - 1 действительный корень + мнимая часть комплексных корней, если 0 (т.е. 2 действительных корня).

3 - 3 действительных корня;

int Cubic (double *x, double a, double b, double c)

{

double q, r, r2, q3;

q = (a * a - 3. * b) / 9.;

r = (a * (2. * a * a - 9. * b) + 27.*c) / 54.;

r2 = r * r;

q3 = q * q * q;

if (r2 < q3)

{

double t = acos(r / sqrt (q3));

a /= 3.;

q = -2. * sqrt (q);

x[0] = q * cos (t / 3.) - a;

x[1] = q * cos ((t + M_2PI) / 3.) - a;

x[2] = q * cos ((t - M_2PI) / 3.) - a;

return (3);

}

else

{

double aa, bb;

if (r <= 0.) r =- r;

aa = -pow (r + sqrt (r2 - q3), 1. / 3.);

if(aa != 0.) bb = q / aa;

else bb = 0.;

a /= 3.;

q = aa + bb;

r = aa - bb;

x[0] = q - a;

x[1] = (-0.5) * q - a;

x[2] = (sqrt (3.) * 0.5) * fabs (r);

if (x[2] == 0.) return (2);

return (1);

}

Пример 14. Алгоритм Краута (нижне-верхняя (LU) декомпозиция матрицы).

Алгоритм Краута - это фактически метод решения систем общего вида, конкурирующий по быстродействию с общеизвестным методом Гаусса-Жордана, но позволяющий более эффективно использовать решение.

Если мы можем разложить матрицу линейной системы A в произведение A=L*U, где L - нижняя, а U - верхняя треугольные матрицы, то решение системы уравнений с произвольной правой частью производится весьма просто, применением двух обратных подстановок. Более того, в отличие от известного метода Гаусса-Жордана, разложенная матрица позволяет быстро решать серии линейных уравнений с различными правыми частями при одной и той же матрице.

Метод Краута позволяет провести LU-декомпозицию матрицы примерно за то же число операций, что и "прямая" часть метода Гаусса-Жордана. Итоговые коэффициенты двух треугольных матриц упаковываются в матрицу того же размера, что и A, и на том же месте в памяти; при этом верхняя матрица U размещается в наддиагональной части и на диагонали; нижняя L в поддиагональной части, а диагональные элементы L считаются все равными 1 (без ограничения общности) и не выводятся.

Алгоритм Краута представляет из себя следующее:

1. Положить L_ii=1 i=0,...,N-1 (здесь ничего не производится, а только имеется в виду для дальнейших шагов;

2. Для каждого j=0,1,...,N-1 проделать:

1. Для всех i=0,...,j вычислить Uij:

U_ij=A_ij - SUM(0<=k<i)(L_ik*U_ik) (при i=0 сумму полагать нулем);

2. Для всех i=j+1,...,N-1 вычислить:

L_ij = A_ij - SUM(0<=k<j)(L_ik*U_ik))/U_ii

Обе процедуры выполняются до перехода к новому j.

Устойчивая работа алгоритма Краута возможно только при условии выбора ведущего элемента для каждого обращения к j из п.2 алгоритма. Выбор производится перед выполнением п.2 для очередного j перестановки необработанных строк матрицы так, чтобы строка, подставляемая на место j, содержала наибольший по модулю элемент в колонке j. Порядок перестановки необходимо запомнить, для чего достаточно использовать дополнительный вектор целых величин. Еще один вектор используется внутри алгоритма для масштабирования.

Алгоритм Краута имеет несколько приложений, одно из которых - решение системы линейных уравнений (обратная подстановка с учетом порядка перестановки строк), другие - вычисление обратной матрицы и вычисление детерминанта. Подробнее вызовы функций, связанных с алгоритмом Краута и реализация деталей алгоритма находятся в комментариях к программе и в ней самой.

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include <string.h>

#define TINY 1.e-30

LU-декомпозиция в соответствии с алгоритмом Краута.

Описание:

int LU_decompos(double **a,int n,int *indx,int *d,double *vv);

Параметры:

a - исходная матрица (n x n) на входе;

n - размер матрицы;

indx - массив целых чисел (размера n) для запоминания перемещений;

d - на выходе, содержит +1 или -1 для четного или нечетного числа перемещений.

vv - временный массив (size n).

Результат:

0 - некорректная входная матрица (непригодна для декомпозиции),

1 - положительный результат.

Обратная замена, использующая LU-декомпозицию матрицы.

Описание:

void LU_backsub(double **a,int n,int *indx,double *b);

Параметры:

a - матрица, разложенная по Крауту;

n - размер матрицы;

indx - порядок перемещения, полученный алгоритмом декомпозиции;

b - вектор (размера n).

Примечание: a и indx не изменяются в этой функции и могут быть использованы повторно.

Инвертирование матрицы с использованием LU-разложенной матрицы.

Описание:

void LU_invert(double **a,int n,int *indx,double **inv,double *col);

Параметры:

a - матрица, разложенная по Крауту;

n - размер матрицы;

indx - порядок перестановки, полученный алгоритмом декомпозиции;

inv - матрица назначения;

col - временный массив (размера n).

Получение матрицы, полученной с использованием LU-разложенной матрицы

Описание:

double LU_determ(double **a,int n,int *indx,int *d);

Параметры:

a - матрица, разложенная по Крауту;

n - размер матрицы;

indx - порядок перемещения, полученный алгоритмом декомпозиции;

d - знак равенства (+1 или -1) полученный при декомпозиции.

Результат: определяющее значение.

/* декомпозиция */

int LU_decompos (double **a, int n, int *indx, int *d, double *vv)

{

double big, sum, temp;

for (i = 0; i < n; i++)

{

big = 0.;

for (j = 0; j < n; j++)

if ((temp = fabs (a[i][j])) > big) big = temp;

if (big == 0.) return (0);

vv[i] = big;

}

/* главный цикл алгоритма Краута */

for (j = 0; j < n; j++)

{ /* это часть a) алгоритма, исключая i==j */

for (i = 0; i < j; i++)

{

sum = a[i][j];

for (k = 0; k < i; k++) sum -= a[i][k] * a[k][j];

a[i][j] = sum;

}

/* инициализация для поиска наибольшего элемента */

big = 0.;

imax = j;

for (i = j; i < n; i++)

{

sum = a[i][j];

for (k = 0; k < j; k++) sum -= a[i][k] * a[k][j];

a[i][j] = sum;

if ((temp = vv[i] * fabs (sum)) >= big)

{

big = temp;

imax = i;

}

/* обмен строк, если нужен, смена равенства и размера множителя */

if (imax != j)

{

for (k = 0; k < n; k++)

{

temp = a[imax][k];

a[imax][k] = a[j][k];

a[j][k] = temp;

}

*d =- (*d);

vv[imax] = vv[j];

}

/* сохраняем индекс */

indx[j] = imax;

if (a[j][j] == 0.) a[j][j] = TINY;

if (j < n - 1)

{

temp = 1. / a[j][j];

for (i = j+1; i < n; i++) a[i][j] *= temp;

}

return (1);

}

/* обратная замена */

void LU_backsub (double **a, int n, int *indx, double *b)

{

double sum;

/* первый шаг обратной замены */

for (i = 0; i < n; i++)

{

ip = indx[i];

sum = b[ip];

b[ip] = b[i];

if (ii >= 0)

for (j = ii; j < i; j++) sum -= a[i][j] * b[j];

else if (sum) ii = i; /* получен ненулевой элемент */

b[i] = sum;

}

/* второй шаг */

for (i = n - 1; i >= 0; i--)

{

sum = b[i];

for (j = i + 1; j < n; j++) sum -= a[i][j] * b[j];

b[i] = sum / a[i][i];

}

/* инвертирование матрицы */

void LU_invert (double **a, int n, int *indx, double **inv, double *col)

{

for (j = 0; j < n; j++)

{

for (i = 0; i < n; i++) col[i] = 0.;

col[j] = 1.;

LU_backsub (a, n, indx, col);

for (i = 0; i < n; i++) inv[i][j] = col[i];

}

/* определяющее вычисление*/

double LU_determ (double **a, int n, int *indx, int *d)

{

double res = (double)(*d);

for (j = 0; j < n; j++) res *= a[j][j];

return (res);

}

Пример 15. Алгоритм поиска интервала нахождения корня нелинейной функции.

Выяснение интервала, на котором корень содержится является важной проблемой поиска корня нелинейной функции действительной переменной. Здесь приведен алгоритм поиска такого интервала и ограничения на его применение. Примем, что корень функции f(x) окружен на интервале [a,b], если f(a) и f(b) имеют противоположные знаки. Чтобы окруженный согласно этому определению корень действительно существовал на этом интервале, достаточно непрерывности f(x), а для его единственности - еще и монотонности. При невыполнении этих свойств возможно отсутствие корня на [a,b] или неопределенность его позиции.

При использовании компьютера, всегда имеем дело с дискретным набором возможных представлений чисел (хотя и достаточно плотным). Кроме того, монотонность вычисленной функции может быть слегка нарушена в пределах точности ее вычисления. Это в ряде случаев усложняет вычисление окруженных корней функции, если к их точности предъявляются завышенные требования.

Окружение корня функции при гарантии ее определения на неограниченном интервале, производится по следующему итерационному алгоритму.

1. Для начального приближения x₀, найти f₀=f(x₀), задать начальный интервал поиска D и его инкремент d>1.

2. Вычислить a=x₀-D, b=x₀+D; f_a=f(a), f_b=f(b).

3. Увеличить интервал поиска: D=D*d. Если интервал превысил некоторый заданный предел, то выйти с индикацией ошибки.

3.a. Если знаки f_a и f₀ отличаются, то считать корень окруженным на [a,x₀] -> выход.

3.b. Если знаки f_b и f₀ отличаются, то считать корень окруженным на [x₀,b] -> выход.

3.c. Если f₀> 0 (случай меньше нуля делается аналогично) алгоритм продолжается:

4. Проверяется, какое из f_a или f_b наименьшее. Если оба одинаковы, то переходим к 4a (двусторонний поиск), если f_b - производим поиск вправо 4b, иначе - поиск влево 4c.

4.a. Находим a=a-D, b=b+D, f_a=f(a), f_b=f(b), идем к пункту 3.

4.b. Присваиваем последовательно a=x₀, x₀=b, f_a=f₀, f₀=f_b; находим b=b+D, f_b=f(b), идем к пункту 3.

4.c. Аналогично 4b, только направление поиска - влево.

Так как интервал поиска постоянно расширяется, то в конце концов используя указанный алгоритм корень будет окружен. Возможны модификации алгоритма в двух направлениях:

1. Увеличивать интервал не в геометрической прогрессии, а в арифметической либо по заданному сценарию;

2. Если область определения функции заведомо ограничена, то расширение интервала поиска также следует ограничивать имеющимися пределами, либо доопределять функцию там, где ее оригинал не определен.

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include <string.h>

Предполагается, что функция существует на бесконечном и непрерывном промежутке.

int BracketRoot (double x0, double *a, double *b, double d0, double di, double dmax, double (*fun)(double));

Параметры:

x0 - начальное приближение;

a - левая граница;

b - правая граница;

d0 - начальный интервал поиска;

di - интервал расширения (в геометрической прогрессии);

dmax - максимальный интервал;

fun - указатель на функцию.

Возвращает:

1 - если корень окружен;

0 - сбой

int BracketRoot (double x0, double *a, double *b, double d0,

double di, double dmax, double (*fun)(double))

{

double fa, fb, f0;

f0 = (*fun)(x0);

*a = x0 - d0;

*b = x0 + d0;

fa = (*fun)(*a);

fb = (*fun)(*b);

while ((d0 *= di) < dmax)

{

if (f0 >= 0.)

{

if (fa < 0.)

{

*b = x0;

return (1);

}

if (fb < 0.)

{

*a = x0;

return (1);

}

if (fa > fb)

{

*a = x0;

x0 = (*b);

*b += d0;

fa = f0;

f0 = fb;

fb = (*fun)(*b);

}

else

if (fa < fb)

{

*b = x0;

x0 = (*a);

*a -= d0;

fb = f0;

f0 = fa;

fa = (*fun)(*a);

}

else

{

*a -= d0;

*b += d0;

fa = (*fun)(*a);

fb = (*fun)(*b);

}

else

if (f0 < 0.)

{

if (fa >= 0.)

{

*b = x0;

return (1);

}

else

if (fb >= 0.)

{

*a = x0;

return (1);

}

if (fa < fb)

{

*a = x0;

x0 = (*b);

*b += d0;

fa = f0;

f0 = fb;

fb = (*fun)(*b);

}

else

if (fa > fb)

{

*b = x0;

x0 = (*a);

*a -= d0;

fb = f0;

f0 = fa;

fa = (*fun)(*a);

}

else

{

*a -= d0;

*b += d0;

fa = (*fun)(*a);

fb = (*fun)(*b);

}

return (0);

}

Пример 16. Алгоритм нахождения Полинома Лагранжа (интерполяция Лагранжа).

Интерполяционный многочлеен Лагранжа - многочлен минимальной степени, принимающий данные значения в данном наборе точек. Для n + 1 пар чисел (x_0, y_0), (x_1, y_1) \ dots (x_n, y_n), где все x_i различны, существует единственный многочлен L(x) степени не более n, для которого L(x_i) = y_i.

В простейшем случае (n = 1) - это линейный многочлен, график которого - прямая, проходящая через две заданные точки.

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include <string.h>

// таблица из учебника по вычислительной математике

float x[6] = {1.5, 1.54, 1.56, 1.60, 1.63, 1.70};

float y[6] = {3.873, 3.924, 3.950, 4.00, 4.037, 4.123};

/* Функция, вычисляющая коэффициенты Лагранжа

x - аргумент функции

n - степень многочлена (или число x-ов)

i - номер узла

float L(float xp, int n, int i)

{

// числитель и знаменатель

float Chesl;

float Znam;

Chesl = 1;

Znam = 1;

int k;

// вычисление числителя

for (k = 0; k != n; k++ )

{

if (k == i) continue;

// убираем множитель x - x(i)

Chesl *= xp - x[k];

}

// вычисление знаменателя

for (k= 0; k!= n; k++)

{

if (x[i] == x[k]) continue;

// убираем, а то ноль в знаменателе

Znam *= x[i] - x[k];

}

return Chesl / Znam;

}

int main()

{

// вычисляем степень полинома

int n = sizeof(y) / sizeof(float);

// начальное значение

float R = 0;

// произвольная точка для проверки

float px = 1.55;

// вычисляем значение интерполяционного многочлена в точке должно получиться 3. 937075

for (int i = 0; i != n; i++)

{

R += y[i] * L(px,n,i);

}

printf("Результат : %f \n",R);

system ("pause");

}

Результат выполнения программы

Результат : 3.937075

Пример 17. Алгоритм вычисление квадратного корня по алгоритму Ньютона.

Для вычисления квадратного корня в этом примере использован метод Ньютона. Рассмотрены машинно-зависимый и машинно-независимый варианты. Перед применением алгоритма Ньютона область определения исходного числа сужается до [0.5,2] (во втором варианте до [1,16]). Второй вариант машинно-независим, но работает дольше.

Вообще, это не самый быстрый вариант, но один из самых быстрых. Основная проблема заключается в том, что нет универсального машинного представления чисел с плавающей точкой, поэтому разделение числа на мантиссу и двоичную экспоненту как составляющих его компьютерного представления нельзя записать единым образом. Имено поэтому подключено описание математической библиотеки, из которой, впрочем, используются только frexp() и ldexp(). Конкретная реализация этих функций очень проста, но машинно-зависима. При разделении числа на мантиссу и экспоненту, мантисса оказывается в пределах [0.5,1). Если при этом экспонента нечетная, мантисса умножается на 2, тем самым область определения становится [0.5,2). К мантиссе применяется алгоритм Ньютона с начальным приближением 1. Окончательный результат получается с помощью применения ldexp() с половинным значением экспоненты.

Для машинно-независимого варианта выделение экспоненты заменяется серией последовательных делений исходного значения на 16, пока аргумент не станет определяться на интервале [1,16]. Если исходный аргумент был меньше 1, то перед серией делений обращаем его. Алгоритм Ньютона применяется с начальным приближением 2. Окончательный результат получается серией последовательных умножений на 4 и дальнейшим обращением, если аргумент обращался.

Сам алгоритм Ньютона для вычисления a = Sqroot(x) представляет быстро сходящуюся (при хорошем начальном приближении) серию итераций: a_i+1=0.5*(a_i+x/a_i), где i - номер итерации.

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include <string.h>

/* для одинарной точности нужны 4 итерации */

#define ITNUM 4

Разложение квадратного корня аргумента в мантиссу и экспонент (быстрый алгоритм):

float Sqroot (float x);

Вычисление квадратного корня без использования внешних фукнций (медленнее, но машинно-независим):

float Sqroot1 (float x);

В случае некорректного промежутка (x<0.) фукнции возвращают 0 и не генерят ошибку.

/* вариант с использованием внешних фукнций разложения/объединения на [0.5,1] */

float Sqroot (float x)

{

int expo, i;

float a, b;

if (x <= 0.F) return (0.F);

/* разложение x в мантиссу на промежутке [0.5,1) и экспонент.

Машинно-зависимые операции представлены вызовами функций. */

x = frexp (x, &expo);

/* нечетный экспонент: умножаем мантиссу на 2 и уменьшаем экспонент, делая его четным.

Теперь мантисса в промежутке [0.5,2.) */

if (expo & 1)

{

x *= 2.F;

expo--;

}

/* начальное приближение */

a = 1.F;

for (i = ITNUM; i > 0; i--)

{

b = x / a;

a += b;

a *= 0.5F;

}

/* делим экспонент на 2 и объединяем результат.

Фукнция ldexp() противоположна frexp. */

a = ldexp(a, expo / 2);

return (a);

}

/* Вариант без использования библиотек. Промежуток уменьшен до [1,16].

Используется 16 повторяющихся делений. */

float Sqroot1 (float x)

{

int sp = 0, i, inv = 0;

float a, b;

if (x <= 0.F) return (0.F);

/* аргумент меньше 1 : инвертируем его */

if (x < 1.F)

{

x = 1.F / x;

inv = 1;

}

/* последовательно делим на 16 пока аргумент не станет <16 */

while (x > 16.F)

{

sp++;

x /= 16.F;

}

/* начальное приближение */

a = 2.F;

/* Алгоритм Ньютона */

for (i = ITNUM; i > 0; i--)

{

b = x / a;

a += b;

a *= 0.5F;

}

while (sp > 0)

{

sp--;

a *= 4.F;

}

/* инвертируем результат для инвертированнго аргумента */

if (inv) a = 1.F / a;

return (a);

}

int main()

{

float x;

printf ("Введите число : ");

scanf ("%f", &x);

printf ("\nРезультат с использованием стандартных библиотек : %f\n", Sqroot(x));

printf ("Результат без использования стандартных библиотек : %f\n", Sqroot1(x));

system ("pause");

}

Результат выполнения программы

Введите число : 35

Результат с использованием стандартных библиотек : 5.916080

Результат без использования стандартных библиотек : 5.916080

⇐ Назад

Далее ⇒