Аппроксимация тригонометрическим полиномом

//apprrrrrrrrrrrrr///////////////////////////////
 
double *apr_b, *apr_x, *apr_y, **apr_sums;
int apr_N;
//N - кол-во шагов
//K - степень
 
void allocmatrix(){
 
int i,j;
apr_a = new double[apr_K+1];
apr_b = new double[apr_K+1];
apr_sums = new double*[apr_K+1];
 
for(i=0; i<apr_K+1; i++){
apr_sums[i] = new double[apr_K+1];
 
}
for(i=0; i<apr_K+1; i++){
apr_a[i]=0;
apr_b[i]=0;
for(j=0; j<apr_K+1; j++){
apr_sums[i][j] = 0;
}
}
 
}
 
 
 
void readmatrix(){
int i=0,j=0, k=0;
 
for(i=0; i<apr_K+1; i++){
for(j=0; j<apr_K+1; j++){
apr_sums[i][j] = 0;
for(k=0; k<apr_N; k++){
apr_sums[i][j] += pow(apr_x[k], i+j); то есть эту строку поменять
}
}
}
for(i=0; i<apr_K+1; i++){
for(k=0; k<apr_N; k++){
apr_b[i] += pow(apr_x[k], i) * apr_y[k]; эту тоже
}
}
}
 
 
void diagonal(){
int i, j, k;
float temp=0;
for(i=0; i<apr_K+1; i++){
if(apr_sums[i][i]==0){
for(j=0; j<apr_K+1; j++){
if(j==i) continue;
if(apr_sums[j][i] !=0 && apr_sums[i][j]!=0){
for(k=0; k<apr_K+1; k++){
temp = apr_sums[j][k];
apr_sums[j][k] = apr_sums[i][k];
apr_sums[i][k] = temp;
}
temp = apr_b[j];
apr_b[j] = apr_b[i];
apr_b[i] = temp;
break;
}
}
}
}
}
 
 
void approxn1(){
int i=0,j=0, k=0;
allocmatrix();
//read data from file
readmatrix();
//check if there are 0 on main diagonal and exchange rows in that case
diagonal();
 
for(k=0; k<apr_K+1; k++){
for(i=k+1; i<apr_K+1; i++){
double M = apr_sums[i][k] / apr_sums[k][k];
for(j=k; j<apr_K+1; j++){
apr_sums[i][j] -= M * apr_sums[k][j];
}
apr_b[i] -= M*apr_b[k];
}
}
 
for(i=(apr_K+1)-1; i>=0; i--){
double s = 0;
for(j = i; j<apr_K+1; j++){
s = s + apr_sums[i][j]*apr_a[j];
}
apr_a[i] = (apr_b[i] - s) / apr_sums[i][i]; 
}
 
}
 
//end apprrrr/////////////////////////////////

apr_x, apr_y - координаты по х и y.
NSteps - кол-во шагов;
double l=(10+10)/2 период;
pa-степень аппроксимации.
const double pi = 3.1415926535897932384626433832795;
double a[12],b[12];
int ia,ja,ka,pa=6;
for (ka=1;ka<pa;ka++)
{
    a[ka]=0;b[ka]=0;
    for(ia=0;ia<NSteps;ia++)
    {   
        a[ka]=a[ka]+apr_y[ia]*cos(2*ka*pi*ia/(NSteps-1));
        b[ka]=b[ka]+apr_y[ia]*sin(2*ka*pi*ia/(NSteps-1));
    }
    a[ka]=2*a[ka]/(NSteps+1);
    b[ka]=2*b[ka]/(NSteps+1);
}
a[0]=0;
for(ia=0;ia<NSteps;ia++)
{
    a[0]=a[0]+apr_y[ia];
}
a[0]=a[0]/(NSteps+1);
 
for (ii = 0; ii< NSteps; ii++)
{
    apr_y_res[ii] = 0; 
    
    apr_y_res[ii] =a[0]/2;
    for (int jj2=1;jj2<pa;jj2++)
        {
            apr_y_res[ii] = apr_y_res[ii] + a[jj2]*cos(jj2*apr_x[ii])+b[jj2]*sin(jj2*apr_x[ii]);
        }
    
}