Точность расчетов

struct factor
{
	double* arg;   // коэффициент при комплексном аргументе s ;
	int max_degree;  // старшая степень комплексного аргумента s;
	double denominator; // модуль корня знаменателя простейшей дроби;
	int denominator_degree; // степень знаменателя простейшей дроби;
};
void Throughput(int, unsigned long**);
float time_delay(VPath* route, int L, float* CLoad);
float integral(int degree_arg, float arg, float intensity);
float factorial(int n);
template <class Type> Type** zeroing(int num, Type**);
int ElementNum (int i, int j, int size);

int Hna = 48;       // длина заголовка уровня сетевого доступа (канальный уровень);
int Hip = 320;      // длина IP заголовка;
unsigned long** V;
int num_knots = 4;
float Td = 0.5;     

int _tmain(int argc, _TCHAR* argv[])
{
	V = zeroing(num_knots,V);
	Throughput(num_knots, V);
	float* CLoad = new float[num_knots*num_knots-num_knots]; // матрица загрузок каналов;
	VPath* route = new VPath;
	route->knot = new int[3];
	route->knot[0] = 4;
	route->knot[1] = 1;
	route->knot[2] = 2;
	route->knot[3] = 3;
	route->lenght = 3;
	CLoad[ElementNum(3,0,num_knots)] = 0.80013329;
	CLoad[ElementNum(0,1,num_knots)] = 0.80013335;
	CLoad[ElementNum(1,2,num_knots)] = 0.80013339;
	cout << time_delay(route,8000,CLoad);
	float a11 = 6.3957348 ;
	float a21 = 6.3957329 ;

	cin.get();
	return 0;
}
float time_delay(VPath* route, int L, float* CLoad) // CLoad (Channel Loading) загрузка канала, L - длина пакета, Ts - время задержки, V - матрица пропускных способностей ЛЦТ-ов, Route - набор виртуальных путей;
{
	int Vb = 64000;                               // скорость вокодера;
	float Tq = Td - (float)(L - Hip) / (float)Vb; // заданное время пребывания пакета в пути (исключая время пакетизации);
	double intensity;                              // intensity - интенсивность потока;
	float result = 0;
	double PTime;
	double summA;
	factor* A;
	do{
		summA = 0;
		PTime = 0;
		A = new factor[route->lenght + 1];   // вектор коэффициентов;
		for(int num_A = 0; num_A< route->lenght + 1; num_A++){
			A[num_A].arg = new double[route->lenght];
			A[num_A].max_degree = 0;
			A[num_A].denominator = 0;
			A[num_A].denominator_degree = 1;
			for(int j = 0; j < route->lenght; j++){
				if(j != 0) A[num_A].arg[j] = 0;
				else A[num_A].arg[j] = 1;
			}
		}

		for(int i=0; i < route->lenght; i++){            // i - номер вершины в пути
			intensity =( V[route->knot[i]-1][route->knot[i+1]-1] / (L + Hna) * (1 - CLoad[ElementNum (route->knot[i]-1, route->knot[i+1]-1, num_knots)]));
			for(int num_A = 0; num_A < route->lenght; num_A++){         // num_A - номер коэффициента
				if(num_A != i){					// проверка умножения №1
					if(A[num_A].denominator != intensity){  // проверка умножения №2
						for(int k = A[num_A].max_degree + 1; k>=0; k--){          // k - степень аргумента S
							if(k == 0)
								A[num_A].arg[k] = A[num_A].arg[k] * intensity;
							else{
								if(k == A[num_A].max_degree)
									A[num_A].arg[k] = A[num_A].arg[k-1] + intensity;
								else{
									if(k == A[num_A].max_degree + 1)
										A[num_A].arg[k] = 1;
									else
										A[num_A].arg[k] = A[num_A].arg[k] * intensity + A[num_A].arg[k-1];
								}
							}
						}
						A[num_A].max_degree++;
					}
					else
						A[num_A].denominator_degree++;
				}
				else
					A[num_A].denominator = intensity;
			}
			A[route->lenght].arg[0] = A[route->lenght].arg[0] * intensity;
		}
		[B]for(int num_A = 0; num_A < route->lenght; num_A++){
			for(int s = 0; s < route->lenght; s++)
				if(s != num_A){
					if (A[num_A].arg[num_A] == 0){
						for (int i = 0; i < route->lenght + 1; i++)
							A[i].arg[num_A] = A[i].arg[num_A] + A[i].arg[num_A+1];
					}
					intensity = A[num_A].arg[s];
					for(int i = 0; i < route->lenght + 1; i++){
						if(i == num_A)
							A[i].arg[s] = 0;
						else
							A[i].arg[s] = A[i].arg[s] - intensity * (A[i].arg[num_A] / A[num_A].arg[num_A]);
					}
				}
		}
[/B]
		
		[B]for(int num_A = 0; num_A < route->lenght; num_A++){
			if(num_A == route->lenght-1)
				A[num_A].arg[num_A] = -summA;
			else{
				A[num_A].arg[num_A] = A[route->lenght].arg[num_A]/A[num_A].arg[num_A];
				summA = summA + A[num_A].arg[num_A];
			}[/B]
			
			if(A[num_A].arg[num_A] != 0){
				intensity = ( V[route->knot[num_A]-1][route->knot[num_A+1]-1] / (float)(L + Hna) * (1 - CLoad[ElementNum (route->knot[num_A]-1, route->knot[num_A+1]-1, num_knots)]));
				PTime = PTime + (A[num_A].arg[num_A] / factorial(A[num_A].denominator_degree - 1)) * integral(A[num_A].denominator_degree - 1,Tq,intensity);
			}
		}
		result = result + (1-PTime);
		//cout << endl;
		//cout << result << endl;  // проверка
		//cout << endl;
		for(int num_A = 0; num_A< route->lenght + 1; num_A++)
			delete [] A[num_A].arg;
		delete[] A;
		A = (factor*)NULL;
		route = route->next_vpath;
	}while(route != (VPath*)NULL);
	return result;
}
float integral(int degree_arg, float arg, float intensity)   // рекурсивная функция поиска интеграла
{
	float result;
	if(degree_arg != 0)
		result = -(pow(arg,degree_arg) * exp(-intensity * arg))/intensity + (degree_arg/intensity)*integral(degree_arg - 1,arg,intensity);
	else
		result = (1-exp(-intensity * arg))/intensity;
	return result;
}
float factorial(int n)
{
	float result = 1;
	if(n != 0)
		for(int i = n; i>0; i--)
			result = result * i;
	return result;
}
void Throughput(int num_knots, unsigned long** result)
{
	for(int i = 0; i < num_knots; i++)
		for(int j = i; j < num_knots; j++){
			if (i != j){
				// cout << "VVedite V["<<i+1<<"]["<<j+1<<"]: ";
				// cin >> result[i][j];
				result[i][j] = 256 * 1024*(i+1);           // заполнение матрицы;
				result[j][i] = 256 * 1024*(i+1);
			}
		}
}
template <class Type> Type** zeroing(int num, Type** result)
{
	result = new Type*[num];
	for(int i = 0; i < num; i++){
		result[i] = new Type[num];
		for(int j = 0; j < num; j++)
			result[i][j] = 0;
	}
	return result;
}
int ElementNum (int i, int j, int size)
{
	if (j > i)
		return i*(size-1) + j - 1;
	else
		return i*(size-1) + j;
}