Некоторые математические приложения
В этой главе рассмотрим различные математические штуки, которые можно использовать благодаря программированию. Некоторое мы уже рассматривали ранее, здесь же будут другие и более сложные вещи.
1. Вычисление производной в точке.
Предположим, у нас есть функция f(x) и нам нужно найти её производную в точке a, т.е. f'(a). С помощью программирования это делается легко и мгновенно:
double a; cin >> a;
double dx= 0.000001;
double pr= ( f(a+dx) - f(a) ) / dx;
cout << pr;
На примере функции x*x. Её производная в точке x будет равна 2*x. Это для проверки. А теперь программа:
#include <iostream>
using namespace std;
double f (double x) { return x*x; }
int main () {
double x; cin >> x; double dx= 0.000001;
double pr= ( f(x+dx) - f(x) ) / dx;
cout << pr << endl;
system ("pause");
}
Прикол в том, что вместо return x*x мы можем вписать любую, даже какую-нибудь очень сложную функцию, например: log(x*sin(x) - cos(pow(x,x*atan(x)))) * sqrt (x*cbrt(x+log2(2*x))); Для компьютера это как 2 пальца об асфальт, он всё равно найдёт ответ за долю секунды (если этот ответ существует). Вручную мы могли бы считать это примерно полчаса, если не больше.
2. Вторая производная в точке.
Теперь посмотрим как считать производную от производной. Например, для функции куба f(x)= x^3 будет производная 3*x^2, а вторая производная: 6*х. Это для проверки.
Сделаем функцию, которая считает производную в точке х:
double ff (double x, const double dx= 0.000001) {
return ( f(x+dx) - f(x) ) / dx;
}
И тогда сама вся программа:
|
|
|
#include <iostream>
#include <cmath> // для функции куба необязательно
using namespace std;
const double dx= 0.000001;
double f (double x) { return x*x*x; }
double ff (double x) { return ( f(x+dx) - f(x) ) / dx; }
int main () {
double x; cin >> x;
double pr= ( ff(x+dx) - ff(x) ) / dx;
cout << pr << endl;
system ("pause");
}
Тут уже заметна погрешность: для х=1 выводит 5.99987 вместо 6. Чтобы уменьшить погрешность, нам нужно уменьшить число dx. Проблема в том, что в вычислениях берётся только 6 знаков после запятой и в связи с этим число dx просто рубится до нуля, если сделать его меньше.
Можно записать программу чуть проще:
#include <iostream>
#include <cmath> // для функции куба необязательно
using namespace std;
const double dx= 0.000001;
double f (double x) {
return x*x*x;
}
double f1 (double x) {
return ( f(x+dx) - f(x) ) / dx;
}
double f2 (double x) {
return ( f1(x+dx) - f1(x) ) / dx;
}
int main () {
double x; cin >> x;
cout << f2(x) << endl;
system ("pause");
}
Если заменить приближение производной на ( f(x+dx) - f(x-dx) ) / 2*dx; то точность вычисления заметно возрастает.
Аналогичным продолжением можно посчитать третью производную и т.д. При каждой новой производной погрешность будет увеличиваться. Рекурсивный метод вычисления производной n-го порядка не найден.
|
|
|
3. Определённый интеграл или площадь криволинейной трапеции.
Здесь существует 3 стандартных численных метода: прямоугольников, трапеций и симпсона (парабол). Первый даёт больше погрешности, последний- сложный и запутанный в реализации. Поэтому будем использовать наиболее интуитивно очевидный- метод трапеций.
Берётся некий шаг dx, который как можно ближе к нулю. Считается интеграл от функции f(x) от A до B. Тогда промежуток равен (B-A) и число разбиений примерно будет равно N= (B-A)/dx. Можно задать число разбиений точно, тогда шаг dx вычислится автоматически: dx= (B-A)/N. Начинаем с точки A и идём считать N слагаемых- площади маленьких трапеций.
double sum= 0;
for (int i=0; i<N; i++) sum+= ( f(A + i*dx) + f(A + (i+1)*dx) ) * dx /2;
Тогда в переменной sum будет значение интеграла.
Полная программа на примере функции x^2. Первообразная будет (x^3)/3, а интеграл (B^3 - A^3) / 3. Можно свериться.
#include <iostream>
#include <cmath> // для f(x)=x*x необязательно
using namespace std;
double f (double x) { return x*x; }
int main () {
double A,B,N; cin>>A>>B>>N;
double dx= (B-A)/N, sum= 0;
for (int i=0; i<N; i++) sum+= ( f(A + i*dx) + f(A + (i+1)*dx) ) * dx /2;
cout << sum << endl;
system ("pause");
}
Вместо f(x) можно в return прописать любую интегрируемую функцию и ответ будет найден так же быстро. Чем больше разбиений (N), тем меньше погрешность, но тем дольше считает программа.
|
|
|
4. Вычисление длины графика.
Чтобы вычислить длину графика на некотором промежутке, надо разбить его на большое количество (N) маленьких ломаных и длины этих ломаных сложить. Пусть у нас есть функция f(x). Тогда мы получаем функцию кусочка ломаной из теоремы Пифагора: L(x)= sqrt (1 + f'(x)) и далее это проинтегрировать. Мы поступим немножечко по другому.
Чтобы получить длину кусочка ломаной, нужно знать 2 точки (концы куска). Пусть промежуток от a до b, разбиваем на N кусков. Тогда шаг dx= (b-a)/N; И первая точка куска ломаной будет a;f(a), вторая a+dx;f(a+dx). Расстояние между 2 точками запишем в функцию:
double q2 (double x) {
return x*x;
}
double dist (double x1, double x2) {
return sqrt ( q2(x2-x1) + q2(f(x2)-f(x1)) );
}
И далее само суммирование в цикле:
double sum=0;
for (int i=0; i<N; i++) sum+= dist (a+i*dx, a+(i+1)*dx);
И в значении sum получим приблизительную длину графика.
Пример полной программы, которая вычисляет длину параболы y=x*x на промежутке от A до B с разбиением на N кусков:
#include <iostream>
#include <cmath> // для f(x)=x*x необязательно
using namespace std;
double f (double x) { return x*x; } // сюда можно подставить любую непрерывную однозначную функцию
|
|
|
double q2 (double x) { return x*x; }
double dist (double x1, double x2) { return sqrt ( q2(x2-x1) + q2(f(x2)-f(x1)) ); }
int main () {
double A,B,N; cin>>A>>B>>N;
double dx= (B-A)/N, sum= 0;
for (int i=0; i<N; i++) sum+= dist (A+i*dx, A+(i+1)*dx);
cout << sum << endl;
system ("pause");
}
5. Вычисление пределов.
Тут совсем всё просто- мы должны просто взять значение, наиболее близкое к предельному. Например, нужно вычислить предел sin(x)/x при x стремящемся к нулю. Как известно, это первый замечательный предел и он равен 1. Но мы можем посчитать его (и вообще любой предел) на компьютере: подставим х=0.000001. Получится 1. В данном случае чем ближе значение x к нулю, тем точнее будет вычислен предел. Можно даже вносить в отдельную функцию, например:
double f (double x) { return pow(2,x)/(x*log2(x)); }
И посчитаем предел при x стремящемся к бесконечности. Для этого надо подставить очень большое число, например x= 100000;
Полная программа:
#include <iostream>
#include <cmath>
using namespace std;
double f (double x) { return pow(2,x)/(x*log2(x)); }
int main() {
cout << f(100000) << endl;
system("pause");
}
Компилятор на devcpp (TDM-GCC 4.9.2) вывел inf, что значит "бесконечность", от слова infinity. Иногда ещё выводится nan (not a number)- неопределённость. В visual studio выводится 1.#INF что тоже значит бесконечность. Ты можешь потестировать с этим: попробуй сначала вывести 1/0, а потом 0/0.
Если x стремится к минус бесконечности, то нужно подставить что-то вроде f(-100000);
Если x стремится к 3, то мы можем подставлять как слева (2.9999), так и справа (3.0001). Например, для предела (x*x-9)/(x-3) при х=3 будет неопределённость 0/0. Имеем f(2.9999)= 5.9999, f(3.0001)= 6.0001. Нетрудно догадаться, что предел равен 6. Ручной способ тут тоже прост: нужно лишь разложить в числителе разность квадратов и сократить на знаменатель. Проблема в том, что не всегда попадаются такие лёгкие и понятные пределы- вот тогда-то на помощь может и прийти компьютер.
6. Метод полного перебора.
Предположим, нам нужно решить квадратное уравнение. Мы уже писали программу для его решения через дискриминант. А кубическое? Ну тут уже нужны новые формулы. Получается, для каждого типа уравнений нужны свои формулы. Но ладно бы это- 99% уравнений вообще не решается через формулы, т.е. аналитически. Поэтому приходится использовать приблизительные методы, которые подходят вообще для любого типа уравнений. Самым простым и интуитивно понятным методом является метод полного перебора.
Для того чтобы решить уравнение методом полного перебора, нужно примерно понимать, в каком промежутке у нас находится корень (решение) и оно должно быть там одно. И тогда мы разбиваем этот промежуток на большое количество кусочков и подставляем значения. Уравнение записываем в виде f(x)= 0. И запоминаем x, при котором было наиболее близкое к нулю значение. Это и будет корнем уравнения.
Посмотрим на примере sin(x) + x = 1. Если записать уравнение в виде sin(x)= 1-x и нарисовать графики, понятно что x находится где-то от 0 до 1. Разобьём на миллион кусков и будем подставлять:
double f (double x) { return sin(x)+x-1; }
double a= 0, b= 1, N= 1000000; double dx= (b-a)/N;
double min= abs(f(a)), root= a;
for (int i=1; i<=N; i++) {
if ( abs(f(a+i*dx)) < min ) {
min= abs(f(a+i*dx));
root= a+i*dx;
}
}
cout << "x= " << root << endl;
Пример программы, которая ищет корень уравнения (x^x)*ln(x)= 10 на заданном промежутке:
#include <iostream>
#include <cmath>
using namespace std;
double f (double x) {
return pow(x,x)*log(x)-10;
}
int main() {
double a,b,r;
while(true) {
system("CLS");
cout << "a= "; cin >> a;
cout << "b= "; cin >> b;
cout << "Accuracy: "; cin >> r;
double N= pow(10,r);
double dx= (b-a)/N;
double min= abs(f(a)), root= a;
for (int i=1; i<=N; i++) {
if ( abs(f(a+i*dx)) < min ) {
min= abs(f(a+i*dx));
root= a+i*dx;
}
}
cout << "x= " << root << endl;
system("pause");
}
}
При a=0, b=10, r=6 программа нашла корень за секунду: х= 2.54218. Если мы получаем на выходе либо a, либо b, то тут либо одна из крайних точек является корнем уравнения (что обычно вряд ли), либо на этом промежутке просто нет корней.
7. Метод монте-карло.
Суть метода можно понять на примере вычисления площади, хотя у него есть и другие применения. Например, посмотрим вычисление площади круга или даже нахождение числа пи через этого: мы будем ставить очень много случайных точек. Потом количество точек, попавших в круг, поделим на площадь квадрата, в который мы заключаем этот круг. И мы получим долю площади, которую занимает круг в квадрате. Если сделать квадрат 1 на 1, то получится площадь самого круга с радиусом 0.5. Если затем мы поделим площадь круга на квадрат радиуса, то получим число пи. Область круга мы можем учитывать через радиус или через уравнение графика. Будем считать площадь 1 четвертинки круга.
Создание случайного числа от 0 до 0.5: double x= rand()%500 / 1000; Точка будет складываться из 2 случайных чисел- координата по х и по y. Пусть у нас будет миллион таких точек:
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;
int main () {
srand(time(NULL));
double inCircle= 0, N=1000000, R= 0.5; double x,y,r;
for (int i=0; i<N; i++) {
x= rand()%500; y= rand()%500;
x/=1000; y/=1000;
r= sqrt(x*x+y*y);
if (r<= R) inCircle++;
}
cout << inCircle << endl;
double S= inCircle/N; double pi= S/(R*R); cout << pi << endl;
system("pause");
}
Погрешность значительна даже при миллионе точек, поэтому площадь находить так не очень эффективно.
Вот более универсальная программа, где можно подставлять уже другую функцию. Правда, и промежутки тоже нужно менять:
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;
double f (double x) {
return sqrt (0.25-x*x);
}
int main () {
srand(time(NULL));
double inCircle= 0, N=1000000; double x,y;
for (int i=0; i<N; i++) {
x= rand()%500; x/=1000;
y= rand()%500; y/=1000;
if (y <= f(x)) inCircle++;
}
cout << inCircle << endl;
double S= inCircle/N; double pi= S*4; cout << pi << endl;
system("pause");
}
Можно этот метод применить и к решению уравнения: подбирать точки случайно, пока не будет меньше некой погрешности. Для уравнения sin(x)+x-1=0 программа:
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;
double f (double x) { return sin(x)+x-1; }
int main () {
srand(time(NULL));
double N= 1000000; double e= 0.0001; double x;
do { x= rand()%1000; x/=1000; } while ( abs(f(x)) > e);
cout << "x= " << x << endl;
system("pause");
}
8. Алгоритм быстрой рекурсивной сортировки.
Данный алгоритм позволяет сортировать массивы в 10 раз большим объёмом, чем мы можем сортировать обычной нерекурсивной сортировкой, рассмотренной ранее. Т.е. эта сортировка более быстрая и эффективная, но алгоритм достаточно громоздкий и тяжёлый в понимании. Лично я потратил полчаса на то, чтобы разобрать что вообще тут происходит- и позже подредактировал и переписал в своём стиле- что предлагаю сделать и тебе. Тут суть в том, что массив делится на 2 куска определённым образом, а для сортировки каждого из кусков используется эта же самая функция. И таким многократным делением массив сортируется.
template <typename Type>
void qSort(Type m[], int size) {
Type terek = 0;
int left= 0, right= 0;
if(size>1) {
Type *Left = new Type[size];
Type *Right = new Type[size];
terek = m[size/2];
for(int i=0;i<size;i++) {
if (i!=size/2) {
if (m[i]<terek) {
Left[left] = m[i]; left++;
}
else {
Right[right] = m[i]; right++;
}
}
}
qSort(Left,left);
qSort(Right,right);
for(int I=0;I<size;I++){
if (I<left) {
m[I]= Left[I];
}
else if (I==left) {
m[I]= terek;
}
else {
m[I]= Right[I-(left+1)];
}
}
}
}
Массив из 100 тыс целочисленных элементов отсортировался на моём компьютере за 0.28 сек. Миллион сортировать отказался- не тянет. 300 тыс элементов- 1.8 сек. 500 тыс- 5.8 сек.
Для сравнения возьмём нерекурсивную сортировку: 100 тыс- 52 сек (почти минута), 10 тыс- 0.52 сек. Т.е. с ростом массива в 10 раз время сортировки растёт в 100 раз. Или по-другому для оценки скорости и эффективности алгоритма пишут: О(N^2). Для быстрой сортировки О(N*log(N)).
9. Нахождение простых чисел решетом Эратосфена.
Ниже будут примеры функций, которые находят простые числа. Простое число- это число, не имеющее делителей (кроме себя самого и 1). Суть алгоритма решета Эратосфена в том, что мы берём какую-то область от 1 до некого числа, где ищем простые числа. Например от 1 до 1000000. Сначала зачёркиваем числа, кратные 2. Затем кратные 3, 4, 5 и т.д. Оставшиеся числа будут простыми.
// вывод простых чисел в промежутке [1;max]
void _SimpleIn (int max) {
bool *z = new bool [max+1];
for(int I=0;I<=max;I++) z[I]=true;
for(int J=2;J<=max;J++) {
if(z[J]) cout << J << " ";
for(int I=J;I<=max;I+=J) z[I]=false;
}
}
// количество простых чисел в промежутке [1;max]
int SimpleIn (int max) {
bool *z = new bool [max+1];
for(int I=0;I<=max;I++) z[I]=true;
int n=0;
for(int J=2;J<=max;J++) {
if(z[J]) n++;
for(int I=J;I<=max;I+=J) z[I]=false;
}
return n;
}
// х-ое простое число
int nSimple (int x) {
int max = x*30; bool *z = new bool [max+1];
for(int I=0;I<=max;I++) z[I]=true;
int n=0;
for(int J=2;J<=max;J++) {
if(z[J]) {
n++; if(n==x) return J;
}
for(int I=J;I<=max;I+=J) z[I]=false;
}
}
// вывести первые х простых чисел
void _nSimple (int x) {
int max = x*30; bool *z = new bool [max+1];
for(int I=0;I<=max;I++) z[I]=true;
int n=0;
for(int J=2;J<=max;J++) {
if(z[J]) {
cout << J << " "; n++;
}
if(n==x) break;
for(int I=J;I<=max;I+=J) z[I]=false;
}
}
10. Проценты.
Я тут даже писать ничего не буду, тут и так всё понятно.
// процент х от y
double proc (double x, double y) { return x*100/y; }
// вывод процента
void _proc (double x, double y) { cout << " (" << proc(x,y) << "%) "; }
Если мы работаем с типом int, то тогда мы сначала умножаем на 100 обязательно, а потом уже делим на что нужно. И получим целочисленный процент. Например, процент 24 от 35: 24*100= 2400, 2400/35= 68 (дробная часть обрежется)
11. Функции из комбинаторики.
Здесь представлены некоторые часто встречающиеся.
// факториал
long long int fact (int x) {
if(x==0) return 1;
else return x*fact(x-1);
}
// n- внизу, k- наверху (число сочетаний)
int cnk (int n, int k) {
return fact(n)/(fact(k)*fact(n-k));
}
// двойной факториал
long long int fact (int x, int k) {
if(k==1) return fact(x);
if(k==2) {
if(x%2==0) return fact(x/2)*pow(2,x/2);
if(x%2==1) return fact(x)/fact(x-1,2);
}
}
// вывести треугольник паскаля
void _pascaltring (int x) {
for(int j=0; j<x+1;j++) {
for(int i=0; i<j+1; i++) cout<<cnk(j,i)<<" ";
cout<<endl;
}
}
// рекурсивный алгоритм х-е число фибоначчи
int fib (int x) {
if((x==0)||(x==1)) return 1;
else return fib(x-1)+fib(x-2);
}
// быстрый нерекурсивный алгоритм для поиска чисел фибоначчи (значения до миллиона)
double fibfast (int x) {
const double f = (1+sqrt(5))/2;
return (pow(f,x) - pow(-f,-x)) / (2*f-1);
}
// числа каталана
int katalan (double n) {
if ( (n==0) || (n==1) ) return 1;
if (n>=2) {
int res= 0;
for (int i=0; i<n; i++) res+= katalan(i)*katalan(n-1-i);
return res;
}
}
12. Вычисление вероятности статистическим методом.
К примеру, мы кидаем 5 кубиков. Нам нужно найти вероятность того что сумма всех очков будет равна 19. Для этого мы можем с помощью компьютера симулировать бросок кубика (генерация псевдослучайных чисел) и сделать миллион бросков. Будем считать те броски, где сумма очков 19. Далее поделим количество этих бросков на миллион (т.е. на все броски вообще). Профит!
#include <iostream>
#include <cstdlib>
#include <ctime>
using namespace std;
int main() {
srand(time(NULL));
int N= 1000000; int x19= 0; int ball;
for (int i=1; i<=N; i++) {
ball= 0; for (int j=1; j<=5; j++) ball+=rand()%6 + 1;
if(ball==19) x19++;
}
cout << (double)x19*100/N << "%" << endl;
system("pause");
}
Программа показала примерную вероятность: 9.4746%. Новые запуски: 9.45%, 9.4945%, 9.457%. Т.е. результат колеблется возле 9.4-9.5%. Для большей точности мы можем сделать больше бросков, например 10 миллионов- надо просто заменить значение int N. Тут программа "думала" уже пару секунд. Имеем: 9.44472%. Можно сделать и 100 миллионов бросков- придётся немного подождать, но результат будет высчитываться чуть точнее. Получили: 9.45464%. Можно также усреднить все колебания: сосчитать вероятность 10 раз и вычислить среднее арифметическое. У меня при 1 млн бросках и 10 повторениях в среднем получилось 9.45626%. Колебания также могут быть связаны с самим алгоритмом генерации случайных чисел. В итоге можно сделать вывод, что вероятность примерно равна 9.45%.
Задание на тренировку: вспомни какую-нибудь свою интересующую тебя ситуацию с кубиками и попробуй посчитать вероятность статистическим методом. Иногда это может быть интересно, особенно если ты любитель подобных азартных игр. В случае азартных игр это даже может помочь принять какое-то решение.
Пример текущей ситуации: посмотрим при 5 кубиках какая сумма баллов наиболее вероятная. Возможно от 5 (1+1+1+1+1) до 30 (6+6+6+6+6) баллов. И посчитаем просто все вероятности и выведем:
#include <iostream>
#include <cstdlib>
#include <ctime>
using namespace std;
int main() {
srand(time(NULL));
for (int k=5; k<=30; k++) {
int N= 1000000; int x= 0; int ball;
for (int i=1; i<=N; i++) {
ball= 0; for (int j=1; j<=5; j++) ball+=rand()%6 + 1;
if(ball==k) x++;
}
cout << k << " --> " << (double)x*100/N << "%" << endl;
}
system("pause");
}
Результат:
5 --> 0.0113%
6 --> 0.0641%
7 --> 0.1883%
8 --> 0.4507%
9 --> 0.8964%
10 --> 1.6043%
11 --> 2.6427%
12 --> 3.9541%
13 --> 5.3723%
14 --> 6.9436%
15 --> 8.386%
16 --> 9.4824%
17 --> 10.0265%
18 --> 10.0391%
19 --> 9.4957%
20 --> 8.3676%
21 --> 6.8555%
22 --> 5.4059%
23 --> 3.9605%
24 --> 2.6356%
25 --> 1.6164%
26 --> 0.9096%
27 --> 0.4466%
28 --> 0.1952%
29 --> 0.0652%
30 --> 0.015%
13. Уравнение 3-й степени по формуле Кардано.
У нас есть уравнение вида A*(x^3) + В*(x^2) + C*x + D = 0. Нужно найти все действительные решения. Это кубическое уравнение (или кубур). Попробуем решить его аналитически в общем виде.
Сначала мы введём эти коэффициенты. Сделаем также функцию для вывода уравнение. Уравнение будем выводить в виде (A, B, C, D) * (x:3) = 0. По сути здесь скалярное произведение вектора коэффициентов на степенной базис.
void print (double a, double b, double c, double d, double x0) {
cout << "(" << a << " ," << b << " ," << c << " ," << d << ") * (";
if(x0==0) cout << "x:3)";
if(x0>0) cout << "x-" << x0 << ":3)";
if(x0<0) cout << "x+" << x0 << ":3)";
cout << " = 0" << endl;
}
cout << "Решатель уравнения вида: A*(x^3) + В*(x^2) + C*x + D = 0" << endl;
cout << "А: "; cin >> A; cout << "B: "; cin >> B;
cout << "C: "; cin >> C; cout << "D: "; cin >> D;
cout << "Исходное уравнение: "; print(A,B,C,D,0);
Если A=0, то уравнение не кубическое- и мы ничего не решаем. Ставим это условие.
Далее делим всё уравнение на старший член: a= B/A; b= C/A; c= D/A; print(1,a,b,c,0);
Далее делаем замену х= у-а/3; И после раскрытия скобок второй член сократится. Получим:
cout << "Замена х= y";
if(a>0) cout << "-" << abs(a/3);
if(a<0) cout << "+" << abs(a/3);
cout << endl;
p= b - a*a/3; q= 2*a*a*a/27 - a*b/3 + c;
print (1,0,p,q,-a/3);
Далее можно сделать замену y= t - p/(3*t), получить квадратное уравнение, а из него уже формулу Кардано. Для упрощения вводится дискриминант: s= q*q/4 + p*p*p/27. От его знака будет зависеть результат решения.
Если дискриминант больше 0, то уравнение имеет 1 решение, которое высчитывается по формуле Кардано. Оставшиеся 2 комплексные корня можно высчитать через схему Горнера и квур, если это нужно.
Если дискриминант равен 0, тогда 2 решения, которые тоже можно формульно задать. Одно решение найдём через формулу Кардано, затем поделим схемой Горнера и получим квадратное уравнение с 1 корнем через выделение полного квадрата.
Если дискриминат меньше 0, тогда у нас 3 корня, причём что интересно- они будут действительные. Просто мнимые части сократятся.
m= -q/2;
if (s>0) { n= sqrt(s); cout << "x= " << cbrt(m-n) + cbrt(m+n) - a/3 << endl; }
На этом первая часть решения нашей задачи закончена. Теперь нужно разобрать ещё случаи когда дискриминант равен 0 или меньше 0.
Если равен 0, то первый корень найти не составит труда: x1= 2*cbrt(m) - a/3; Далее по схеме Горнера мы поделим на y^3 + p*y + q на y-y1 и получим:
y^2 + 2*cbrt(m) + cbrt(m)^2 = (y+cbrt(m))^2 = 0, отсюда y2= -cbrt(m) или x2= -cbrt(m) - a/3;
Т.к. компьютер считает с погрешностью, то иногда дискриминант может вместо нуля получиться очень маленьким числом. Поэтому пропишем так:
if ( (s==0) || (abs(s)<pow(10,-15)) )
Если дискриминант меньше 0, придётся обратиться в теорию комплексных чисел и тригонометрическую запись комплексного числа.
Число вида a+b*i, где i= sqrt(-1) можно представить в тригонометрическом виде, а именно- вынести модуль sqrt(a*a+b*b) за скобки.
M * (cos(F) + i*sin(F)), где F- тоже некое число, которое называют аргументом.
Тогда кубический корень из этого комплексного числа будет равен:
cbrt(M) * ( cos((F+2*pi*k)/3) + i*sin((F+2*pi*k)/3) ), где k принимает 3 значения: 0, 1, 2.
И далее подставив это в формулу Кардано, получим: 2*cbrt(M)*cos((F+2*pi*k)/3);
Здесь M= sqrt(m*m+n*n), n= sqrt(-s).
С аргументом формула чуть сложнее: если q<0, то F= atan(n/m). Если q>0, то надо ещё прибавить "пи", если q=0, то F= pi/2.
Таким образом, вся программа:
#include <iostream>
#include <cmath>
using namespace std;
const double pi= 3.1415926535;
void print (double a, double b, double c, double d, double x0) {
cout << "(" << a << " ," << b << " ," << c << " ," << d << ") * (";
if(x0==0) cout << "x:3)";
if(x0>0) cout << "x-" << abs(x0) << ":3)";
if(x0<0) cout << "x+" << abs(x0) << ":3)";
cout << " = 0" << endl;
}
int main() {
setlocale(0,""); double A,B,C,D,a,b,c,p,q,s,m,n,x,x1,x2,_x[3],M,F;
while(true) {
system("CLS");
cout << "Решатель уравнений вида: A*(x^3) + Â*(x^2) + C*x + D = 0" << endl;
cout << "À: "; cin >> A;
cout << "B: "; cin >> B;
cout << "C: "; cin >> C;
cout << "D: "; cin >> D;
cout << "Исходное уравнение: "; print(A,B,C,D,0);
//
cout << "Деление на старший член: ";
a= B/A; b= C/A; c= D/A; print(1,a,b,c,0);
//
cout << "Замена х= y";
if(a>0) cout << "-" << abs(a/3);
if(a<0) cout << "+" << abs(a/3);
cout << endl;
p= b - a*a/3; q= 2*a*a*a/27 - a*b/3 + c;
print (1,0,p,q,-a/3);
//
cout << "Дискриминант: "; s= q*q/4 + p*p*p/27; cout << s << endl;
//
m= -q/2;
if (s>0) {
n= sqrt(s); x= cbrt(m-n) + cbrt(m+n) - a/3;
cout << "x= " << x << endl;
}
if ( (s==0) || (abs(s)<pow(10,-15)) ) {
x1= 2*cbrt(m) - a/3; x2= -cbrt(m) - a/3;
if (x1==x2) cout << "x= " << x1 << endl;
else cout << "x1= " << x1 << endl << "x2= " << x2 << endl; }
if (s<0) {
n= sqrt(-s); M= sqrt(m*m+abs(s));
if (q<0) F= atan(n/m);
if (q>0) F= atan(n/m)+ pi;
if (q==0) F= pi/2;
for (int i=0; i<=2; i++) {
_x[i]= 2*cbrt(M)*cos((F+2*pi*i)/3) - a/3;
cout << "x" << i+1 << "= " << _x[i] << endl;
}
}
cout << endl; system("pause");
}
}
14. Обратные неэлементарные функции.
Как мы можем сосчитать арифметический корень, если б у нас не было такой функции? Предположим, надо сосчитать корень из 7. Ну мы считаем 1*1=1. Если 1<7, то продолжаем. 2*2=4. 4<7? Продолжаем. 3*3=9. 9<7? Нет. Тогда теперь считаем 2.1*2.1= 4.41 и т.д. по такому принципу. Так же мы можем считать и другие обратные функции. Ниже программа на примере обратной функции от f(x)= x*(e^x), тут e- число "е" или 2.71828 примерно. Введём обратную функцию как kat(x), т.е. если f(x)= y, то x= kat(y). Например, kat(2.71828)= 1 (приблизительно). Тут считаются значения только больше нуля, т.к. только на этом промежутке функция монотонна и однозначна.
#include <iostream>
#include <cmath>
using namespace std;
const double e= exp(1);
double f (double x) { return x * pow (e,x) ; }
double kat (double x) {
if(x==0) return 0;
double y=0; double k=0;
while(k>(-7)) {
if((f(y)<x)&&(f(y+pow(10,k))>x)) k--;
else y+=pow(10,k);
}
return y;
}
int main() {
setlocale(0,""); double x;
while(true) {
system("CLS"); cout << "x= "; cin >> x;
if(x>=0) cout << "Подбор со степенным шагом: kat(x)= " << kat(x) << endl;
system("pause");
}
}
15. Синтезатор линейных уравнений с рациональным решением.
Привожу сразу саму программу:
#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cmath>
using namespace std;
void info (int win, int lose) {
system("CLS");
cout << "Решено: " << win << "/" << win+lose << endl << endl << endl;
}
template <typename Type>
Type solve (Type A[]) {
return (A[3] - A[1]) / (A[0] - A[2]);
}
double* syntesis () {
double *res= new double [4];
int *eq= new int [4];
do {
for(int i=0;i<4;i++) {
eq[i]= -9 + rand()%19;
res[i]= eq[i];
}
} while ( (res[0]==res[2]) || (1000*solve(res)!=1000*solve(eq)) );
return res;
}
void print (double A[]) {
cout << "Решите уравнение: " << endl;
if(A[0]!=0) {
if(abs(A[0])!=1) cout << A[0] << "*x";
if(A[0]==1) cout << "x";
if(A[0]==-1) cout << "-x";
}
if((A[1]>0)&&(A[0]!=0)) cout << " + ";
if(A[1]<0) {
if(A[0]!=0) cout << " ";
cout << "-";
if(A[0]!=0) cout << " ";
}
if(A[1]!=0) cout << abs(A[1]) << " ";
if((A[0]==0)&&(A[1]==0)) cout << 0 << " ";
if((A[1]==0)&&(A[0]!=0)) cout << " ";
cout << "= ";
if((A[2]==0)&&(A[3]==0)) cout << 0;
if(A[2]!=0) {
if(abs(A[2])!=1) cout << A[2] << "*x";
if(A[2]==1) cout << "x";
if(A[2]==-1) cout << "-x";
}
if((A[3]>0)&&(A[2]!=0)) cout << " + ";
if(A[3]<0) {
if(A[2]!=0) cout << " ";
cout << "-";
if(A[2]!=0) cout << " ";
}
if(A[3]!=0) cout << abs(A[3]);
cout << endl;
cout << "x= ";
}
int main() {
srand(time(NULL)); setlocale(0,"");
int win=0, lose=0; bool stop= false; int k=0;
int t1= clock();
while(!stop) {
info(win,lose);
double *A= syntesis(); print(A);
double x; cin >> x;
if(x==solve(A)) { cout << "Верно"; win++; }
else { cout << "Неверно"; lose++; }
cout << endl << endl;
system("pause");
info(win,lose);
cout << "(0) Продолжить\n(1) Закончить" << endl;
cin >> stop;
system("CLS");
}
info(win,lose);
cout << "Решено " << win*100/(win+lose) << "%" << endl;
double x= clock()-t1; x/=60000;
cout << "Скорость: " << win/x << " уравнений в минуту" << endl;
cout << endl;
system("pause");
}
16. Симулятор так иль так (ботва онлайн)
Суть игры в том, что игрок делает ставку валютой и выбирает тип ставки- больше, меньше или равно. Далее он бросает 3 кубика и считает сумму очков. Далее противник (бармен) бросает тоже 3 кубика и у него будет своя сумма очков. Если игрок ставил на больше/меньше и его сумма очков получилась больше/меньше чем у бармена, то ставка удваивается. Т.е. поставил 1000 и вернули 2000. Если же нет, то вся ставка сгорает. В случае победы игры на равно ставка увеличивается в 5 раз.
Программа:
#include <iostream>
#include <cstdlib>
#include <ctime>
using namespace std;
void _() { cout << endl; }
void __(int k) {
if(k==0) system("CLS");
if(k==1) {
_(); system("pause");
}
if(k==2) {
__(1); __(0);
}
}
void res (int a, int k[]) {
cout << k[0] << " " << k[1] << endl;
if(a==0) cout << "Ты проиграл" << endl;
if(a==1) cout << "Ты выиграл" << endl;
}
int main() {
setlocale(0,""); srand(time(NULL));
double stavka; int type;
cout << "Начальный капитал: "; double $; cin >> $;
bool passtype= ( (type==0) || (type==1) || (type==2) );
int kub, kub_en;
while ($>0) {
__(0);
cout << "Твои деньги: " << $ << endl; _();
do {
cout << "Ставка: "; cin >> stavka;
if(stavka<=0) cout << "Ставка не может быть меньше или равна нулю!\n";
if(stavka>$) cout << "Ставка не может быть больше чем у тебя есть денег!\n";
} while ( (stavka<=0)||(stavka>$) );
do {
cout << "(0) Равно\n(1) Меньше\n(2) Больше\n"; cin >> type;
if(!passtype) cout << "Ошибка! Неправильный выбор типа ставки!\n";
} while (!passtype);
int kub[2]= {0,0};
for (int i=0; i<=1; i++) for (int j=1; j<=3; j++) kub[i]+=rand()%6+1;
if ( (type==0)&&(kub[0]==kub[1]) ) {
$+=stavka*4; res(1,kub);
}
if ( ((type==0)&&(kub[0]!=kub[1])) || ((type==1)&&(kub[0]>=kub[1])) || ((type==2)&&(kub[0]<=kub[1])) ) {
$-=stavka; res(0,kub);
}
if ( ((type==1)&&(kub[0]<kub[1])) || ((type==2)&&(kub[0]>kub[1])) ) {
$+=stavka; res(1,kub);
}
__(1);
}
__(0); cout << "Игра закончена! Ты проиграл все деньги!\n"; __(1);
}
17. Ещё из НОД, НОК и т.п.
Мы уже рассмотрели функцию НОД (наибольший общий делитель), как нерекурсивную, так и рекурсивную. Наименьшее общее кратное находится просто: nok(a,b)= a*b/nod(a,b). Ниже приведены все функции.
int nod (int x, int y) { int z; while(y!=0) { z=x; x=y; y=z%x; } return x; }
int nok (int x, int y) { return x*y/nod(x,y); }
// количество итераций при нахождении НОД
int nnod (int x, int y) {
int k=0, z;
while(y!=0) {
z=x; x=y; y=z%x;
k++;
}
return k;
}
// проверка на взаимно простые
bool isInterSimple (int x, int y) {
if (nod(x,y)==1) return true;
else return false;
}
18. Определитель матрицы и системы линуров.
#include <iostream>
#include <cmath>
using namespace std;
// функция для генерации минора из матрицы
double* minor (int p, double A[], int n) {
double *M= new double [(n-1)*(n-1)]; int s=0;
for (int i=0; i<(n*n); i++) {
if((i%n==p%n)||(i/n==p/n)) continue;
M[s]=A[i]; s++;
}
return M;
}
// функция для высчитывания определителя
double det (double A[], int n) {
if(n==1) return A[0];
else {
double res=0;
for(int i=0;i<n;i++) {
res+= A[i] * det(minor(i,A,n),n-1) * pow(-1,i);
}
return res;
}
}
// главная функция
int main() {
setlocale(0,"");
while(true) {
system("CLS");
cout << "Размерность матрицы: "; int R; cin >> R;
double *a= new double [R*R];
// ввод
for (int i=0; i<R*R; i++ ) {
cout << "Parca " << i+1 << ": "; cin >> a[i];
}
cout << endl;
// вывод
for (int i=0; i<R*R; i++) {
cout << a[i] << " ";
if(i%R==R-1) cout << endl;
}
cout << endl;
// результат
cout << "det("<<R<<")= " << det(a,R) << endl;
cout<<endl; system("pause");
}
}
Теперь давайте подробно разберёмся что тут происходит. Рассмотрим на примере матрицы 3 порядка.
Предположим, у нас есть матрица с такими элементами:
(7) (-3) (5)
(2) (0) (8)
(-1) (6) (5)
В скобках- значения элементов матрицы. Как мы считаем определитель: берём элементы первой строки и умножаем на алгебраическое дополнение. Тут полезно бы обратиться к теории, как этот определитель вообще считается, и всё станет понятно что происходит в программе.
Т.е. определитель равен 7 * алгебр.доп ячейки 0 + (-3)* а.д.ячейки 1 + 5*а.д.ячейки 2. Мы всю матрицу загоним в одномерный массив и пронумеруем ячейки от 0 до 8. Алгебрическое дополнение- это минор, но с учётом знака. Знаки чередуются. В рекурсивной функции det ты можешь это видеть.
Теперь как сделать минор- он зависит от самой матрицы, её размерности и от позиции. Минор- это определитель новой матрицы, которая образуется путём зачёркивания элементов по той же вертикали и по той же горизонтали. Т.е. для ячейки 0 минором будет определитель матрицы:
(0) (8)
(6) (5)
Для ячейка 1 получится минорная матрица:
(2) (8)
(-1) (5)
И определитель такой матрицы считается по такому же принципу. Определителем матрицы 1 на 1 будет само значение элемента. Так мы получаем рекурсию.
Нужно лишь занести эти значения в минорную матрицу- надо пропустить ненужные значения. Для ячейки 1 мы пропускаем 0, 1, 2 по горизонтали и 1, 4, 7 по вертикали. По горизонтали- это целочисленное деление на 3, по вертикали- остаток при делении на 3. Теперь можешь взглянуть на функцию минора и сам уже в ней разобраться.
Теперь как решать системы линуров: объясню на простейшем примере 2х2:
7 *х1 + 5*х2 = -2
-3*х1 + 1*х2 = 0
Мы делаем матрицу системы:
(7) (5)
(-3) (1)
И далее матрицы путём замещения значениями справа. Для х1 будет матрица:
(-2) (5)
(0) (1)
И для х2 аналогично. Чтобы получить х1, нужно поделить определитель матрицы для х1 на определитель матрицы системы. Это называется методом Крамера. И его можно обобщить на системы любой размерности, где кол-во уравнений совпадает с кол-вом неизвестных.
Чтобы написать программу для решения системы линейных уравнений, нужно добавить новую функцию для создания замещённой матрицы и можно ещё сделать функцию для вычисления икса. Полная программа:
#include <iostream>
#include <cmath>
using namespace std;
// матрица из минора на позиции p
double* minor (int p, double A[], int n) {
double *M= new double [(n-1)*(n-1)];
int s=0;
for (int i=0; i<(n*n); i++) {
if((i%n==p%n)||(i/n==p/n)) continue;
M[s]=A[i]; s++;
}
return M;
}
// определитель
double det (double A[], int n) {
if(n==1) return A[0];
else {
double res=0;
for(int i=0;i<n;i++) {
res+= A[i] * det(minor(i,A,n),n-1) * pow(-1,i);
}
return res;
}
}
// матрица с замещением вектора-столбца на столбце p
double *zam (double A[], int n, int p, double B[]) {
double *M= new double [n*n];
for(int i=0; i<n*n; i++) {
if (i%n==p%n) M[i]= B[i/n];
else M[i]= A[i];
}
return M;
}
double X (double A[], int n, int p, double B[]) { return det(zam(A,n,p,B),n) / det(A,n); }
int main() {
setlocale(0,"");
while(true) {
system("CLS");
cout << "Размерность системы: "; int R; cin >> R;
double *a= new double [R*R]; double *b= new double [R];
// ввод матрицы системы
for (int i=0; i<R*R; i++) { cout << "Parca " << i+1 << ": "; cin >> a[i]; }
cout << endl;
// ввод вектора справа от знака равно
for (int i=0; i<R; i++) { cout << "Droit " << i+1 << ": "; cin >> b[i]; }
// вывод
for (int i=0; i<R*R; i++) {
cout << a[i] << " ";
if(i%R==R-1) cout << " = " << b[i/R] << endl;
} cout << endl;
// результат
for (int i=0; i<R; i++) { cout << "x" << i+1 << "= " << X(a,R,i,b) << endl; }
cout<<endl; system("pause");
}
}
19. Дифференциальные уравнения.
С помощью простого программирования можно решать даже дифференциальные уравнения первого порядка с заданным начальным условием. Общий вид решаемого уравнения: y' = f(x,y), где y(x0)= y0. Мы будем использовать самый простой метод- метод Эйлера.
Суть в том, что мы заменяем функцию на ломаную линию. И чем больше изломов, тем точнее результат. При этом куски ломаных будут параллельны соответствующим касательным в точках графика. Здесь мы находим именно координаты этих изломов, поэтому решением будет таблица чисел x, y.
На примере: dy/dx= x*x - 2*y, причём y(0)= 1. Решение ищем на отрезке [0;1].
Считаем f(x0,y0)= f(0,1)= -2
Вводим шаг dx= (B-A)/N, где [A;B]- заданный отрезок, N- число разбиений. Пусть dx= 0.1 для наглядности.
Тогда: dx*f(0,1)= -0.2
И далее:
x1= 0.1, y1= y0 + dx*f(x0,y0);
x2= 0.2, y2= y1 + dx*f(x1,y1);
Вносим всё в цикл:
_y= y0;
for (int i=1; i<=N; i++) {
x= A + dx*i; y= _y + dx*f(A+dx*(i-1), _y);
cout << "x= " << x << " --> y= " << y << endl;
_y= y;
}
А теперь вся программа:
#include <iostream>
#include <cmath> // для функции x*x-2*y необязательно
using namespace std;
double f (double x, double y) { return x*x-2*y; }
int main() {
double A, B, N; A=0, B=1, N=100; double dx= (B-A)/N;
double x0,y0; x0=0, y0=1;
double x,y, _y= y0;
for (int i=1; i<=N; i++) {
x= A + dx*i; y= _y + dx*f(A+dx*(i-1), _y);
cout << "x= " << x << " --> y= " << y << " --> f(x,y)= " << f(x,y) << endl;
_y= y;
}
cout << endl; system("pause");
}
Метод можно усовершенствовать- чтобы касательные были не на левых краях, а посерединке каждого кусочка ломаной.
Тогда изменится следующее:
double x,y, _y= y0, _x= x0;
for (int i=1; i<=N; i++) {
x= _x + dx; y= _y + dx*f(_x + dx/2, _y + (dx/2)*f(_x,_y));
cout << "x= " << x << " --> y= " << y << " --> f(x,y)= " << f(x,y) << endl;
_y= y; _x= x;
}
Ну и напоследок- более точный, но и более сложный метод: Рунге-Кутты.
double x,y, _y= y0, _x= x0; double k1,k2,k3,k4;
for (int i=1; i<=N; i++) {
x= _x + dx;
k1= f(_x,_y);
k2= f(_x+dx/2, _y + dx*k1/2);
k3= f(_x+dx/2, _y + dx*k2/2);
k4= f(_x+dx, _y+dx*k3);
y= _y + (dx/6)*(k1+2*k2+2*k3+k4);
cout << "x= " << x << " --> y= " << y << " --> f(x,y)= " << f(x,y) << endl;
_y= y; _x= x;
}
20. Бонус! Коллекция функций с массивами.
// v 1.3 [13.01.2017]
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <string>
using namespace std;
const double pi= 3.141592653;
const double kgrad= 180 / pi;
//////////////////////////////////////////////////////////////////////////////////////
// ПЕРЕРАБОТКА САМИХ МАССИВОВ //
// создание копии массива
template <typename Type>
Type* ArrayCopy (Type m[], int size) {
Type *M= new Type [size];
for(int i=0; i<size; i++) M[i]= m[i];
return M;
}
// частичная копия с n1 до n2 элемента
template <typename Type>
Type* ArrayCopy (Type m[], int size, int n1, int n2) {
Type *M= new Type [n2-n1+1];
for(int i=n1; i<=n2; i++) M[i-n1]= m[i];
return M;
}
// максимальный элемент массива
template <typename Type>
Type ArrayMax (Type m[], int size) {
Type max= m[0];
for (int i=1; i<size; i++) {
if ( m[i] > max ) max= m[i];
}
return max;
}
template <typename Type>
Type ArrayMin (Type m[], int size) {
Type min= m[0];
for (int i=1; i<size; i++) {
if ( m[i] < min ) min= m[i];
}
return min;
}
// сортировка по убыванию
template <typename Type>
Type* ArraySorte (Type m[], int size) {
Type *M= ArrayCopy(m,size);
Type *y= new Type [size]; // буферный массив
Type max, poz; // максимальные места
int i,j; // счётчики
for(j=0;j<size-1;j++) {
max=M[j]; poz=j;
for(i=j+1;i<size;i++) {
if(M[i]>=max) {
max=M[i]; poz=i;
}
}
y[j]=max;
for(i=j+1;i<=poz;i++) y[i]=M[i-1];
for(i=poz+1;i<size;i++) y[i]=M[i];
for(i=0;i<size;i++) M[i]=y[i];
}
delete []y;
return M;
}
// поменять порядок элементов наоборот
template <typename Type>
Type* Reverse (Type m[], int size) {
Type *M= new Type [size];
for (int i=0; i<size; i++) M[i]= m[size-1-i];
return M;
}
// сортировка по возрастанию
template <typename Type>
Type* ArraySorteReverse (Type m[], int size) {
Type *M= Reverse( ArraySorte (m,size) ,size );
return M;
}
// слияние 2 массивов
template <typename Type>
Type* Merge (Type m[], int msize, Type n[], int nsize) {
Type *M= new Type [msize+nsize];
for(int i=0; i<msize; i++) M[i]= m[i];
for(int i=0; i<nsize; i++) M[i+msize]= n[i];
return M;
}
// удаление элемента n из массива
template <typename Type>
Type* nDelete (Type m[], int size, int n) {
Type *M= new Type [size-1];
for(int i=0; i<n; i++) M[i]= m[i];
for(int i=n+1; i<size; i++) M[i-1]= m[i];
return M;
}
// вставка элемент в массив после n
template <typename Type>
Type* nInsert (Type m[], int size, int n, Type info) {
Type *M= new Type [size+1];
if(n!= -1) for(int i=0; i<=n; i++) M[i]= m[i];
M[n+1]= info;
for(int i=n+1; i<size; i++) M[i+1]= m[i];
return M;
}
// есть ли элемент в массиве
template <typename Type>
bool isInArray (Type m[], int size, Type checker) {
bool res= false;
for(int i=0; i<size; i++) {
if (m[i]==checker) res= true;
}
return res;
}
// поиск позиции элемента в массиве
template <typename Type>
int PozitionInArray (Type m[], int size, Type checker) {
if(!isInArray(m,size,checker)) return -1;
else {
int poz;
for(int i=0; i<size; i++) {
if (m[i]==checker) poz= i;
}
return poz;
}
}
// количество таких элементов в массиве
template <typename Type>
int QuantityInArray (Type m[], int size, Type checker) {
int num= 0;
for(int i=0; i<size; i++) {
if(m[i]==checker) num++;
}
return num;
}
// позиция максимального элемента в массиве
template <typename Type>
int ArrayMaxPozition (Type m[], int size) {
return PozitionInArray (m, size, ArrayMax(m,size));
}
template <typename Type>
int ArrayMinPozition (Type m[], int size) {
return PozitionInArray (m, size, ArrayMin(m,size));
}
// вывод массива
template <typename Type>
void ArrayPrint (Type m[], int size) {
cout << "(";
for(int i=0; i<size; i++) {
cout << m[i];
if(i!= size-1) cout << ", ";
}
cout << ")";
}
//////////////////////////////////////////////////////////////////////////////////////
// МАТЕМАТИКА
// сумма всех элементов массива
template <typename Type>
Type ArraySum (Type m[], int size) {
Type res= m[0];
for (int i=1; i<size; i++) res+=m[i];
return res;
}
template <typename Type>
Type ArraySum (Type m[], int size, int n) {
Type res= 0;
for (int i=n; i<size; i++) res+=m[i];
return res;
}
// произведение всех элементов массива
template <typename Type>
Type ArrayProd (Type m[], int size) {
Type res= m[0];
for (int i=1; i<size; i++) res*=m[i];
return res;
}
// скалярное произведение
template <typename Type>
Type Scalar (Type m[], Type n[], int size) {
Type res= 0;
for (int i=0; i<size; i++) res+= m[i]*n[i];
return res;
}
// модуль вектора
template <typename Type>
double Module (Type m[], int size) {
Type res= 0;
for (int i=0; i<size; i++) res+= m[i]*m[i];
return sqrt(res);
}
// нормирование (деление всех элементов на модуль) или направляющие косинусы
template <typename Type>
Type* ArrayNorm (Type m[], int size) {
Type *M= new Type [size];
double *module= new double;
*module= Module(m,size);
for(int i=0; i<size; i++) M[i]= m[i] / *module;
delete module;
return M;
}
// углы направляющих косинусов в градусах
template <typename Type>
double* ArrayNormAngle (Type m[], int size) {
Type *M= new Type [size];
double *module= new double;
*module= Module(m,size);
for(int i=0; i<size; i++) M[i]= acos (m[i] / *module) * kgrad;
delete module;
return M;
}
// в радианах
template <typename Type>
double* ArrayNormAngleRad (Type m[], int size) {
Type *M= new Type [size];
double *module= new double;
*module= Module(m,size);
for(int i=0; i<size; i++) M[i]= acos (m[i] / *module) ;
delete module;
return M;
}
// угол между 2 векторами
template <typename Type>
double BothAngle (Type m[], Type n[], int size) {
return Scalar(m,n,size) * kgrad / ( Module(m,size) * Module(n,size) ) ;
}
// в радианах
template <typename Type>
double BothAngleRad (Type m[], Type n[], int size) {
return Scalar(m,n,size) / ( Module(m,size) * Module(n,size) ) ;
}
// сумма массивов
template <typename Type>
Type* SumOfArrays (Type m[], Type n[], int size) {
Type *M= new Type [size];
for(int i=0; i<size; i++) M[i]= m[i]+n[i];
return M;
}
// произведение массивов
template <typename Type>
Type* ProdOfArrays (Type m[], Type n[], int size) {
Type *M= new Type [size];
for(int i=0; i<size; i++) M[i]= m[i]*n[i];
return M;
}
// массив из разности элементов
template <typename Type>
Type* SubstArray (Type m[], int size) {
Type *M= new Type [size-1];
for(int i=0; i<size-1; i++) M[i]= abs(m[i+1]-m[i]);
return M;
}
// массив из сумм элементов
template <typename Type>
Type* SumArray (Type m[], int size) {
Type *M= new Type [size-1];
for(int i=0; i<size-1; i++) M[i]= m[i+1]+m[i];
return M;
}
// противоположный по знаку массив
template <typename Type>
Type* Versus (Type m[], int size) {
Type *M= new Type [size];
for(int i=0; i<size; i++) M[i]= (-1) * m[i];
return M;
}
// разность массивов
template <typename Type>
Type* SubstOfArrays (Type m[], Type n[], int size) {
Type *M= new Type [size];
for(int i=0; i<size; i++) M[i]= m[i]-n[i];
return M;
}
// расстояние между 2 точками
template <typename Type>
double ArrayDistance (Type m[], Type n[], int size) {
double dist= 0;
for(int i=0; i<size; i++) dist+= (m[i]-n[i])*(m[i]-n[i]);
return sqrt(dist);
}
// умножение на число
template <typename Type>
Type* kArray (Type m[], int size, double k) {
Type *M= new Type [size];
for(int i=0; i<size; i++) M[i]= k*m[i];
return M;
}
// сумма взаимодействий 2 элементов
template <typename Type>
Type Interaction2 (Type m[], int size) {
if(size>=2) {
Type res= 0;
for(int i=0; i<size-1; i++) for(int j=i+1; j<size; j++) res+= m[i]*m[j];
return res;
}
}
template <typename Type>
Type Interaction2 (Type m[], int size, int n) {
if(size>=2) {
Type res= 0;
for(int i=n; i<size-1; i++) for(int j=i+1; j<size; j++) res+= m[i]*m[j];
return res;
}
}
template <typename Type>
Type Interaction2 (Type m[], int size, int n1, int n2) {
if(size>=2) {
Type res= 0;
for(int i=n1; i<n2-1; i++) for(int j=i+1; j<n2; j++) res+= m[i]*m[j];
return res;
}
}
// cумма взаимодействия p элементов
template <typename Type>
Type Interaction (int p, Type m[], int size, int n) {
if(size>=p) {
if(p==1) return ArraySum(m,size,n);
Type res= 0;
for(int i=n; i<size-p+1; i++) res+= m[i] * Interaction(p-1,m,size,i+1) ;
return res;
}
}
template <typename Type>
Type Interaction (int p, Type m[], int size) {
int n= 0;
if(size>=p) {
if(p==1) return ArraySum(m,size,n);
Type res= 0;
for(int i=n; i<size-p+1; i++) res+= m[i] * Interaction(p-1,m,size,i+1) ;
return res;
}
}
// среднее арифметическое
template <typename Type>
double Average (Type m[], int size) {
return ArraySum(m,size) / size ;
}
// среднее геометрическое
template <typename Type>
double AverageG (Type m[], int size) {
return pow (ArrayProd(m,size) , 1/size) ;
}
// среднее квадратичное
template <typename Type>
double AverageQ (Type m[], int size) {
return Module(m,size);
}
// дисперсия
template <typename Type>
double Dispers (Type m[], int size) {
double res= 0;
for(int i=0; i<size; i++) res+= pow(m[i] - Average(m,size) , 2);
return sqrt(res);
}
// среднее гармоническое
template <typename Type>
double AverageH (Type m[], int size) {
if(ArrayProd(m,size)!=0) {
double res= 0;
for(int i=0; i<size; i++) res+= 1 / m[i];
return 1/res;
}
}
//////////////////////////////////////////////////////////////////////////////////////
// ФОРМИРОВАНИЕ УДОБНЫХ МАССИВОВ
// из нулей
int* Zeros (int size) {
int *M= new int [size];
for(int i=0; i<size; i++) M[i]= 0;
return M;
}
// из единиц
int* Ones (int size) {
int *M= new int [size];
for(int i=0; i<size; i++) M[i]= 1;
return M;
}
// из порядков
int* Pozitions (int size) {
int *M= new int [size];
for(int i=0; i<size; i++) M[i]= i;
return M;
}
// из порядков +1
int* Pozitions1 (int size) {
int *M= new int [size];
for(int i=0; i<size; i++) M[i]= i+1;
return M;
}
// из случайных чисел
int* Randoms (int size, int a, int b) {
int *M= new int [size];
for(int i=0; i<size; i++) M[i]= rand()%(b-a+1)+a;
return M;
}
int* Randoms (int size) {
int *M= new int [size];
for(int i=0; i<size; i++) M[i]= rand()%19-9;
return M;
}
________________________________________________
РАЗДЕЛ 2
Дата добавления: 2018-09-22; просмотров: 201; Мы поможем в написании вашей работы! |
Мы поможем в написании ваших работ!