Коллекция алгоритмов от Johna Smith - C++ - Ответ 4579895
19.05.2013, 15:58. Показов 4289. Ответов 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
| //////////////////////////////////////////////////////////////////////////////
//
// 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);
} |
|
Вернуться к обсуждению: Коллекция алгоритмов от Johna Smith C++
1
|