#include<stdio.h>
#include<stdlib.h>
long int getoff(int index,const int * dims,const int * stride, int n) 
{ 	
	long int ret = 0; 
	for(int i = 0; i < n; i++)
	{
		int ii = index % dims[i];
		ret += ii * stride[i];
		index/= dims[i];
	} 	return ret;
} 

void transpose(int n, const double *A, double *B,  const int *lda,  int *perm)
{
	int i = 0, j;
	int r_perm[100], lda_s[100], ldb_s[100],temp[100] ;
	lda_s[0] = 1;
	ldb_s[0] = 1;
	long size = 1;
	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
		{
			if(perm[j] == i)
			{
				r_perm[i] = j; 
				break;
			}
		}
		size *= lda[i];
	}	
	for(i = 1; i < n; i++)
	{
		lda_s[i] = lda_s[i-1] * lda[i];
		ldb_s[i] = ldb_s[i-1] * lda[perm[i-1]];
	} 
       for(i = 0; i < n; i++)
        {
                temp[i] = ldb_s[r_perm[i]];
        }

        for(i=0; i < size; i++)
        {
                        B[getoff(i, lda, temp, n)] = A[i];
        }
}

int transpose_equal(const double *A, const double*B, int total_size){ 
	int error = 0; 
	const double *Atmp= A; 
	const double *Btmp= B; 
	for(int i=0;i < total_size ; ++i){ 
		double Aabs = (Atmp[i] < 0) ? -Atmp[i] : Atmp[i]; 
		double Babs = (Btmp[i] < 0) ? -Btmp[i] : Btmp[i]; 
		double max =  (Aabs < Babs) ? Babs : Aabs; 
		double diff = (Aabs - Babs); 
		diff = (diff < 0) ? -diff : diff; 
		if(diff > 0){ 
			double relError = (diff/max); 
			if(relError > 1.000000e-10){ 
				error += 1; 
			} 
		} 
	} 
	printf("Number of errors: %d\t",error); 
	return (error > 0) ? 0 : 1; 
}
int main(int argc, char* argv[])
{
	FILE * f;
	int dim[10];
	int permutation[10];
	int i, j;
	if(argc < 2)
	{
		fprintf(stderr, "Please input filename to read data\n");
		exit(0);
	}
	f = fopen(argv[1], "r");
	char line[255];
	char temp[50];
	int ndim;
	while(1)
	{
		size_t total = 1;
		fscanf(f, "%d", &ndim);
		if(feof(f)){
			break;
		}
		for(i = 0; i < ndim; i++)
		{
			fscanf(f, "%d", &dim[i]);
			if(feof(f)){
				break;
			}
			printf("%d ", dim[i]);
			total *= dim[i];
		}
		printf("\t");
		for(i = 0; i < ndim; i++)
		{
			fscanf(f, "%d", &permutation[i]);
			printf("%d ", permutation[i]);
		}

		size_t totalsize = total * sizeof(double);
		double* A = (double*) malloc(totalsize);
		double* B = (double*) malloc(totalsize);
		double*	B_trans = (double *) malloc(totalsize);
		for(i=0; i < total ; ++i){
			A[i] = (double)i;
		}
		int r_perm[10], r_dim[10];
		for(i = 0; i < ndim; i++)
		{
			for(j = 0; j < ndim; j++)
			{
				if(permutation[j] == i)
				{
					r_perm[i] = j; 
					//r_dim[i] = dim[j];
					break;
				}
			}
			r_dim[i] = dim[permutation[i]];
		}	
		transpose(ndim, A, B_trans, dim, permutation);
		transpose(ndim, B_trans, B, r_dim, r_perm);
		transpose_equal(A, B, total);
		free(A);
		free(B);
		free(B_trans);
		printf("\n");
	}
}