Форум программистов, компьютерный форум CyberForum.ru

Коллекция алгоритмов от Johna Smith - C++

Войти
Регистрация
Восстановить пароль
Другие темы раздела
C++ Работа с графом (Требуется по заявке клиента предложить способы обмена жилплощади) http://www.cyberforum.ru/cpp-beginners/thread870941.html
В файле записаны предложения по обмену жилплощадью. Имеются варианты размена одной квартиры на две других либо на квартиру и комнату. Требуется по заявке клиента предложить способы обмена. Предусмотреть возможность нахождения обменов, в которых участвуют более двух сторон.
C++ Определить, является ли текст десятичной записью числа, кратного 9 Является ли текст записью десятичного числа,кратного 9 В заданный непустой текст входят только цифры и буквы. Определить. удовлетворяет ли он следующему свойству: 1) текст является десятичной записью числа, кратного 9; http://www.cyberforum.ru/cpp-beginners/thread870934.html
Удаление первых n элементов из vector C++
Почему, к примеру, если k=3 а pop=2, то студия выдаст ошибку(итератор вне допустимого диапазона) при запуске функции erase. По моей логике, необходимо было удалить первый элемент. #include <cstdio> #include <iostream> #include <vector> #define pb push_back #define ull unsigned long long using namespace std; vector<int> t;
Найти сумму квадратов элементов матрицы C++
Помогите пожалуйста!)
C++ Получить целочисленную матрицу http://www.cyberforum.ru/cpp-beginners/thread870928.html
Задание ниже: Nastik23, оформите тему в соответствии с правилами форума: текстовые задания набирайте от руки
C++ Считывание структур из файла #include <stdio.h> #include <string.h> #include <malloc.h> #define Lmax 20 struct student { struct { char name,surname,patronymic; } fio; подробнее

Показать сообщение отдельно
LK
Заблокирован
19.05.2013, 15:58  [ТС]     Коллекция алгоритмов от Johna Smith
Вычитание
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (substraction)
//  (c) Johna Smith, 1996
//
//  Given: A,B - matrixes M x N
//  This algorithm substract these matrixes and display the result.
//
//////////////////////////////////////////////////////////////////////////////
#define M   5
#define N   6
#include <stdio.h>
int a[N][M]={{4,4,4,4,4},{5,4,3,2,1},{1,2,3,4,5},{5,4,3,2,1},{1,2,3,4,5},{0,1,0,1,0}};
int b[N][M]={{3,4,4,4,4},{5,3,3,2,1},{1,2,2,4,5},{5,4,3,1,1},{1,2,3,4,4},{0,1,0,1,0}};
int c[N][M];
void show_matrix(int *matrix) // this functions displays matrix
{
 for(int i=0;i<N;i++)
 {
    for(int j=0;j<M;j++)
      printf("%d ",*(matrix+i*M+j));
    printf("\n");
 }
 printf("\n");
}
void main(void)
{
  // display matrixes A and B
  show_matrix(a[0]);
  show_matrix(b[0]);
  // substract these matrixes
  for(int i=0;i<N;i++)
    for(int j=0;j<M;j++)
      c[i][j]=a[i][j]-b[i][j];
  // display the difference
  show_matrix(c[0]);
}

Транспонирование
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (transponating)
//  (c) Johna Smith, 1996
//
//  Given: A (N x N),                                 T
//  This algorithm transponates matrix A and returns A  .
//          T
//   A(ji)=A (ij)
//
//////////////////////////////////////////////////////////////////////////////
#define N  3
#include <stdio.h>
float a[N][N]={{4,8,0},{8,8,8},{2,0,1}};
void Transponate(float *matrix)
{
  float swp;
  for (int i=0;i<N;i++)
    for (int j=i+1;j<N;j++)
    {
      swp=*(matrix+i*N+j);
      *(matrix+i*N+j)=*(matrix+j*N+i);
      *(matrix+j*N+i)=swp;
    }
}
void show_matrix(float *matrix) // this functions displays matrix
{
 for(int i=0;i<N;i++)
 {
   for(int j=0;j<N;j++)
     printf("%f ",*(matrix+i*N+j));
   printf("\n");
 }
 printf("\n");
}
void main(void)
{
  // display matrix A
  show_matrix(a[0]);
  // Transponate it
  Transponate(a[0]);
  // display transponated matrix
  show_matrix(a[0]);
  // Transponate it again
  Transponate(a[0]);
  // display the transponation of the transponated matrix
  show_matrix(a[0]);
}


Операции с комплексными величинами

Сложение и вычитание в алгебраической форме
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (addition and substraction in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1+z2=(a1+a2)+i*(b1+b2)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
struct complex
{
  float re;
  float im;
};
void show_complex(complex c) // this functions displays complex value
{
  printf("%f%+fi",c.re,c.im);
}
complex Add(complex a, complex b)
{
  complex c;
  c.re=a.re+b.re;
  c.im=a.im+b.im;
  return c;
}
complex Substract(complex a, complex b)
{
  complex c;
  c.re=a.re-b.re;
  c.im=a.im-b.im;
  return c;
}
complex a={2,3},b={4,6};
void main(void)
{
  // addition
  show_complex(a);
  printf(" + ");
  show_complex(b);
  printf(" = ");
  show_complex(Add(a,b));
  // substraction
  printf("\n");
  show_complex(b);
  printf(" - ");
  show_complex(a);
  printf(" = ");
  show_complex(Substract(b,a));
}

Деление в алгебраической форме
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (division in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1/z2=(a1*a2-b1*b2)/(a2^2+b2^2)+i*(b1*a2-a1*b2)/(a2^2+b2^2)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
struct complex
{
  float re;
  float im;
};
void show_complex(complex c) // this functions displays complex value
{
  printf("%f%+fi",c.re,c.im);
}
complex Div(complex a, complex b)
{
  complex c;
  c.re=(a.re*b.re+a.im*b.im)/(b.re*b.re+b.im*b.im);
  c.im=(a.im*b.re-a.re*b.im)/(b.re*b.re+b.im*b.im);
  return c;
}
complex a={2,3},b={-1,2};
void main(void)
{
  show_complex(a);
  printf(" / ");
  show_complex(b);
  printf(" = ");
  show_complex(Div(a,b));
}

Умножение в алгебраической форме
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (multiplication in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1*z2=(a1*a2-b1*b2)+i*(a1*b2+a2*b1)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
struct complex
{
  float re;
  float im;
};
void show_complex(complex c) // this functions displays complex value
{
  printf("%f%+fi",c.re,c.im);
}
complex Mul(complex a, complex b)
{
  complex c;
  c.re=a.re*b.re-a.im*b.im;
  c.im=a.im*b.re+a.re*b.im;
  return c;
}
complex a={2,3},b={-1,2};
void main(void)
{
  show_complex(a);
  printf(" * ");
  show_complex(b);
  printf(" = ");
  show_complex(Mul(a,b));
}

Деление в экспоненциальной форме
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (division in M*exp(i*phi) form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=M1*exp(i*phi1), z2=M2*exp(i*phi2)
//  z1/z2=M1/M2*exp(i*(phi1-phi2))
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define Pi      3.1415926536
struct exp_complex
{
  float M;
  float phi;
};
void show_exp_complex(exp_complex c) // this functions displays complex value
{
  printf("%f*exp(%f*i)",c.M,c.phi);
}
exp_complex Div(exp_complex a,exp_complex b)
{
  exp_complex c;
  c.M=a.M/b.M;
  c.phi=a.phi-b.phi;
  return c;
}
exp_complex a={1,0}, b={-1,-Pi/2};
void main(void)
{
  show_exp_complex(a);
  printf(" / ");
  show_exp_complex(b);
  printf(" = ");
  show_exp_complex(Div(a,b));
}

Умножение в экспоненциальной форме
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (multiplication in M*exp(i*phi) form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=M1*exp(i*phi1), z2=M2*exp(i*phi2)
//  z1*z2=M1*M2*exp(i*(phi1+phi2))
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define Pi      3.1415926536
struct exp_complex
{
  float M;
  float phi;
};
void show_exp_complex(exp_complex c) // this functions displays complex value
{
  printf("%f*exp(%f*i)",c.M,c.phi);
}
exp_complex Mul(exp_complex a,exp_complex b)
{
  exp_complex c;
  c.M=a.M*b.M;
  c.phi=a.phi+b.phi;
  return c;
}
exp_complex a={-1,Pi/2}, b={-1,-Pi/2};
void main(void)
{
  show_exp_complex(a);
  printf(" * ");
  show_exp_complex(b);
  printf(" = ");
  show_exp_complex(Mul(a,b));
}


Математическая статистика

Интервальная оценка среднего (рассмотрены случаи значений доверительной вероятности 0.95 и 0.99)
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating <x>, (x-dx,x+dx)
//  (c) Johna Smith, 1996
//
//  Method description:
//    dx=t(alpha,n)*Sn/sqrt(n)
//    <x>=Sum x/n
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N   6 // number of experiments
float x[N]={3,5,4,5,3,4}; // experimental data
double Student095(unsigned int n)
{
  double result;
  result=1.96+2.4/(n-1)+5.901610*pow(n-1,-2.372);
  return result;
}
double Student099(unsigned int n)
{
  double result;
  result=2.576+5.0/(n-1)+29.12178*pow(n-1,-2.591843);
  return result;
}
// this function calculates average x
double Average(void)
{
  double sum=0;
  for (int i=0;i<N;i++) sum+=x[i];
  return sum/N;
}
// this function calculates delta x
double DeltaX(void)
{
  double ave,sum=0,Sn;
  ave=Average();
  for (int i=0;i<N;i++) sum+=x[i]*x[i];
  Sn=sqrt(((double)N/(N-1))*((double)sum/N-ave*ave));
  return Student095(N)*Sn/sqrt(N); // you can use Student099 if you want
}
void main(void)
{
  printf("Data: ");
  for (int i=0;i<N;i++) printf("%f ",x[i]);
  printf("\nx=%f +/- %f",Average(),DeltaX());
}

Вычисление коэффициента Стьюдента
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating Student's coefficients t(alpha,n), if alpha=0.95 or 0.99
//  (c) Johna Smith, 1996
//
//  Method description:
//   To calculate Student's coefficients we use the following
//  approximation:
//    t(n)=C0+C1/(n-1)+C2(n-1)^(-C3)
//  Coefficients Ci for alpha=0.95 and alpha=0.99 is given
//  n - number of experiments
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
double Student095(unsigned int n)
{
  double result;
  result=1.96+2.4/(n-1)+5.901610*pow(n-1,-2.372);
  return result;
}
double Student099(unsigned int n)
{
  double result;
  result=2.576+5.0/(n-1)+29.12178*pow(n-1,-2.591843);
  return result;
}
void main(void)
{
  printf("alpha=0.95 n=5 t5= %g\n",Student095(5));
  printf("alpha=0.99 n=15 t15= %g",Student099(15));
}

Теоретический коэффициент корреляции для двух величин
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating theoretical corellation coefficient for two values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between two values
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N       3 // number of experiments
double x[N]={1,2,3};
double y[N]={2,4,5};
double TheoreticalCorellation(void)
{
  double sxy=0,sx2=0,sy2=0;
  for (int i=0;i<N;i++)
  {
    sxy+=x[i]*y[i];
    sx2+=x[i]*x[i];
    sy2+=y[i]*y[i];
  }
  return sxy/sqrt(sx2*sy2);
}
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%f %f\n",x[i],y[i]);
  printf("\nCorellation coefficient r= %f\n",TheoreticalCorellation());
}

Сглаживание экспериментальных данных по 3 точкам
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 3 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 3-points groups
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N 10      // number of experiments
float data[N]={0.9,2.12,2.92,4.15,4.9,6.1,6.92,8.15,9.05,9.8};
float smoothed[N];
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%.4f ",data[i]);
  //smoothing
  smoothed[0]=(5*data[0]+2*data[1]-data[2])/6;
  for (i=1;i<N-1;i++) smoothed[i]=(data[i-1]+data[i]+data[i+1])/3;
  smoothed[N-1]=(5*data[N-1]+2*data[N-2]-data[N-3])/6;
  printf("\nSmoothed data:\n");
  for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}

Сглаживание экспериментальных данных по 5 точкам
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 5 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 5-points groups
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N  10      // number of experiments
float data[N]={0.9,2.12,2.92,4.15,4.9,6.1,6.92,8.15,9.05,9.8};
float smoothed[N];
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%.4f ",data[i]);
  //smoothing
  smoothed[0]=(3*data[0]+2*data[1]+data[2]-data[4])/5;
  smoothed[1]=(4*data[0]+3*data[1]+2*data[2]+data[3])/10;
  for (i=2;i<N-2;i++) smoothed[i]=(data[i-2]+data[i-1]+data[i]+
                                   data[i+1]+data[i+2])/5;
  smoothed[N-2]=(data[N-4]+2*data[N-3]+3*data[N-2]+4*data[N-1])/10;
  smoothed[N-1]=(3*data[N-1]+2*data[N-2]+data[N-3]-data[N-5])/5;
  printf("\nSmoothed data:\n");
  for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}

Сглаживание экспериментальных данных по 7 точкам
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 7 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 7-points groups
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N 10      // number of experiments
float data[N]={0.7,1.08,1.39,1.64,1.76,1.99,2.04,2.22,2.28,2.42};
float smoothed[N];
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%.4f ",data[i]);
  //smoothing
  smoothed[0]=(39*data[0]+8*data[1]-4*(data[2]+data[3]-data[5]) + data[4]-2*data[6])/42;
  smoothed[1]=(8*data[0]+19*data[1]+16*data[2]+6*data[3]-4*data[4] - 7*data[5]+4*data[7])/42;
  smoothed[2]=(-4*data[0]+16*data[1]+19*data[2]+12*data[3]+2*data[4] - 4*data[5]+data[6])/42;
  for (i=3;i<N-3;i++)
    smoothed[i]=(7*data[i]+6*(data[i-1]+data[i+1]) + 
          3*(data[i-2]+data[i+2])-2*(data[i-3]+data[i+3]))/21;
  smoothed[N-3]=(data[N-7]-4*data[N-6]+2*data[N-5]+12*data[N-4] + 
          9*data[N-3]+16*data[N-2]-4*data[N-1])/42;
  smoothed[N-2]=(4*data[N-7]-7*data[N-6]-4*data[N-5]+6*data[N-4] + 
          16*data[N-3]+19*data[N-2]+8*data[N-1])/42;
  smoothed[N-1]=(-2*data[N-7]+4*data[N-6]+data[N-5]-4*data[N-4] - 
          4*data[N-3]+8*data[N-2]+39*data[N-1])/42;
  printf("\nSmoothed data:\n");
  for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}

Эмпирический коэффициент корреляции для двух величин
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating empirical corellation coefficient for two values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between two values
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N 3 // number of experiments
double x[N]={1,2,3};
double y[N]={2,4,5};
double EmpiricalCorellation(void)
{
  double sxy=0,sx2=0,sy2=0,sx=0,sy=0;
  for (int i=0;i<N;i++)
  {
    sxy+=x[i]*y[i];
    sx2+=x[i]*x[i];
    sy2+=y[i]*y[i];
    sx+=x[i];
    sy+=y[i];
  }
  return (sxy-sx*sy/N)/(sqrt(sx2-sx*sx/N)*sqrt(sy2-sy*sy/N));
}
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%f %f\n",x[i],y[i]);
  printf("\nCorellation coefficient r= %f\n",EmpiricalCorellation());
}

