#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. Функция решение кубического уравнения с действительными коэффициентами методом Виета-Кардано. Корни могут быть комплексными.
Кубическое уравнение записывается в виде:
x3+a*x2+b*x+c=0.
Для нахождения его корней, в случае действительных коэффициентов, вначале вычисляются:
Q = (a2-3b)/9
R = (2a3-9ab+27c)/54.
Далее, если R2 < Q3, то уравнение имеет три действительных корня, вычисляющихся по формулам (Виета):
е = acos(R/sqrt(Q3))/3,
x1 = -2*sqrt(Q)cos(t)-a/3,
x2 = -2*sqrt(Q)cos(t+(2*pi/3))-a/3,
x3 = -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.
Действительный корень будет:
x1 = (A+B)-a/3.
Комплексно-сопряженные корни:
x2,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. Положить Lii=1 i=0,...,N-1 (здесь ничего не производится, а только имеется в виду для дальнейших шагов;
2. Для каждого j=0,1,...,N-1 проделать:
1. Для всех i=0,...,j вычислить Uij:
Uij=Aij - SUM(0<=k<i)(Lik*Uik) (при i=0 сумму полагать нулем);
2. Для всех i=j+1,...,N-1 вычислить:
Lij = Aij - SUM(0<=k<j)(Lik*Uik))/Uii
Обе процедуры выполняются до перехода к новому 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)
{
register int i, imax, j, k;
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)
{
register int i, j, ip, ii = -1;
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)
{
register int i, j;
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)
{
register int j;
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. Для начального приближения x0, найти f0=f(x0), задать начальный интервал поиска D и его инкремент d>1.
2. Вычислить a=x0-D, b=x0+D; fa=f(a), fb=f(b).
3. Увеличить интервал поиска: D=D*d. Если интервал превысил некоторый заданный предел, то выйти с индикацией ошибки.
3.a. Если знаки fa и f0 отличаются, то считать корень окруженным на [a,x0] -> выход.
3.b. Если знаки fb и f0 отличаются, то считать корень окруженным на [x0,b] -> выход.
3.c. Если f0 > 0 (случай меньше нуля делается аналогично) алгоритм продолжается:
4. Проверяется, какое из fa или fb наименьшее. Если оба одинаковы, то переходим к 4a (двусторонний поиск), если fb - производим поиск вправо 4b, иначе - поиск влево 4c.
4.a. Находим a=a-D, b=b+D, fa=f(a), fb=f(b), идем к пункту 3.
4.b. Присваиваем последовательно a=x0, x0=b, fa=f0, f0=fb; находим b=b+D, fb=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), где все xi различны, существует единственный многочлен L(x) степени не более n, для которого L(xi) = yi.
В простейшем случае (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) представляет быстро сходящуюся (при хорошем начальном приближении) серию итераций: ai+1=0.5*(ai+x/ai), где 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
