#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <immintrin.h> //int _rdrand64_step
#include <iostream> //std::cin.get()
const int matrixNumber = 4131;//3 oder 4131
double cputime()
{
	return ((double)clock() * 1000.0/ CLOCKS_PER_SEC);
}
 
int main()
{
	double* m[matrixNumber];
        double merk[matrixNumber];//gemerkte letzte Zeile
	for (int i = 0; i < matrixNumber; i++)
		m[i]=(double*)malloc((matrixNumber + 1) * sizeof(double));
	double dF,dD;
	int IR, IC,NR = matrixNumber;//Zeilen  len(m)
	int	NC = matrixNumber+1;     //Spalten
	unsigned __int64 val;
	int iOK = _rdrand64_step(&val);
	for (int IT = 0; IT < NR; IT++)
	{
		for (IR = 0; IR < NC; IR++) 
		{
			iOK = _rdrand64_step(&val); m[IT][IR] = ((double)val)*1e-19;
		} 
	}
        for (IR = 0; IR < matrixNumber; IR++)  merk[IR] = m[matrixNumber - 1][IR];
	printf("die ersten 5 Ergebnisse einer Matrix mit %ld Zeilen:\n", matrixNumber);
	for (int i = 0; i < 5; i++) printf(" %6.15f ", m[i][matrixNumber]);
        printf("\nsoll= %6.15f \n", m[matrixNumber - 1][matrixNumber]);
	printf("nun berechnen...\n");
	double begin = cputime();
	for (int IT = 0; IT < NR; IT++) //IT in range(0, NR)
	{
		for (IR = 0; IR < IT; IR++) 
		{
			dF = m[IR][IT] / m[IT][IT];
			for ( IC = 0; IC < NC; IC++) 
				m[IR][IC] = m[IR][IC] - dF * m[IT][IC];
		}
		for (IR = 1 + IT; IR < NR; IR++)
		{
			dF = m[IR][IT] / m[IT][IT];
			for (IC = 0; IC < NC; IC++)
				m[IR][IC] = m[IR][IC] - dF * m[IT][IC];
		}
	}
	for (IR = 0;IR < NR;IR++)
	{
		dD = m[IR][IR];
		for (IC = 0; IC < NC;IC++)
				m[IR][IC] = m[IR][IC] / dD;
	}
	begin = (cputime() - begin) / 1000.0;	printf("\nin %6.3f s\n", begin);
	for (int i = 0; i < 5; i++) printf(" %6.15f ", m[i][matrixNumber]);
	double dSum = 0.0;
	for (int i = 0; i < matrixNumber; i++) dSum += merk[i] * m[i][matrixNumber];
	printf("ist = %6.15f \n", dSum);
	for (int i = 0; i < matrixNumber; i++) 	free(m[i]);
	std::cin.get();//mit Compiler AVX2: 44s
}