Эмпирический коэффициент корреляции для трех величин
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating corellation coefficient for three values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between three values
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N   3 // number of experiments
double x[N]={1,2,3};
double y[N]={2,4,5};
double z[N]={1,7,0};
double EmpiricalCorellation2D(double *i, double *j)
{
  double sij=0,si2=0,sj2=0,si=0,sj=0;
  for (int a=0;a<N;a++)
  {
    sij+=*(i+a)**(j+a);
    si2+=*(i+a)**(i+a);
    sj2+=*(j+a)**(j+a);
    si+=*(i+a);
    sj+=*(j+a);
  }
  return (sij-si*sj/N)/(sqrt(si2-si*si/N)*sqrt(sj2-sj*sj/N));
}
double Corellation3D(void)
{
  double rxz,ryz,rxy;
  // calculating partial corellation coefficients
  rxy=EmpiricalCorellation2D(x,y);
  rxz=EmpiricalCorellation2D(x,z);
  ryz=EmpiricalCorellation2D(y,z);
  return sqrt((rxz*rxz-ryz*ryz-2*rxy*rxz*ryz)/(1-rxy*rxy));
}
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%f %f %f\n",x[i],y[i],z[i]);
  printf("\nCorellation coefficient R= %f\n",Corellation3D());
}

Частный коэффициент корреляции для трех величин
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating partial corellation coefficient for three values
//  (c) Johna Smith, 1996
//
//  Method description:
//   Partial corellation coefficient shows dependencies between two values
//   using set of three values (we exclude influence of one of these values)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// excluding influence of z, recieving corellation
// coefficients between x,y  y,z and x,z
double PartialCorellation3D(double rxy,double rxz,double ryz)
{
  return (rxy-rxz*ryz)/sqrt((1-rxz*rxz)*(1-ryz*ryz));
}
void main(void)
{
  double rxy=0.998753,rxz=0.69753,ryz=0.7473;
  printf("Data:\n");
  printf("%f %f %f\n",rxy,rxz,ryz);
  printf("\nPartial corellation coefficient r= %.10f\n",
  PartialCorellation3D(rxy,rxz,ryz));
}


Добавлено через 8 минут
Интегрирование функций

Метод трапеций
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (trapecias method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We approximate area under function with trapecias, using
//   the same intervals by the X axe.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
double f(double x)
{
  // !!! you should specify your own function here
  return x;
}
double Integrate(double x1, double x2, double step)
{
  double x=x1;
  double I=0; // integral value
  while (x<x2-step)
  {
    // calculating trapecia's area as f(x)+f(x+step) * halfheight of it
    I += (f(x)+f(x+step))*step/2;
    x += step;
  }
  return I;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e-4));
}

Метод Симпсона
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (Simpson's method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We approximate area under function with parabolas, using
//    three points. We repeat calculation with less step if
//    the diffference between I(h) and I(h/2) is greater than eps
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
double f(double x)
{
  // you should specify your own function here
  return x*x*sin(x);
}
double Integrate(double x1, double x2, double step, double precision)
{
  double x=x1;
  double I1=0,I2=0; // integral value
  double s;
  do
  {
    I1=I2;
    s=(f(x1)-f(x2))/2;
    x=x1+step;
    while (x<x2)
    {
      s+=2*f(x)+f(x+step);
      x+=2*step;
    }
    I2=2*step*s/3;
    step/=2.0;  // try once more using less step
  } while (fabs(I1-I2)>precision);
  return I2;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e-4,1e-6));
}

Метод прямоугольников
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (rectanges method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We approximate area under function with rectangles, using
//    the same intervals by the X axe.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
double f(double x)
{
  // you should specify your own function here
  return x;
}
double Integrate(double x1, double x2, double step)
{
  double x=x1;
  double I=0; // integral value
  
  while (x<x2-step)
  {
    // calculating function at the halfpoint of interval
    I+=step*f(x+step/2);
    x+=step;
  }
  
  return I;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e-4));
}

Метод Монте-Карло
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (Monte-Karlo method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We generate random numbers from interval [a;b] where we want to
//   integrate our function and sum function values in this random points.
//   Here we use approximation of the integral:
//   I=<y>*sigma, where <y> is average value of the function f on [a;b]
//   and sigma is volume of integrating area (if we integrate by interval it
//   is its length)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
double f(double x)
{
  // you should specify your own function here
  return x*x*sin(x);
}
double Integrate(double x1, double x2, long int steps)
{
  double I=0; // integral value
  randomize();  // initializing random numbers generator
  for (long int i=0;i<steps;i++)
    I+=f((random(10000)/10000.0)*(x2-x1)+x1);
  I*=(x2-x1)/steps;
  return I;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e5));
}


Календарь

Число дней между двумя датами
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating number of days between two dates
//  (c) Johna Smith, 1996
//
//  Method description:
//    We calculte two Julian dates and substract them
//    Result will be number of days between two dates
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
enum Month {January=1,February,March,April,May,June,July,August,
            September,October,November,December};
struct Date
{
  float day;
  Month month;
  int year;
};
double JulianDate(Date date)
{
  long int A,B,Y,M;
  if (date.month==January || date.month==February)
  {
    Y=date.year-1;
    M=date.month+12;
  } else
  {
    Y=date.year;
    M=date.month;
  }
  if ((date.year>1581) ||
      (date.year==1581 && date.month>October) ||
      (date.year==1581 && date.month==October && date.day>=15))
  {
    A=Y/100;
    B=2-A+A/4;
  } else B=0;
  return B+(long int)(365.25*Y)+(long int)(30.6001*(M+1))+date.day+1720994.5;
}
Date Revolution={17.5,October,1917};
Date Today={11,October,1996};
void main(void)
{
  printf("Days between %f/%d/%d and %f/%d/%d = %f",
         Revolution.day,Revolution.month,Revolution.year,
         Today.day,Today.month,Today.year,
         JulianDate(Today)-JulianDate(Revolution));
}

Определение дня недели по дате
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating day of week by date
//  (c) Johna Smith, 1996
//
//  Method description:
//  Here we use special formulas different for January & February
//  and other monthes. Use this algorithm only for dates after
//  October, 15, 1528
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
enum Month {January=1,February,March,April,May,June,July,August,
            September,October,November,December};
struct Date
{
  float day;
  Month month;
  int year;
};
int DayOfWeek(Date date)
{
  float F;
  if (date.month<March)
    F=365*date.year+date.day+31*(date.month-1)+
      (int)((date.year-1)/4)-(int)(3*(int)((date.year-1)/100+1)/4);
  else
    F=365*date.year+date.day+31*(date.month-1)-(int)(0.4*date.month+2.3)+
      (int)(date.year/4)-(int)(3*(int)(date.year/100+1)/4);
  return  F-7*(int)(F/7)-1;
}
Date Today={11,October,1996};
void main(void)
{
  printf("%f/%d/%d is ",Today.day,Today.month,Today.year);
  switch (DayOfWeek(Today))
  {
    case -1:printf("Sunday\n");break;
    case 0:printf("Monday\n");break;
    case 1:printf("Tuesday\n");break;
    case 2:printf("Wednesday\n");break;
    case 3:printf("Thursday\n");break;
    case 4:printf("Friday\n");break;
    case 5:printf("Saturday\n");break;
  }
}

Вычисление даты по юлианскому календарю
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
//////////////////////////////////////////////////////////////////////////////
//
//  Julian date calculating
//  (c) Johna Smith, 1996
//
//  Method description:
//    Julian date is the number of days since afternoon (GMT) of
//    January, 1, 4713 BC.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
enum Month {January=1,February,March,April,May,June,July,August,
            September,October,November,December};
struct Date
{
  float day;
  Month month;
  int year;
};
double JulianDate(Date date)
{
  long int A,B,Y,M;
  if (date.month==January || date.month==February)
  {
    Y=date.year-1;
    M=date.month+12;
  } else
  {
    Y=date.year;
    M=date.month;
  }
  if ((date.year>1581) ||
      (date.year==1581 && date.month>October) ||
      (date.year==1581 && date.month==October && date.day>=15))
  {
    A=Y/100;
    B=2-A+A/4;
  } 
  else 
    B=0;
  return B+(long int)(365.25*Y)+(long int)(30.6001*(M+1))+date.day+1720994.5;
}
Date Revolution={17.5,October,1917};
void main(void)
{
  printf("%f/%d/%d = %f",Revolution.day,Revolution.month,
                         Revolution.year,JulianDate(Revolution));
}


Решение систем линейных уравнений

Метод квадратного корня
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (square root method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    From given matrix A we build two auxulary triangular matrixes S & D
//    Then solving equations Gy=b (where g[I][j]=s[j][I]*d[j]) and Sx=y
//    we obtain vector x.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N    4     // size of matrix
#define N1   N+1
float matrix[N][N1]=
    {{10.9,1.2,2.1,0.9,23.2},
  {1.2,11.2,1.5,2.5,38.1},
  {2.7,1.5,9.8,1.3,40.3},
  {0.9,2.5,1.3,12.1,58.2}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  register char i,j,k;
  float s[N][N],x[N],y[N],d[N],sum;
  // Printing given matrix
  ShowMatrix();
  // Building matrixes S and D
  for (j=0;j<N;j++)
    for (i=0;i<=j;i++)
    {
      sum=matrix[i][j];
      for (k=0;k<i;k++)
        sum-=s[k][i]*d[k]*s[k][j];
      if (i==j)
      {
        d[i]=(sum>0?1:-1);
        s[i][j]=sqrt(fabs(sum));
      }
      else s[i][j]=sum/(d[i]*s[i][i]);
    }
  // Solving equation Gy=b (G: g[I][j]=s[j][I]*d[j])
  for (i=0;i<N;i++)
  {
    y[i]=matrix[i][N];
    for (j=0;j<i;j++) y[i]-=s[j][i]*d[j]*y[j];
    y[i]/=s[i][i]*d[i];
  }
  // Solving equation Sx=y
  for (i=N-1;i>=0;i--)
  {
    x[i]=y[i];
    for (j=i+1;j<N;j++) x[i]-=s[i][j]*x[j];
    x[i]/=s[i][i];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
  printf("x%d=%f\n",i+1,x[i]);
}

Метод Некрасова
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Nekrasov method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This is an iterative method of solving systems of linear equations.
//    Criteria for end of iterations is ||xn||-||x(n-1)||<epsilon, where
//    ||x|| is norm of vector x.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N     4     // size of matrix
#define N1    N+1
float matrix[N][N1]=
      {{10.9,1.2,2.1,0.9,23.2},
       {1.2,11.2,1.5,2.5,38.1},
       {2.7,1.5,9.8,1.3,40.3},
       {0.9,2.5,1.3,12.1,58.2}
     };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
// this function calculates ||x||
float norma(float *x)
{
  float sum=0;
  for (unsigned i=0;i<N;i++) sum+=fabs(*(x+i));
  return sum;
}
void main(void)
{
  // Variables declaration
  register char i,j;
  float x_current[N],x_next[N],sum1,sum2,epsilon;
  epsilon=1e-7;
  // Printing given matrix
  ShowMatrix();
  // Solving
  do
  {
    for (i=0;i<N;i++) x_current[i]=x_next[i];
    for (i=0;i<N;i++)
    {
      sum1=0;sum2=0;
      for(j=0;j<i;j++) sum1+=matrix[i][j]*x_next[j];
      for(j=i+1;j<N;j++) sum2+=matrix[i][j]*x_current[j];
      x_next[i]=(matrix[i][N]-sum1-sum2)/matrix[i][i];
    }
  } while (fabs(norma(x_current)-norma(x_next))>epsilon);
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x_current[i]);
}

Метод 3-диагональных матриц
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (3-diagonal matrix)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is developed for 3-diagonal matrixes:
//    (a11 a12 0 ..............)
//    (a21 a22 a23 0 ..........)
//    (0   a32 a33 a34 0 ......) = M
//    (........................)
//    (.......... 0 an(n-1) ann)
//    Cause all elements except three diagonals is 0 we'll store
//    only non zero elements in the A,B,C arrays. Array D stores
//    vector B (M*X=B). We build two auxulary vectors P and Q and
//    then calculate vector X using the following recurrent formula:
//      X(k-1)=P(k)*X(k)+Q(k)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#define N    20
#define N1   N+1
void main(void)
{
  // Variables declaration
  register int i,i1,k;  // using processor registers for better performance
  float a[N]={0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
  float b[N]={2,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,2};
  float c[N]={-1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
  float d[N]={0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
  float x[N1],p[N1],q[N1],tmp;
  // Some initializations
  p[0]=0;
  q[0]=0;
  x[N]=0;
  c[N-1]=0; // to avoid overflow
  // Building vectors P and Q
  for (i=0;i<N;i++)
  {
    i1=i+1;
    if ((tmp=b[i]+a[i]*p[i])==0)
    {
      printf("Wrong method for such system...\n");
      return;
    };
    p[i1]=-c[i]/tmp;
    q[i1]=(d[i]-a[i]*q[i])/tmp;
  }
  // Building vector X
  for (k=N;k>0;k--) x[k-1]=p[k]*x[k]+q[k];
  // Printing solution
  printf("Solution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

Метод Гаусса
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This methos based on excluding variables from equations.
//    Using first equation we exclude x1 from all other equations
//    Using second equation equation we exclude x2 from all equations
//   (excluding first) and so on. So we have triangular matrix from which
//    we can easily get vector X. (AX=B)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N     4     // size of matrix
#define N1    N+1
float matrix[N][N1]=
                  {{2,5,4,1,28},
                   {1,3,2,1,17},
                   {2,10,9,7,77},
                   {3,8,9,2,54}
                  };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float tmp,x[N1];
  register short int i,j,k;
  // Printing given matrix
  ShowMatrix();
  // Main loop
  for (i=0;i<N;i++)
  {
    // Excluding variable x[i] from equations
    tmp=matrix[i][i];
    for (j=N;j>=i;j--) matrix[i][j]/=tmp;
    for (j=i+1;j<N;j++)
    {
      tmp=matrix[j][i];
      for (k=N;k>=i;k--)
        matrix[j][k]-=tmp*matrix[i][k];
    }
  }
  // Calculating vector x
  x[N-1]=matrix[N-1][N];
  for (i=N-2;i>=0;i--)
  {
    x[i]=matrix[i][N];
    for (j=i+1;j<N;j++) x[i]-=matrix[i][j]*x[j];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

Метод Гаусса с выбором максимального элемента
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method with max element selecting)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This methos based on excluding variables from equations.
//    Using first equation we exclude x1 from all other equations
//    Using second equation equation we exclude x2 from all equations
//    (excluding first) and so on. So we have triangular matrix from which
//    we can easily get vector X. (AX=B)
//    'Selecting max element' means that on each step we select equation
//    with max element in current column and use it as next equation to exclude
//    variable.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N      4     // size of matrix
#define N1     N+1
float matrix[N][N1]=
    {{2,5,4,1,28},
  {1,3,2,1,17},
   {2,10,9,7,77},
  {3,8,9,2,54}
        };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float max,swp,tmp,x[N1];
  register short int row_with_max_element,i,j,k;
  // Printing given matrix
  ShowMatrix();
  // Main loop
  for (i=0;i<N;i++)
  {
    // Searching for max element in the current column (i)
    max=fabs(matrix[i][i]);
    row_with_max_element=i;
    for (j=i+1;j<N;j++)
      if (fabs(matrix[j][i])>max)
      {
        row_with_max_element=j;
        max=fabs(matrix[j][i]);
      }
    // Swapping 2 lines of matrix - row_with_max_element & i
    if (row_with_max_element!=i)
    for (j=0;j<N1;j++)
    {
      swp=matrix[row_with_max_element][j];
      matrix[row_with_max_element][j]=matrix[i][j];
      matrix[i][j]=swp;
    }
    // Excluding variable x[i] from equations
    tmp=matrix[i][i];
    for (j=N;j>=i;j--) matrix[i][j]/=tmp;
    for (j=i+1;j<N;j++)
    {
      tmp=matrix[j][i];
      for (k=N;k>=i;k--)
        matrix[j][k]-=tmp*matrix[i][k];
    }
  }
  // Calculating vector x
  x[N-1]=matrix[N-1][N];
  for (i=N-2;i>=0;i--)
  {
    x[i]=matrix[i][N];
    for (j=i+1;j<N;j++) x[i]-=matrix[i][j]*x[j];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

Метод итераций
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (iterations method)
//  (c) Johna Smith, 1996
//
//  Method description:
//  This program automatically transforms system of equations to
//  iterative form and solves it.
//  Iterative form:
//     x1=a11*x1+a12*x2+... + a1n*xn+b1
//     x2=a21*x1+a22*x2+... + a2n*xn+b2
//     ...
//     xn=an1*x1+an2*x2+... + ann*xn+bn
//  Setting first iteration of vector X we can get next iteration from
//  such form of system of equations. Iterations will be finished when
//  difference between two consequitive iterations will be less than epsilon
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N      3     // size of matrix
#define N1     N+1
float matrix[N][N1]=
    {{14.38,-2.41,1.39,5.86},
  {1.84,25.36,-3.31,-2.28},
  {2.46,-3.49,16.37,4.47}
    };
int maxiterations=10;  // maximum number of iterations
float epsilon=0.0001;  // required accuracy
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float x[N],y[N],t;
  register short int i,j,k;
  int iterations=0;
  // Printing given matrix
  ShowMatrix();
  // setting first iteration of vector X
  for (i=0;i<N;i++) x[i]=matrix[i][N];
  do
  {
    for (i=0;i<N;i++)
    {
      t=-matrix[i][N];
      for (j=0;j<N;j++)
        t+=matrix[i][j]*x[j];
        y[i]=(-t+matrix[i][i]*x[i])/matrix[i][i];
    }
    iterations++;
    k=0;
    // checking solution
    while (fabs(x[k]-y[k])<epsilon && k<N) k++;
    // new iteration becomes old
    for(i=0;i<N;i++) x[i]=y[i];
  } while (k!=N && iterations<maxiterations);
  if (iterations==maxiterations)
  {
    printf("Iterations are very slow...");
  } else
  {
    // Printing solution
    printf("\nSolution:\n");
    for (i=0;i<N;i++)
      printf("x%d=%f\n",i+1,x[i]);
    printf("%d iterations were made",iterations);
  }
}

Метод Монте-Карло
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Monte-Karlo method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We should translate given system of linear equations to special form:
//      x1=a(11)x1+a(12)x2+...+a(1n)xn+a(1n+1)
//      x2=a(21)x2+a(22)x2+...+a(2n)xn+a(2n+1)
//      ....
//      xn=a(n1)x1+a(n2)x2+...+a(nn)xn+a(nn+1)
//            n
//     where SUM |aij|<1 (i=1, 2, ... ,n)
//           j=1
//    Lets divide interval [0;1] into N+1 smaller intervals.
//    Imagine a particle that is moving randomly at [0;1] interval.
//    Remember its moves until it reachs one of the bounds of the interval.
//    So we'll get a trajectory of the particle s(i),s(j),s(k),s(l),s(m),...
//    If particle moves from interval s(i) to s(j) we'll write v(ij)
//    y[i]=v(ij)*v(jk)*...*v(tm)*w(m), v(ij)=sign(aij), v(jk)=sign(ajk)
//                    n
//   w(m)=a(mn+1)/(1-SUM a(mj)),   x[i]=y/M, where M is number of particle runs
//                   j=1
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N     3     // size of matrix
#define N1    N+1
#define M     1000 // number of tries
float matrix[N][N1]=
    {{3,1,-1,-2},
  {1,4,-2,-3},
  {2,-1,5,-15}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    printf("x%d=",i+1);
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],j+1);
    printf("%+f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float b[N][N],w[N],y,c,x[N],l;
  register short int i,j,v;
  unsigned char finish;
  int row,   // currently analysed row
      tries; // number of tries to calculate true root by generating randoms
  // Transforming matrix
  for (i=0;i<N;i++)
  {
    l=0;
    for (j=0;j<N1;j++) l+=fabs(matrix[i][j]);
    v=(matrix[i][i]>0?-1:1);
    for (j=0;j<N1;j++) matrix[i][j]=v*(matrix[i][j])/l+(i==j?1:0);
  }
  // Printing transformed matrix
  ShowMatrix();
  // filling matrix B
  for (i=0;i<N;i++)
  {
    b[i][0]=fabs(matrix[i][0]);
    for(j=1;j<N;j++)
      b[i][j]=b[i][j-1]+fabs(matrix[i][j]);
  }
  // filling matrix w
  for(i=0;i<N;i++) w[i]=matrix[i][N]/(1-b[i][N-1]);
  // for each equation in the system
  for(i=0;i<N;i++)
  {
    tries=0;
    row=i;
    y=0;
    v=1;
    while (tries<M)
    {
      // generating random number
      c=random(32767)/32767.0;
      j=N-1;
      finish=0;
      while(j>=0 && !finish)
      {
        if (c>b[row][j])
        {
          if (j==N-1)
          {
            tries++;
            y+=v*w[row];
            row=i;v=1;
          }
          else
          {
            v*=(matrix[row][j+1]>=0?1:-1);
            row=j+1;
          }
          finish=1;
        }
        j--;
      }
      if (!finish)
      {
        v*=(matrix[row][j+1]>=0?1:-1);
        row=0;
      }
    }
    // calculating root
    x[i]=y/M;
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

Метод ортогонализации
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (orthogonalization method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given system of linear equations:
//      a11 x1 + a12 x2 + ... + a1n xn + a(1n+1) = 0
//      a21 x1 + a22 x2 + ... + a2n xn + a(2n+1) = 0
//
//      an1 x1 + an2 x2 + ... + ann xn + a(nn+1) = 0
//
//    Left part of each equation is a result of scalar multiplication of
//    two vectors: ai=(ai1,ai2,...,ain,a(in+1)) and x=(x1,x2,..,xn,1)
//    So to solve this system we need only to build vector x which will be
//    orthogonal to each of ai vectors.
//                            n+1
//    Let u1=a1,   b1=u1/sqrt(SUM u(1j)^2)
//                            j=1
//
//    Other ui we can get from the following iterative formula:
//                      k
//       u(k+1)=u(k+1)-SUM {u(k+1), bj}bj, where {..} means scalar multiplication
//                     j=1
//                          n+1
//       b(k+1)=u(k+1)/sqrt(SUM u(k+1,j)^2)
//                          j=1
//    And finally we can obtain roots of this system from this formula:
//
//    xi=b(n+1,i)/b(n+1n+1)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N     3     // size of matrix
#define N1    N+1
double matrix[N1][N1]=
    {{3,1,-1,-2},
  {1,4,-2,-3},
  {2,-1,5,-15}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("%+f=0\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  double c[N],x[N],s;
  int i,j,k,l,m;
  ShowMatrix();
  // expanding system by vector (0,0,0,...,0,1)
  for (i=0;i<N;i++) matrix[N][i]=0;
  matrix[N][N]=1;
  // for all equations
  for (i=0;i<N1;i++)
  {
    if (i>0)
    {
      k=0;
      // make some iterations to increase accuracy of calculations
      while (k<=3)
      {
        for (l=0;l<i;l++)
        {
          c[l]=0;
          for(m=0;m<N1;m++) c[l]+=matrix[l][m]*matrix[i][m];
        }
        for (j=0;j<N1;j++)
        {
          s=0;
          for(l=0;l<i;l++) s+=c[l]*matrix[l][j];
          matrix[i][j]-=s;
        }
        k++;
      }
    }
     // normalizing vector
     s=0;
     for (k=0;k<N1;k++) s+=matrix[i][k]*matrix[i][k];
     s=sqrt(s);
     for (j=0;j<N1;j++) matrix[i][j]/=s;
  }
  s=matrix[N][N];
  for (j=0;j<N;j++) x[j]=matrix[N][j]/s;
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

Метод Зайделя
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Zaidel method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: system of linear equations
//       a11 x1 + a12 x2 + ... + a1n xn = b1
//       a21 x1 + a22 x2 + ... + a2n xn = b2
//       ...
//       an1 x1 + an2 x2 + ... + ann xn = bn
//
//    We use following iterative formula to solve this system:
//
//                     1    i-1                N
//     xi(j+1)=xi(j)- --- ( SUM aik xk(j+1) + SUM aik xk(j) -bi)
//                    aii   k=1               k=i
//
//     where xi(j+1) is j+1-th iteration, xi(j) is j-th iteration
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N    3     // size of matrix
#define N1   N+1
float matrix[N][N1]=
    {{10,1,1,12},
  {2,10,1,13},
  {2,2,10,14}
    };
float z[N]={0,0,0};  // first iteration
float epsilon=0.0001;  // required accuracy
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float x[N];
  short int i,j;
  int iterations=0,finish=0;
  // Printing given matrix
  ShowMatrix();
  while (!finish)
  {
    finish=1;
    for (i=0;i<N;i++)
    {
      x[i]=-matrix[i][N];
      for (j=0;j<N;j++) x[i]+=matrix[i][j]*z[j];
      // don't stop iterations until required accuracy is reached
      if (fabs(x[i]/matrix[i][i])>=epsilon) finish=0;
      x[i]=z[i]-x[i]/matrix[i][i];
      // next iteration
      z[i]=x[i];
    }
    iterations++;
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
  printf("%d iterations were made",iterations);
}


Решение нелинейных уравнений

Метод хорд
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (chords method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    given: f(x)=0, [a;b], f(a)*f(b)<0
//    find x0: f(x0)=0
//    We approximate curve by chord between its borders and get
//    point where the chord crosses X axe as a new border of
//    interval with the root. When interval becomes small enough
//    we can take one of its borders as a solution.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// this function returns value of f(x)
double f(double x)
{
  // sin 2x - ln x = 0
  return sin(2*x)-log(x);
}
double Solve(double a, double b, double epsilon)
{
  double u,v,x,y;
  u=f(a); v=f(b);
  do
  {
    // x - new border
    x=(a*v-b*u)/(v-u);
    y=f(x);
    // determine whether x is left or right border
    // f(a)*f(b)<0 condition must be true
    if (y*u<0)
    {
      b=x;
      v=y;
    } else
    {
      a=x;
      u=y;
    }
  } while (b-a>epsilon);
  return x;
}
void main(void)
{
  printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
  printf("\nx0=%f",Solve(1.3,1.5,1e-9));
}

Метод дихотомии
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (dichotomia method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    given: f(x)=0, [a;b], f(a)*f(b)<0
//    find x0: f(x0)=0
//    We split interval into two parts and select that part that
//    contains root of the equation by checking condition
//    f(a)*f(b)<0 We take middle point of the interval as a new border of
//    interval with the root. When interval becomes small enough
//    we can take one of its borders as a solution.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// this function returns value of f(x)
double f(double x)
{
  // sin 2x - ln x = 0
  return sin(2*x)-log(x);
}
double Solve(double a, double b, double epsilon)
{
  double x;
  do
  {
    // x - new border
    x=(b+a)/2;
    // determine whether x is left or right border
    // f(a)*f(b)<0 condition must be true
    if (f(x)*f(a)<0) b=x;
    else a=x;
  } while (b-a>epsilon);
  return x;
}
void main(void)
{
  printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
  printf("\nx0=%f",Solve(1.3,1.5,1e-9));
}

Метод Эйткена-Стеффенсона
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (Eitken-Steffenson method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is used for solving equations like x=f(x)
//    (see also iterations method)
//    Here we use the following iterations:
//     x(n+1)=(x0*x2-x1^2)/(x0-2x1+x2)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// this function returns value of f(x)
double f(double x)
{
  // x = x + 0.37(sin 2x - ln x)
  return x+0.37*(sin(2*x)-log(x));
}
double Solve(double x, double epsilon)
{
  double x1,x2,xold,tmp;
  do
  {
    xold=x;
    x1=f(x);
    x2=f(x1);
    tmp=x-2*x1+x2;
    x=(x*x2-x1*x1)/tmp;
  } while (fabs(xold-x)>epsilon && tmp!=0);
  return x;
}
void main(void)
{
  printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
  printf("\nx0=%f",Solve(1.4,1e-6));
}

Метод итераций
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (iterations method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is used for solving equations like x=f(x)
//    We should set first approximation of the root and then
//    make some iterations x=f(x). When the difference between
//    two consequental iterations will be less than epsilon we
//    can say that we found the root.
//    To apply this method to equations like F(x)=0 we should
//    transform it into iterational form x=f(x), where f(x):
//      1) defined on [a;b]
//      2) for every point on [a;b] exists f'(x)
//      3) for every x from [a;b] f(x) is in [a;b]
//      4) exists number q : |f'(x)|<=q<1 for all x from [a;b]
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// this function returns value of f(x)
double f(double x)
{
  // x = x + 0.37(sin 2x - ln x)
  return x+0.37*(sin(2*x)-log(x));
}
double Solve(double x, double epsilon)
{
  double x1;
  do
  {
    x1=x;
    x=f(x);
  } while (fabs(x1-x)>epsilon);
  return x;
}
void main(void)
{
  printf("sin 2x - ln x = 0");
  printf("\nx0=%f",Solve(1.4,1e-6));
}

Метод Ньютона
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (Newton method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Here we use the following iterations:
//     x(n+1)=xn-f(xn)/f'(xn)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// this function returns value of f(x)
double f(double x)
{
  // sin 2x - ln x = 0
  return sin(2*x)-log(x);
}
double Solve(double x, double epsilon)
{
  double x1;
  do
  {
    x1=x;
    x=x-f(x)*epsilon/(f(x+epsilon)-f(x));
  } while (fabs(x1-x)>epsilon);
  return x;
}
void main(void)
{
  printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
  printf("\nx0=%f",Solve(1.4,1e-9));
}

Метод Лобачевского
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (Lobachevsky method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   Given: a0+a1x+a2x^2+...+anx^n=0
//   This method allows to find modulus of the greatest root of this equation
//   even if it's complex. But in last case there can appear several messages
//   about impossibilty of calculation root of negative number.
//   The main idea of this method is to change given equation to other
//   equation which roots equals to powered roots of given equation. For example
//   if roots of the given equation are x0,x1,.. xn then roots of new equation
//   will be x0^2, x1^2, ..., xn^2. Repeating this operation we get an equation
//   where one root is much greater than other ones. So we can easily
//   obtain modulus of the greatrest root of the given equation.
//   To obtain other roots of equation we need to divide given equation
//   by (x-x0) (where x0 is found root) and apply this method to result.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N   4
#define N1  N+1
#define Iterations  15  // number of iterations
double a[N1]={24,-50,35,-10,1};
void main(void)
{
  double r,b[N1],c[N1],g,bi,d;
  int z,k;
  // printing given equation
  printf("%f",a[0]);
  for(int i=1;i<N1;i++) printf("%+fx^%d",a[i],i);
  printf("=0\n\n");
  // preparing auxiliary arrays b and c
  for (i=0;i<N1;i++)
  {
    b[i]=a[i]/a[N];
    c[i]=0;
  }
  // setting required parameters
  r=1/2.0;
  g=1;
  // make all iterations
  for(int y=0;y<Iterations;y++)
  {
    // calculate coefficients c[i] (coefficients of new equation)
    z=1;
    for(i=0;i<N1;i++)
    {
      bi=z*b[i];
      k=(i+1)/2;
      for(int j=i%2;j<N1;j+=2)
      {
        c[k]+=bi*b[j];
        k++;
      }
      z=-z;
    }
    d=z*c[N-1];
    // check whether we could calculate root of d
    if(d>0)
    {
      // calculating and printing new iteration
      g*=powl(d,r);
      printf("%f\n",g);
      for (i=0;i<N1;i++)
      {
        // preparing data for next iteration
        b[i]=c[i]/powl(d,N-i);
        c[i]=0;
      }
      b[N-1]=z;
      b[N]=-z;
    } else
    {
      // d is negative - can't calculate root
      for(i=0;i<N1;i++)
      {
        // preparing data for next iteration
        b[i]=c[i];
        c[i]=0;
      }
      printf("no iteration (can't calculate root from negative number)\n");
    }
    r/=2.0;
  }
}


Интерполяция

Метод Эйткена
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
//////////////////////////////////////////////////////////////////////////////
//
//  Interpolation (using Eitken method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: several points of the unknown function
//    We should find values of this function between given points
//    Interpolation polynom isn't calculated directly:
//    To calculate values of function between given points we should
//    use the following formulas of linear interpolation:
//                  1     |y0  x0-x|
//     y(x,x0,x1)=------- |        |
//                (x1-x0) |y1  x1-x|
//
//                  1     |y0  x0-x|
//     y(x,x0,x2)=------- |        |
//                (x2-x0) |y2  x2-x|
//
//                      1   |y(x,x0,x1)  x1-x|
//     y(x,x0,x1,x2)=-------|                | and so on...
//                   (x2-x1)|y(x,x0,x2)  x2-x|
//
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#define N  4
float x[N]={-0.12, 1.68, 3.41, 5.62};
float y[N]={0.324, -0.6, -0.4, -0.21};
float xvalues[N-1],f[N-1],tmp[N];
float xa;
void main(void)
{
  // initializing variables
  for(int i=0;i<N-1;i++)
  {
    xvalues[i]=(x[i+1]+x[i])/2; // we'll calculate function between given points
    f[i]=0;
  }
  // calculating function values between given points
  for (int z=0;z<N-1;z++)
  {
    for(int j=0;j<N;j++) tmp[j]=y[j];
    for(j=0;j<N-1;j++)
      for(i=j+1;i<N;i++)
        tmp[i]=((xvalues[z]-x[j])*tmp[i]-(xvalues[z]-x[i])*tmp[j])/(x[i]-x[j]);
    f[z]=tmp[N-1];
  }
  for (i=0;i<N-1;i++)
    printf("F(%f)= %f\n",xvalues[i],f[i]);
}

Полином Лагранжа
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
//////////////////////////////////////////////////////////////////////////////
//
//  Interpolation (using Lagrange polynom)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: several points of the unknown function
//    We should find values of this function between given points and
//    approximate  function by polynom.
//    Interpolation polynom calculated by the following formula:
//             N       (x-x0)...(x-x(i-1))(x-x(i+1))...(x-xn)
//       L(x)=SUM y(i)----------------------------------------
//            i=0     (xi-x0)..(xi-x(i-1))(xi-x(i+1))..(xi-xn)
//    To calculate values of function between given points we should
//    calculate value of interpolation polynom in these points.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#define N  4
float x[N]={-0.12, 1.68, 3.41, 5.62};
float y[N]={0.324, -0.6, -0.4, -0.21};
float xvalues[N-1],f[N-1];
int i,j,z;
float l,k;
float pow[N],temp[N],tempaux[N];
void main(void)
{
  // initializing variables
  for(z=0;z<N-1;z++)
  {
    xvalues[z]=(x[z+1]+x[z])/2; // calculating function between given points
    f[z]=0;
  }
  for(z=0;z<N;z++) pow[z]=0;
  // calculating function values
  for(z=0;z<N-1;z++)
    for(i=0;i<N;i++)
    {
      l=1;
      for(j=0;j<N;j++)
        if(i!=j) l*=(xvalues[z]-x[j])/(x[i]-x[j]);
      l*=y[i];
      f[z]+=l;
    }
  // Determining interpolation polynom
  for(z=0;z<N;z++)
  {
    k=y[z];
    for(i=0;i<N;i++)
      if(i!=z) k/=(x[z]-x[i]);
    for(i=0;i<N;i++) temp[i]=0;
    temp[0]=k;
    for(i=0;i<N;i++)
      if(i!=z)
      {
        for(j=0;j<N-1;j++) tempaux[j+1]=temp[j];
        tempaux[0]=0;
        for(j=0;j<N;j++) tempaux[j]-=temp[j]*x[i];
        for(j=0;j<N;j++) temp[j]=tempaux[j];
      }
    for(i=0;i<N;i++) pow[i]+=temp[i];
  }
  // printing results
  for(z=0;z<N-1;z++) printf("F(%f)=%f\n",xvalues[z],f[z]);
  printf("\nPolynom: F(X)=");
  for(z=N-1;z>0;z--) printf("%fX^%d+\n",pow[z],z);
  printf("%f\n",pow[0]);
}

Метод Ньютона
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
//////////////////////////////////////////////////////////////////////////////
//
//  Interpolation (using Newton polynom)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: several points of the unknown function
//    We should find values of this function between given points and
//    approximate  function by polynom.
//    Interpolation polynom is searched by the following formulas:
//
// P(x)=y1+(x-x1)f(x1,x2)+(x-x1)(x-x2)f(x1,x2,x3)+...              n
//                       +(x-x1)(x-x2)...(x-x(n-1))f(x1,x2,..,xn)=SUM Ak fi
//     i            i      i-j             k  k                   i=0
// fi=MUL (x-xj) = SUM aj x    ;   ak=(-1)  MUL xj
//    j=1          j=0                       j=1
//
//         1  l      m+1          k         m
//  a(k-l)=- SUM (-1)   a(k-l+m) SUM (-1/xp); a0=1
//         l m=1                 p=1
//
//
//    To calculate values of function between given points we should
//    calculate value of interpolation polynom in these points.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N  4
float x[N]={-0.12, 1.68, 3.41, 5.62};
float Y[N]={0.324, -0.6, -0.4, -0.21};
float xvalues[N-1],f[N],a[N],y[N],pw,w,s;
int i,j,z,k,l,m,p;
void main(void)
{
  // initializing variables
  for(z=0;z<N-1;z++)
  {
    xvalues[z]=(x[z+1]+x[z])/2; // calculating function between given points
    f[z]=0;
  }
  for(z=0;z<N;z++) y[z]=Y[z];
  a[0]=1;
  f[N-1]=y[0];
  //           n-1
  // P   (x) = SUM A[k]
  //  n-1      k=0
  for(k=1;k<N;k++)
  {
    for(i=1;i<=N-k;i++) y[i-1]=(y[i]-y[i-1])/(x[i+k-1]-x[i-1]);
    //          k  k
    // a[k]=(-1)  MUL x[j]
    //            j=1
    pw=1;
    for(j=1;j<=k;j++) pw*=x[j-1];
    a[k]=(k%2==0?1:-1)*pw;
    //        1  l      m+1          k            m
    // a[k-l]=- SUM (-1)   a[k-l+m] SUM  (-1/x[p])
    //        l m=1                 p=1
    if (k!=1)
    {
      for(l=1;l<=k;l++)
      {
        w=0;
        for(m=1;m<=l;m++)
        {
          s=0;
          for(p=1;p<=k;p++) s+=pow(-1/x[p-1],m);
          w+=(m%2==0?-1:1)*a[k-l+m]*s;
        }
        a[k-l]=w/l;
      }
    }
    for(j=N;j>=N-k;j--) f[j-1]+=a[j-N+k]*y[0];
  }
  // printing results
  for(z=0;z<N-1;z++)
  {
    s=f[0];
    for(i=1;i<N;i++) s=s*xvalues[z]+f[i];
    printf("F(%f)=%f\n",xvalues[z],s);
  }
  printf("\nPolynom: F(X)=");
  for(z=0;z<N-1;z++) printf("%fX^%d+\n",f[z],N-z-1);
  printf("%f",f[N-1]);
}


Операции с полиномами

Вычисление производной
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating derivative of the polynom (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    calculating P'(x)
//    P'(x) = n a0 x^(n-1) + (n-1) a1 x^(n-2) +...
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
void PrintPolynom(float *p, int n)
{
  for(int i=0;i<n;i++)
    printf("%+f x^%d",*(p+i),n-i-1);
  printf("\n");
}
#define N   4  // polynom degree
float Polynom[N]={1,2,3,4}; // x^3 + 2x^2 + 3x + 4 = P(x)
void main(void)
{
  // printing given polynom
  PrintPolynom(Polynom,N);
  // calculating derivative
  for(int i=0;i<N;i++)
    Polynom[i]*=(N-i-1);
  // printing derivative
  PrintPolynom(Polynom,N-1);
}

Деление полиномов
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
//////////////////////////////////////////////////////////////////////////////
//
//  Dividing polynoms (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    Q(x) - polynom,  deg Q = M, M<=N
//
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    Q(x) = b0 x^m + b1 x^(m-1) + ... + bm
//
//    calculating P(x)/Q(x)=S(x)+R(x)/Q(x), deg R < deg Q
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
void PrintPolynom(float *p, int n)
{
  for(int i=0;i<n;i++)
    printf("%+f x^%d",*(p+i),n-i-1);
  printf("\n");
}
#define N    4  // deg P-1
#define M    2  // deg Q-1
float P[N]={1,1,1,-3}; // x^3 + x^2 + x - 3 = P(x)
float Q[M]={1,-1};     // x - 1 = Q(x)
float S[N-M+1];        // S(x)
                
// 0 here means constant that will be in          
// the integral I(P)=I+C
void main(void)
{
  // printing given polynoms
  printf("P(x):\n");
  PrintPolynom(P,N);
  printf("Q(x):\n");
  PrintPolynom(Q,M);
  for (int i=0;i<N+M-1;i++) S[i]=0;
  // dividing
  for(i=0;i<N-M+1;i++)
  {
    S[i]=P[i]/Q[0];
    for(int j=0;j<M;j++)
      P[i+j]-=S[i]*Q[j];
  }
  // printing result
  printf("\nresult:\n");
  PrintPolynom(S,N-M+1);
  // and remainder
  printf("remainder:\n");
  PrintPolynom(P+(N-M+1),M-1);
}

Интеграл от полинома
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
//////////////////////////////////////////////////////////////////////////////
//
//  Calculating integral of the polynom (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    calculating Integral(P(x)dx)
//    Integral(P(x)dx) =  a0/(n+1) x^(n+1) + a1/n x^n +...
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
void PrintPolynom(float *p, int n)
{
  for(int i=0;i<n;i++)
    printf("%+f x^%d",*(p+i),n-i-1);
  printf("\n");
}
#define N  4  // polynom degree
float Polynom[N+1]={8,6,4,1,0}; // 8x^3 + 6x^2 + 4x + 1 = P(x)
// 0 here means constant that will be in
// the integral I(P)=I+C
void main(void)
{
  // printing given polynom
  PrintPolynom(Polynom,N);
  // calculating integral
  for(int i=0;i<N;i++)
    Polynom[i]/=(N-i);
  // printing integral
  PrintPolynom(Polynom,N+1);
}

Умножение полиномов
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
//////////////////////////////////////////////////////////////////////////////
//
//  Multiplication of polynoms (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    Q(x) - polynom,  deg Q = M
//
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    Q(x) = b0 x^m + b1 x^(m-1) + ... + bm
//
//    calculating S(x)=P(x)*Q(x)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
void PrintPolynom(float *p, int n)
{
  for(int i=0;i<n;i++)
    printf("%+f x^%d",*(p+i),n-i-1);
  printf("\n");
}
#define N       4  // deg P-1
#define M       2  // deg Q-1
float P[N]={2,0,4,1}; // 2x^3 + 4x + 1 = P(x)
float Q[M]={1,-2};    // x - 2 = Q(x)
float S[M+N-1];       // S(x)
// 0 here means constant that will be in
// the integral I(P)=I+C
void main(void)
{
  // printing given polynoms
  PrintPolynom(P,N);
  PrintPolynom(Q,M);
  for (int i=0;i<N+M-1;i++) S[i]=0;
  // multiplying
  for(i=0;i<N;i++)
    for(int j=0;j<M;j++)
      S[i+j]+=P[i]*Q[j];
  // printing results
  PrintPolynom(S,N+M-1);
}


Сжатие и шифрование

Компрессия/декомпрессия RLE
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
//////////////////////////////////////////////////////////////////////////////
//
//  Compressing/decompressing file (RLE technology)
//  (c) Johna Smith, 1996
//
//  Method description:
//  This compression method is based on coding sequences like AAA..A
//  We can write this sequence as two bytes N A, where A is
//  character and N is amount of its repetitions. So this method is
//  very useful for *.bmp files. But if there aren't such sequences
//  then each separate character will turn into two - 1 A 1 B ...
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <string.h>
#include <errno.h>
void usage(void)
{
  // show how to use this program if there isn't enough parameters
  printf("Usage: compress.exe filename1 COMPRESS|DECOMPRESS filename2\n");
  printf("filename1 - name of source file\n");
  printf("filename2 - name of destination file\n");
  printf("\nExample: compress.exe compress.bmp COMPRESS compress.rle\n");
}
FILE *input, *output;
void compress(void)
{
  unsigned char c1,c,length;
  unsigned long int rcounter=0,wcounter=0;
  fread(&c1,sizeof(char),1,input);
  do
  {
    c=c1;
    length=0;
    do
    {
      fread(&c1,sizeof(char),1,input);
      length++;
      rcounter++;
    } while (!feof(input) && c1==c && length!=255);
    // putting char and number of its repetitions to output
    fwrite(&c,sizeof(char),1,output);
    fwrite(&length,sizeof(char),1,output);
    wcounter++;
  } while (!feof(input));
  printf("%ld byte(s) compressed\n",rcounter);
  printf("ratio: %d\n",wcounter/rcounter);
}
void decompress(void)
{
  unsigned char c,length;
  unsigned long int rcounter=0,wcounter=0;
  fread(&c,sizeof(char),1,input);
  do
  {
    fread(&length,sizeof(char),1,input);
    for(int i=0;i<length;i++)
    {
      fwrite(&c,sizeof(char),1,output);
      wcounter++;
    }
    fread(&c,sizeof(char),1,input);
    rcounter+=2;
  } while (!feof(input));
  printf("%ld byte(s) decompressed\n",wcounter);
  printf("ratio: %d\n",wcounter/rcounter);
}
void main(int argc, char **argv)
{
  // copyleft;)_
  printf("Compress. copyleft {c} 1996 Johna Smith. freeware\n\n");
  // analysing parameters
  if (argc!=4)
  {
    usage();
    return;
  }
  if ((input=fopen(argv[1],"rb"))==NULL)
  {
    // error opening file
    printf("! ERR03: can't open source file :(\n");
    return;
  };
  if ((output=fopen(argv[3],"wb"))==NULL)
  {
    // error opening file
    fclose(input);
    printf("! ERR04: can't open destination file :(\n");
    return;
  };
  if (stricmp(argv[2],"COMPRESS")==0) compress();else
  if (stricmp(argv[2],"DECOMPRESS")==0) decompress();else
  {
    // must be COMPRESS or DECOMPRESS magic word there to know what to do
    printf("! ERR02: serious bug was found in second parameter :)\n");
    fclose(input);
    fclose(output);
    return;
  }
  fclose(input);
  fclose(output);
}

Шифрование/дешифрование сдвигом
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
//////////////////////////////////////////////////////////////////////////////
//
//  Encrypting/decrypting file (shift method)
//
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Set the base for randomizer (this is the main code for decrypting)
//    2) Get random number using the following formula: x=fract(11*x+Pi)
//    3) Crypt character using formula ch'=ch+rnd()*255
//    4) Repeat steps 3 & 4 for all chars in the text
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <errno.h>
#define Pi  3.14159265358979323846
void usage(void)
{
  // show how to use this program if there isn't enough parameters
  printf("Usage: crypt.exe filename1 DECRYPT|ENCRYPT filename2 key\n");
  printf("filename1 - name of source file\n");
  printf("filename2 - name of destination file\n");
  printf("key - key code, must be a number\n");
  printf("\nExample: crypt.exe crypt.cpp ENCRYPT crypt.crp 123.45678\n");
}
float rnd(float _base=0)  // generates random number
{
  static float base;
  float tmp;
  // reinitializing randomizer if _base!=0
  if (_base!=0) base=_base;
  // generating random number
  tmp=11*base+Pi;
  base=tmp-(int)tmp;
  return base;
}
FILE *input, *output;
void encrypt(void)
{
  unsigned char ch;
  int tmp;
  unsigned long int counter=0;
  fread(&ch,sizeof(char),1,input);
  do
  {
    // encrypting char
    tmp=rnd()*255;
    if (ch+tmp>255) ch+=(tmp-255);
    else ch+=tmp;
    // putting char to output
    fwrite(&ch,sizeof(char),1,output);
    counter++;
    // getting char from input
    fread(&ch,sizeof(char),1,input);
  } while (!feof(input));
  printf("%ld byte(s) encrypted\n");
}
void decrypt(void)
{
  unsigned char ch;
  int tmp;
  unsigned long int counter=0;
  fread(&ch,sizeof(char),1,input);
  do
  {
    // decrypting char
    tmp=rnd()*255;
    if ((int)ch<tmp) ch-=(tmp-255);
    else ch-=tmp;
    // putting char to output
    fwrite(&ch,sizeof(char),1,output);
    counter++;
    // getting char from input
    fread(&ch,sizeof(char),1,input);
  } while (!feof(input));
  printf("%ld byte(s) decrypted\n");
}
void main(int argc, char **argv)
{
  // copyleft;)_
  printf("Crypt. copyleft {c} 1996 Johna Smith. freeware\n\n");
  // analysing parameters
  if (argc!=5)
  {
    usage();
    return;
  }
  rnd(atof(argv[4])); // setting up randomizer base
  if (errno==ERANGE)
  {
    // wrong key code
    printf("! ERR01: key code is out of range\n");
    return;
  }
  if ((input=fopen(argv[1],"rb"))==NULL)
  {
    // error opening file
    printf("! ERR03: can't open source file :(\n");
    return;
  };
  if ((output=fopen(argv[3],"wb"))==NULL)
  {
    // error opening file
    fclose(input);
    printf("! ERR04: can't open destination file :(\n");
    return;
  };
  if (stricmp(argv[2],"ENCRYPT")==0) encrypt();else
  if (stricmp(argv[2],"DECRYPT")==0) decrypt();else
  {
    // must be ENCRYPT or DECRYPT magic word there to know what to do
    printf("! ERR02: serious bug was found in second parameter :)\n");
    fclose(input);
    fclose(output);
    return;
  }
  fclose(input);
  fclose(output);
}

Шифрование/дешифрование ХОR-ом
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
//////////////////////////////////////////////////////////////////////////////
//
//  Encrypting/decrypting file (xor method)
//
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Get the key code (this is the main code for crypting)
//    2) Crypt part of text using formula str'=str xor keycode
//    3) Repeat step 2 for all parts of text
//
//   (a xor b) xor b = a -> encrypting==decrypting
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <errno.h>
void usage(void)
{
  // show how to use this program if there isn't enough parameters
  printf("Usage: xorcrypt.exe filename1 filename2 key\n");
  printf("filename1 - name of source file\n");
  printf("filename2 - name of destination file\n");
  printf("key - key code, can be up to 20 characters (without spaces)\n");
  printf("\nExample: xorcrypt.exe xorcrypt.cpp xorcrypt.crp SAMPLEKEYCODE\n");
}
FILE *input, *output;
char keycode[20];
char keycodelength;
void encrypt(void)
{
  unsigned char ch;
  char j=0;
  unsigned long int counter=0;
  fread(&ch,sizeof(char),1,input);
  do
  {
    // processing char
    ch=(ch | keycode[j])-(ch & keycode[j]);  // a xor b = a or b-a and b
    j++;
    if (j==keycodelength) j=0;
    // putting char to output
    fwrite(&ch,sizeof(char),1,output);
    counter++;
    // getting char from input
    fread(&ch,sizeof(char),1,input);
  } while (!feof(input));
  printf("%ld byte(s) processed\n");
}
void main(int argc, char **argv)
{
  // copyleft;)_
  printf("XORCrypt. copyleft {c} 1996 Johna Smith. freeware\n\n");
  // analysing parameters
  if (argc!=4)
  {
    usage();
    return;
  }
  if (strlen(argv[3])>20)
  {
    // key code is too long
    printf("! ERR01: key code is too long\n");
    return;
  }
  strcpy(keycode,argv[3]);
  keycodelength=strlen(keycode);
  if ((input=fopen(argv[1],"rb"))==NULL)
  {
    // error opening file
    printf("! ERR03: can't open source file :(\n");
    return;
  };
  if ((output=fopen(argv[2],"wb"))==NULL)
  {
    // error opening file
    fclose(input);
    printf("! ERR04: can't open destination file :(\n");
    return;
  };
  encrypt();
  fclose(input);
  fclose(output);
}


Решение дифференциальных уравнений

Метод Эйлера
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Eiler method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y=h*f(x,y(i))
//
//    In this example y'=cos(y)+3x, [a;b]=[0;1], x0=0, y0=1.3
//
//////////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <stdio.h>
const float a=0,b=1; // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=1.3; // initial conditions
int M=1;
float f(float x, float y)
{
  return (cos(y)+3*x); // y'=cos(y)+3*x
}
void calculate(int m,float *y)
{
  float x,yi,h;
  h=(b-a)/(num_points*m);
  yi=y0; x=x0;
  for (int i=0;i<num_points;i++)
  {
    for (int k=0;k<m;k++)
    {
      yi+=h*f(x,yi);
      x+=h;
    }
    *(y+i)=yi;
  }
}
void main(void)
{
  float yh[num_points],yh2[num_points];
  calculate(M,yh);
  calculate(2*M,yh2);  // doubled step for better accuracy
  // epsilon is difference between solutions with single
  // and double steps
  printf("X\t\tYH\t\tYH2\t\tEPSILON\n");
  for (int i=0;i<num_points;i++)
    printf("%f\t%f\t%f\t%f\n",(x0+((i+1)*(b-a))/num_points),
                                      yh[i],yh2[i],yh2[i]-yh[i]);
}

Метод Адамса
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Adams method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    We use Runge-Kutta method to process first 'num_starting_points'
//    points and then use Adams extrapolation to process other points.
//    In this example y'=x+y, [a;b]=[0;2], x0=0, y0=1
//
//////////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <stdio.h>
const float a=0,b=2;             // bounds of the interval
const int num_points=10,         // number of points to solve
          num_starting_points=4; // number of points to solve with Runge-Kutta method
float x0=0,y0=1;                 // starting conditions
float f(float x, float y)
{
  return x+y;  // y'=x+y
}
// this function realises Runge-Kutta method for n starting points
void calculate(float *y)
{
  float k1,k2,k3,k4,x,yi,h;
  h=(b-a)/num_points;  // step
  yi=y0; x=x0;
  for (int i=0;i<num_starting_points;i++)
  {
    k1=h*f(x,yi);
    k2=h*f(x+h/2,yi+k1/2);
    k3=h*f(x+h/2,yi+k2/2);
    k4=h*f(x+h,yi+k3);
    yi+=(k1+2*k2+2*k3+k4)/6;
    x+=h;
    *(y+i+1)=yi;
  }
}
void main(void)
{
  float y[num_points+1],h;
  y[0]=y0;
  // apply Runge-Kutta method
  calculate(y);
  h=(b-a)/num_points;
  // extrapolating
  for (int i=num_starting_points;i<num_points;i++)
    y[i] = y[i-1]+h/24*(55*f(x0+(i-1)*h,y[i-1])-
             59*f(x0+(i-2)*h,y[i-2])+
             37*f(x0+(i-3)*h,y[i-3])-
             9*f(x0+(i-4)*h,y[i-4]));
  printf("X\t\tY\t\tExact solution\n");
  for (i=0;i<num_points;i++)
    printf("%f\t%f\t%f\n",(x0+i*h),y[i],(2*exp(x0+i*h)-(x0+i*h)-1));
}

Метод Эйлера-Коши
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Eiler-Coshi method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y=h*(f(x,y(i))+f(x,z))/2, z=y(i)+h*(f(x(i),y(i))
//
//    In this example y'=cos(y)+3x, [a;b]=[0;1], x0=0, y0=1.3
//
//////////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <stdio.h>
const float a=0,b=1;     // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=1.3;       // initial conditions
int M=1;
float f(float x, float y)
{
  return (cos(y)+3*x); // y'=cos(y)+3*x
}
void calculate(int m,float *y)
{
  float x,yi,h,z;
  h=(b-a)/(num_points*m);
  yi=y0; x=x0;
  for (int i=0;i<num_points;i++)
  {
    for (int k=0;k<m;k++)
    {
      z=yi+h*f(x,yi);
      yi+=h*(f(x,yi)+f(x,z))/2;
      x+=h;
    }
    *(y+i)=yi;
  }
}
void main(void)
{
  float yh[num_points],yh2[num_points];
  calculate(M,yh);
  calculate(2*M,yh2);  // doubled step for better accuracy
  // epsilon is difference between solutions with single
  // and double steps
  printf("X\t\tYH\t\tYH2\t\tEPSILON\n");
  for (int i=0;i<num_points;i++)
    printf("%f\t%f\t%f\t%f\n",(x0+((i+1)*(b-a))/num_points),
                                       yh[i],yh2[i],yh2[i]-yh[i]);
}

Метод Рунге-Кутта
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Runge-Kutta method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y =1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4
//    k1=hf(x,y)
//    k2=hf(x+h/2,y+k1/2)
//    k3=hf(x+h/2,y+k2/2)
//    k4=hf(x+h,y+k3)
//
//    In this example y'=cos(x+y)+x-y, [a;b]=[0;1], x0=0, y0=0
//
//////////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <stdio.h>
const float a=0,b=1;     // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=0;         // initial conditions
int M=1;
float f(float x, float y)
{
  return (cos(x+y)+x-y); // y'=cos(x+y)+x-y
}
void calculate(int m,float *y)
{
  float k1,k2,k3,k4,x,yi,h;
  h=(b-a)/(num_points*m);
  yi=y0; x=x0;
  for (int i=0;i<num_points;i++)
  {
    for (int k=0;k<m;k++)
    {
      // calculatin coefficients k
      k1=h*f(x,yi);
      k2=h*f(x+h/2,yi+k1/2);
      k3=h*f(x+h/2,yi+k2/2);
      k4=h*f(x+h,yi+k3);
      yi+=(k1+2*k2+2*k3+k4)/6;
      x+=h;
   }
   *(y+i)=yi;
  }
}
void main(void)
{
  float yh[num_points],yh2[num_points];
  calculate(M,yh);
  calculate(2*M,yh2);  // doubled step for better accuracy
  // epsilon is difference between solutions with single
  // and double steps
  printf("X\t\tYH\t\tYH2\t\tEPSILON\n");
  for (int i=0;i<num_points;i++)
    printf("%f\t%f\t%f\t%f\n",(x0+((i+1)*(b-a))/num_points),
                                    yh[i],yh2[i],(yh2[i]-yh[i]));
}


Вычисление CRC

CRC-16 DOS
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
//////////////////////////////////////////////////////////////////////////////
//
//  Calculacting CRC-16 (DOS) of data block
//  (c) Johna Smith, 1996
//
//  Method description:
//    CRC is used to check data sequences for errors
//    CRC-16 DOS is DOS modification of 16-bit control code (CRC)
//    To calculate CRC-16 DOS
//    1) CRC=0
//    2) Call CRC16DOS for each byte of data block
//    3) Call CRC16DOS(0) twice
//    Function CRC16DOS recieves current CRC value and returns updated value
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
int CRC16DOS(char c, unsigned int crc)
{
  unsigned int CRC_MASK=0xa001;
  asm
  {
    mov     al,c      // register AL contains current byte value (8 bits)
    mov     dx,crc    // register DX contains current CRC value (16 bits)
    push    cx        // store CX register in stack
    mov     cx,8      // put number 8 into CX register (CX=8)
  CRC_Loop:
    ror     al,1
    //  ROR  AL,1 - shifting bits of AL to the right, putting lowest
    //  bit to highest bit and to carry flag
    //  Example: 10010110b -> 01001011b and Carry Flag is set to 0
    rcr     dx,1
    //  RCR  DX,1  -  shifting  bits  of  DX to the right (like >> operation)
    //  and putting value from Carry Flag to the highetst bit
    //  Example: 10010110b (Carry Flag=1) -> 11001011 (Carry Flag=0)
    jnc     CRC2  // Goto CRC2 if Carry Flag is 0
    xor     dx,CRC_MASK // DX=DX xor CRC_MASK=(DX | CRC_MASK) - (DX & CRC_MASK)
  CRC2:
    loop    CRC_Loop // CX--; if (CX!=0) goto CRC_Loop
    pop     cx       // restore CX register from stack
    mov     crc,dx   // CRC=DX
  }
  return crc;
}
void main(void)
{
  int crc=0;   // initially CRC=0
  char data[]="This is data sample.";
  // Call CRC16DOS for each byte in the block
  for (int i=0;i<sizeof(data)/sizeof(char);i++) crc=CRC16DOS(data[i],crc);
  // Calculate CRC(0) twice
  crc=CRC16DOS(0,crc);
  crc=CRC16DOS(0,crc);
  // Print calculated value
  printf("%u\n",crc);
}

CRC-16 CCITT
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
//////////////////////////////////////////////////////////////////////////////
//
//  Calculacting CRC-16 (CCITT) of data block
//  (c) Johna Smith, 1996
//
//  Method description:
//    CRC is used to check data sequences for errors
//    CRC-16 CCITT is CCITT modification of 16-bit control code (CRC)
//    To calculate CRC-16 CCITT
//    1) CRC=0
//    2) Call CRC16CCITT for each byte of data block
//    3) Call CRC16CCITT(0) twice
//    Function CRC16CCITT recieves current CRC value and returns updated value
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
int CRC16CCITT(char c, unsigned int crc)
{
  unsigned int CRC_MASK=0x1021;
  asm
  {
    mov     al,c     // register AL contains current byte value (8 bits)
    mov     dx,crc   // register DX contains current CRC value (16 bits)
    push    cx       // store CX register in stack
    mov     cx,8     // put number 8 into CX register (CX=8)
  CRC_Loop:
    rol     al,1
    //  ROL  AL,1 - shifting bits of AL to the left, putting highest
    //  bit to the lowest bit and to carry flag
    //  Example: 11010110b -> 10101101b and Carry Flag is set to 1
    rcl     dx,1
    //  RCL  DX,1  -  shifting  bits  of  DX to the left (like << operation)
    //  and putting value from Carry Flag to the lowest bit
    //  Example: 01010110b (Carry Flag=1) -> 10101101b (Carry Flag=0)
    jnc     CRC2     // Goto CRC2 if Carry Flag is 0
    xor     dx,CRC_MASK // DX=DX xor CRC_MASK=(DX | CRC_MASK) - (DX & CRC_MASK)
  CRC2:
    loop    CRC_Loop // CX--; if (CX!=0) goto CRC_Loop
    pop     cx       // restore CX register from stack
    mov     crc,dx   // CRC=DX
  }
  return crc;
}
void main(void)
{
  int crc=0;   // initially CRC=0
  char data[]="This is data sample.";
  for (int i=0;i<sizeof(data)/sizeof(char);i++) crc=CRC16CCITT(data[i],crc);
  // Calculate CRC(0) twice
  crc=CRC16CCITT(0,crc);
  crc=CRC16CCITT(0,crc);
  // Print calculated value
  printf("%u\n",crc);
}

CRC-32
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
//////////////////////////////////////////////////////////////////////////////
//
//  Calculacting CRC-32 of data block
//  (c) Johna Smith, 1996
//
//  Method description:
//    CRC is used to check data sequences for errors
//    CRC-32 is 32-bit control code (CRC)
//    To calculate CRC-32:
//    1) CRC=-1
//    2) Call CRC32 for each byte of data block
//    3) CRC=~CRC
//    Function CRC32 recieves current CRC value and returns updated value
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
long int CRC32(char c, long int crc)
{
  int CRC_MASK=0x8320; // CRC-32 mask lo word
  int CRC_MASK1=0xEDB8; // CRC-32 mask hi word
  char CRC_FF=CRC_MASK & 0xFF;
  int d,b;
  d=(int)crc;  // d is lo word of CRC-32
  b=crc>>16;   // b is hi word of CRC-32
  asm
  {
    mov al,c   // register AL contains current byte value (8 bits)
    mov bx,b   // register BX contains hi word of current CRC value (16 bits)
    mov dx,d   // register BX contains lo word of current CRC value (16 bits)
    push cx    // store CX register in stack
    push ax    // store AX register in stack
    xor al,dl  // AL=AL xor DL=(AL | DL) - (AL & DL)
    mov cx,8   // put number 8 into CX register (CX=8)
  CRC_Loop:
    shr bx,1   // bx>>=1 (shifting bits of BX to the right)
    rcr dx,1
    //  RCR  DX,1  -  shifting  bits  of  DX to the right (like >> operation)
    //  and putting value from Carry Flag to the highetst bit
    //  Example: 10010110b (Carry Flag=1) -> 11001011 (Carry Flag=0)
    shr al,1   // al>>=1 (shifting bits of AL to the right)
    jnc CRC2   // Goto CRC2 if Carry Flag is 0
    xor dx,CRC_MASK  // DX=DX xor CRC_MASK=(DX | CRC_MASK) - (DX & CRC_MASK)
    xor bx,CRC_MASK1 // BX=BX xor CRC_MASK1=(BX | CRC_MASK1) - (BX & CRC_MASK1)
    xor al,CRC_FF // AL=AL xor CRC_MFF=(AL | CRC_FF) - (AL & CRC_FF)
  CRC2:
    loop CRC_Loop // CX--; if (CX!=0) goto CRC_Loop
    pop ax        // restore AX register from stack
    pop cx        // restore CX register from stack
    mov b,bx      // b=BX
    mov d,dx      // d=DX
  }
  // combining lo word and hi word to unsigned long crc
  crc=( (unsigned long)( ((unsigned)(d)) | (((unsigned long)((unsigned)(b)))<<16)) );
  return crc;
}
void main(void)
{
  unsigned long int crc32=-1; // initially CRC=-1
  char data32[]={1,2,3,4,5,6,7,8,9,0};
  // Call CRC32 for each byte of data block
  for (int i=0;i<sizeof(data32)/sizeof(char);i++) crc32=CRC32(data32[i],crc32);
  crc32=~crc32;
  // Print calculated value
  printf("%lu",crc32);
}
 
Текущее время: 06:48. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Рейтинг@Mail.ru