//                block lower Hessenberg sparse matrix


//      number of double multiplications and test operations
//                  compared with theoric bounds 							





# include <iostream.h>

# include <fstream.h>

# include <conio.h>

# include <stdlib.h>

# include <math.h>

# include <iomanip.h>



const int NNZTIL =10000;

const int N =1024;

const int D =4;

const int NN =4096;





typedef double mat[NN][NN];

typedef double dmat[D][D];

typedef double dvet[D];

typedef double xmat[NN][D];

typedef int init[NN];

typedef double reda[NN];





struct Sparsatil {

	double a[NNZTIL] ;

	int b[NNZTIL] ;

	int r[NN+1] ;

   };



Sparsatil s;

long seed;











void costra();

void crccs(int vet[],int vetc[],int vetr[], int p);

void crss(int vetc[],int vetr[],int test,int p);

void crcc(int vet[],int vetc[],int vetr[],int p);

void crs(int vetc[],int vetr[],int p);



void initrand(int v);

double rrand();

int irand(int n);

void scambia(int& a,int& b);

void quicksort(int vet[],int l, int r);

int contdop(int vet[],int p);



double infnorm(double vet[],int lun);

double norm2(double vet[],int lun);







int pow(int a,int b);

void base(int val,int l);

void recA(int val,int l,dmat xx);

void recC(int val,int l,dmat xx);



void inv(dmat xx);

int pos(int i);

void acc(int val,int l,xmat v);

void accT(int val,int l,xmat z,reda yyv);

void comX(int val,int l);



void spa_mv(reda q );

void stspa_pie(ofstream& out);















int p,n,d,nn;//n=2^p

int spn,bouv;

int m;

int tests,test;

xmat x,ww,yy,zz,au;

reda v,bb,b1,dy,yyv;



init ind;



int ncc=0,nca=0;

int asca=0,asca1=0;





int nmulacctv=0,nmulcomxv=0,nmulcorv=0,nmulm=0;

int nmulacct=0,nmulcomx=0,nmulinv=0;

int nmulcor=0,nmulxg=0;



int main() {

 ofstream out_stream;

 out_stream.open("spsp.txt");

 reda vet;

 cout.setf(ios::scientific);cout.precision(3);

 costra();

 stspa_pie(out_stream);





 for(int k=0;k<nn;k++){

 	bb[k]=k+1;b1[k]=k+1;

   dy[k]=0.;

 	for(int l=0;l<d;l++)x[k][l]=0.;

   }



for(int k=0;k<n;k++) for(int h=0;h<d;h++)ind[k*d+h]=(k+1)*d;



 base(0,0);



 double max;

 spa_mv(vet );

 out_stream.setf(ios::scientific);out_stream.precision(3);



 for(int i=0;i<nn;i++)dy[i]=vet[i]-b1[i];

 max=infnorm(dy,nn);

 cout<<" residual infnorm1= "<<max<<endl;;



 max=norm2(dy,nn);

 cout<<" residual 2norm= "<<max<<endl;

 cout<<endl;



out_stream<<" solution "<<endl;

for(int i=1;i<nn+1;i++){out_stream<<v[i-1]<<' ';if(!(i%6))out_stream<<endl;}

out_stream<<endl;



 cout<<" number of test operations in accT = "<<ncc<<endl;

 cout<<" number of test operations in acc  = "<<nca<<endl;

 cout<<" total number of test operations = "<<(ncc+nca)<<endl;

 cout<<" theoric bound  2 #SA= "<<2*spn<<endl<<endl;

 cout<<" number of double assignment operations in acc  =
"<<(asca+asca1)<<endl;

 cout<<" theoric bound d^2 (n-1) + #SA= "<<(spn+d*d*(n-1))<<endl<<endl;

 int nmulacctt=nmulxg+nmulacct+nmulacctv;

 int nmulcomxt=nmulcomx+nmulcomxv;

 int nmulinvmt=nmulinv+nmulm;

 int nmulcort=nmulcor+nmulcorv;

 cout<<" total number of double multiplications (factorization+solution)   = ";

 cout<<(nmulacctt+nmulcort+nmulcomxt+nmulinvmt)<<endl;

 cout<<" theoric bound : d^3 n/2 logn  + (d^2 n log n + #SA + d^2 n) (d+1) = ";

 cout<<(d*d*d*n/2*p+(d+1)*(spn + d*d*n*(1+p)))<<endl;





 out_stream.close();

 getch();

 return 0;

}









int pow(int a,int b){

 int i,c=1;for(i=1;i<=b;i++) c=c*a;

 return c;

}



int pos(int i){

 if(i<d)return 0; else return d;

}







void acc(int val,int l,xmat v){



   int nwb=n/pow(2,l);int nw=nwb*d;

	int ofm=val*nw;

   for(int k=0;k<nn;k++)for(int h=0;h<d;h++) v[k][h]=0.;

   for(int h=0;h<d;h++)for(int q=0;q<d;q++){asca++;

   	v[ofm+q][h]=s.a[s.r[ofm+h]+pos(ofm+h)+q];}

  	for(int h=0;h<d;h++) {

  	int im=ofm+h;int is=s.r[im]+pos(im)+d;

   int alf=s.b[is];

    while((alf<=(ofm+nw-1))&&(is<=(s.r[ofm+1+h]-1))){

    	v[alf][h]=s.a[is];

    	is++;asca1++;

      alf=s.b[is];

      nca++;

      }

    }

}



void accT(int val,int l,xmat z,reda yyv){

dmat as;

int ofm,ofs;

int nwb=n/pow(2,l);int nw0b=nwb/2;int nw=nwb*d;int nw0=nw0b*d;

ofm=val*nw;

int ofm1=(2*val+1)*nw0;

recC(val,l,as);

for(int i=0;i<nw0;i++)

	for(int j=0;j<d;j++){

   	au[i+ofm][j]=0.;

      for(int k=0;k<d;k++){

         nmulxg++;

         au[i+ofm][j]=au[i+ofm][j]+x[i+ofm][k]*as[k][j];}

      }

for(int k=0;k<nn;k++){

	yyv[k]=0.;

	for(int h=0;h<d;h++) z[k][h]=0.;}

int h;int alf;

for(int jj=0;jj<nw0b-1;jj++){

	for(int q=0;q<d;q++){

   	int j=ofm+jj*d+q;

      alf=ind[j];

      int alf1=s.r[j+1]; h=alf;int beta=s.b[h];

      while((h<alf1)&&(beta<(ofm1))){

      ncc++;h++;beta=s.b[h];}

      while((h<alf1)&&(beta<(ofm+nw))){

         for(int pp=0;pp<d;pp++){

            nmulacct++;

            z[beta][pp]=z[beta][pp]+s.a[h]*au[j][pp];}

         nmulacctv++;

         yyv[beta]=yyv[beta]+s.a[h]*v[j];

         h++;

         beta=s.b[h];

         ncc++;

         }

      ind[j]=h;

      }

   }

for(int q=0;q<d;q++){

	int j=ofm+(nw0b-1)*d+q;

   ofs=s.r[j]+pos(j)+d;

   alf=ofs;

   int alf1=s.r[j+1]; h=alf;int beta=s.b[h];

   while((h<alf1)&&(beta<(ofm1))){

   ncc++;h++;beta=s.b[h];}

  	while((h<alf1)&&(beta<(ofm+nw))){

  		for (int pp=0;pp<d;pp++){

         nmulacct++;

         z[beta][pp]=z[beta][pp]+s.a[h]*au[j][pp];}

      nmulacctv++;

      yyv[beta]=yyv[beta]+s.a[h]*v[j];

      h++;

      beta=s.b[h];

      ncc++;

      }

   ind[j]=h;

   }

}



void recC(int val,int l,dmat xx){

 int nwb=n/pow(2,l);int nw0b=nwb/2;int nw0=nw0b*d;

 int ofm1=(2*val+1)*nw0;

for(int j=0;j<d;j++){

	int rec=ofm1+j;

	for(int i=0;i<d;i++){

   	xx[i][j]=s.a[s.r[rec]+i];

   	}

	}

}



void comX(int val,int l){

	dmat cs,cs1;

   reda cuv;

   dvet cs1v;

   xmat cu;

 	int nwb=n/pow(2,l);int nw0b=nwb/2;int nw=nwb*d;int nw0=nw0b*d;

 	int ofm=val*nw;int ofm1=(2*val+1)*nw0;

	recC(val,l,cs);

	for(int i=0;i<d;i++){

		cs1v[i]=0.;

		for(int j=0;j<d;j++){

   		cs1[i][j]=0.;

         nmulcomxv++;

      	cs1v[i]=cs1v[i]+cs[i][j]*v[j+ofm1];

      	for(int k=0;k<d;k++){

            nmulcomx++;

         	cs1[i][j]=cs1[i][j]+cs[i][k]*x[k+ofm1][j];}

      	}

   	}

	for(int i=0;i<nw0;i++){

   	cuv[i]=0.;

		for(int j=0;j<d;j++){

   		cu[i][j]=0.;

         nmulcomxv++;

      	cuv[i]=cuv[i]+x[i+ofm][j]*cs1v[j];

			for(int k=0;k<d;k++){

            nmulcomx++;

         	cu[i][j]=cu[i][j]+x[i+ofm][k]*cs1[k][j];}

		}

}

for(int i=0;i<nw0;i++){

	v[ofm+i]=v[ofm+i] -cuv[i];

	for(int j=0;j<d;j++)x[ofm+i][j]= -cu[i][j]; }

}





void recA(int val,int l,dmat xx){

	int nw=n*d/pow(2,l);int ofm=val*nw;



	int poso=pos(ofm);

	for(int j=0;j<d;j++){

		int recc=ofm+j;

		for(int i=0;i<d;i++){

   		xx[i][j]=s.a[s.r[recc]+i+poso];

         }

      }

}



void inv(dmat xx){

	dmat b;double t;

	for(int i=0;i<d;i++)for(int k=0;k<d;k++)b[i][k]=0.;

	for(int i=0;i<d;i++)b[i][i]=1.;

		for(int k=0;k<d;k++){

			for( int i=(k+1);i<d;i++) {

 nmulinv++;

 t=xx[i][k]/xx[k][k];

for( int j=k+1;j<d;j++){

	nmulinv++;

	xx[i][j]=xx[i][j]-t*xx[k][j];}



for( int h=0;h<k;h++){

   nmulinv++;

   b[i][h]=b[i][h]-t*b[k][h];}

b[i][k]=b[i][k]-t;

            }

         }

	for(int  k=0;k<d;k++){

      nmulinv++;

   	b[d-1][k]=b[d-1][k]/xx[d-1][d-1];

   	for(int i=1;i<d;i++){

			int di=d-1-i;

      	for(int j=0;j<i;j++){

               nmulinv++;

               b[di][k]=b[di][k]-xx[di][di+j+1]*b[di+j+1][k];}

        nmulinv++;

        b[di][k]=b[di][k]/xx[di][di];

         	}

   	}

  	for(int i=0;i<d;i++)for(int k=0;k<d;k++)xx[i][k]=b[i][k];

}









void base(int val,int l){

 dmat xx;

 dmat ras1;dvet vv;

 int nwb=n/pow(2,l);int nw0b=nwb/2;int nw=nwb*d;int nw0=nw0b*d;

 int ofm=val*nw;

 int ofm1=(2*val+1)*nw0;

  if (l==p){

 	if (val==0){

   	recA(val,l,xx);

      inv(xx);

      for(int k=0;k<d;k++){vv[k]=0.;for(int j=0;j<d;j++){

         nmulm++;

         vv[k]=xx[k][j]*bb[j]+vv[k];}}

      for(int k=0;k<d;k++){

      	v[ofm+k]=vv[k];

      	for(int j=0;j<d;j++)x[ofm+k][j]=xx[k][j];

         }

            }

	else { //val!=0

      for(int k=0;k<d;k++){

      	for(int j=0;j<d;j++)xx[k][j]=ww[ofm+k][j];}

      inv(xx);

      for(int k=0;k<d;k++){vv[k]=0.;for(int j=0;j<d;j++){

         nmulm++;

         vv[k]=xx[k][j]*bb[j+ofm]+vv[k];}}

      for(int k=0;k<d;k++){

      	v[ofm+k]=vv[k];

      	for(int j=0;j<d;j++)x[ofm+k][j]=xx[k][j];

         }

      }

	}

 else{

	base(2*val,l+1);

   accT(val,l,yy,yyv);

	if (val){

    for(int h=0;h<d;h++){for(int k=0;k<d;k++)

         ras1[k][h]=au[k+ofm][h];

         }

    for(int k=0;k<nw0;k++){

         for(int h=0;h<d;h++){

            nmulcorv++;

				yyv[k+nw0+ofm]=yyv[k+nw0+ofm]-zz[k+nw0+ofm][h]*v[h+ofm];

         		for(int q=0;q<d;q++){

                  nmulcor++;

               	yy[k+nw0+ofm][h]=yy[k+nw0+ofm][h]-zz[k+nw0+ofm][q]*ras1[q][h];

            	}

         }

      }

		for(int k=0;k<nw0;k++){

      	for(int h=0;h<d;h++)

         zz[k+ofm1][h]=yy[k+ofm1][h];

         }

      }

   else for(int k=0;k<nw0;k++){

   	for(int h=0;h<d;h++)

   		zz[k+ofm1][h]=yy[k+ofm1][h];

      }

   acc(2*val+1,l+1,ww);

   for(int k=0;k<nw0;k++){

   	bb[k+ofm1]=bb[k+ofm1]-yyv[ofm1+k];

      for(int h=0;h<d;h++)

   	ww[k+ofm1][h]=ww[k+ofm1][h]-zz[ofm1+k][h];

      }

	base(2*val+1,l+1);

   comX(val,l);

	}

}





void spa_mv(reda q ){

	int i,j,count,lc;

	for( i=0;i<nn;i++){

		q[i]=0.;

		for (j=0;j<(nn-1);j++){

      	count=s.r[j]; lc=s.r[j+1]-1;

      	while ((s.b[count]<i)&&(count<lc))count++;

      	if (s.b[count]==i)q[i]=q[i]+s.a[count]*v[j];

         }

     j=nn-1;

     count=s.r[j]; lc=m-1;

     while ((s.b[count]<i)&&(count<lc))count++;

     if (s.b[count]==i){

     	q[i]=q[i]+s.a[count]*v[j];

      }

     }

   }



void crccs(int vet[],int vetc[],int vetr[],int p){

for( int i=0;i<p;i++) vetc[i]=0;

for( int i=0;i<p;i++) {

   int lb=0, u=n-1,j=0,ofm=d,l=0,uv=0;

   while(vet[i]>=u*d*d){uv=u;lb=u;u=u+n-2-l;j=j+d;ofm=ofm+d;l++;}

   int le=lb*d*d, ue=uv*d*d+(n-l-1)*d,k=n-l-1;

   while(vet[i]>=ue){le=ue;ue=ue+k*d;j=j++;

   	}

	vetc[i]=j;vetr[i]=vet[i]-le+ofm;

	}

}







void crss(int vetc[],int vetr[],int test,int p){

int k;int j;

s.r[0]=0;int jk=0;k=0;

if(test==1){

	for(int j=0;j<d;j++){

		s.r[j]=jk;

   	for(int h=0;h<d;h++){

   		s.a[jk]=0.;if((j==0)&&(h==0))s.a[jk]=0.01;if((j==1)&&(h==1))s.a[jk]=1.;

   		s.b[jk]=h;jk++;}

		int alf=vetc[k];

		while((alf==j)&&(k<p)){s.a[jk]=1.;s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

   	}

 		int pp;

 		for(int jj=1;jj<n;jj++)for(int q=0;q<d;q++){

   		pp=(jj-1)*d;j=d*jj+q;

			s.r[j]=jk;

			for(int h=0;h<d;h++){s.a[jk]=0.;if(q==0)s.a[jk]=1.;s.b[jk]=h+pp;jk++;}

   		for(int h=0;h<d;h++){s.a[jk]=0.;

            if((q==0)&&(h==0))s.a[jk]=0.01;if((q==1)&&(h==1))s.a[jk]=1.;

         	s.b[jk]=h+pp+d;jk++;}

   		int alf=vetc[k];

			while((alf==j)&&(k<p)){s.a[jk]=1.;s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

		}

 	}



if(test==2){

	for(int j=0;j<d;j++){

		s.r[j]=jk;

   	for(int h=0;h<d;h++){s.a[jk]=0.;

   		if(j==h)s.a[jk]=1.;else s.a[jk]=0.01;

   		s.b[jk]=h;jk++;}

		int alf=vetc[k];

		while((alf==j)&&(k<p)){s.a[jk]=0.99;s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

   	}

 		int pp;

 		for(int jj=1;jj<n;jj++)for(int q=0;q<d;q++){

   		pp=(jj-1)*d;j=d*jj+q;

			s.r[j]=jk;

			for(int h=0;h<d;h++){s.a[jk]=0.;

				if((q==0)&&(h==1))s.a[jk]=0.99;

         	s.b[jk]=h+pp;jk++;}

   		for(int h=0;h<d;h++){

   			s.a[jk]=0.;

   			if(q==h)s.a[jk]=1.;else s.a[jk]=0.01;

         	s.b[jk]=h+pp+d;jk++;}

   		int alf=vetc[k];

			while((alf==j)&&(k<p)){s.a[jk]=0.99;s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

		}

 	}

if(test==3){

	for(int j=0;j<d;j++){

		s.r[j]=jk;

   	for(int h=0;h<d;h++){

   		s.a[jk]=0.;

   		if(j==h)s.a[jk]=1.;else {if(h==1)s.a[jk]=-0.00049;else s.a[jk]=-490.;}

   		s.b[jk]=h;jk++;}

		int alf=vetc[k];

		while((alf==j)&&(k<p)){s.a[jk]=-490.;

   		s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

   	}

 		int pp;

 		for(int jj=1;jj<n;jj++)for(int q=0;q<d;q++){

   		pp=(jj-1)*d;j=d*jj+q;

			s.r[j]=jk;

			for(int h=0;h<d;h++){

				s.a[jk]=0.;

         	if((q==0)&&(h==1))s.a[jk]=-0.00049;

         	s.b[jk]=h+pp;jk++;}

   		for(int h=0;h<d;h++){

   			s.a[jk]=0.;

            if(q==h)s.a[jk]=1.;else {if(h==1)s.a[jk]=-0.00049;else
s.a[jk]=-490.;}

         	s.b[jk]=h+pp+d;jk++;}

   		int alf=vetc[k];

			while((alf==j)&&(k<p)){s.a[jk]=-490.;s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

		}

 	}

}



double infnorm(double vet[],int lun){

double max, au;

 max=fabs(vet[0]);

 for(int i=1;i<lun;i++){au=fabs(vet[i]);if (max<au)max=au;}

 return max;}



double norm2(double vet[],int lun){

double n2=0.;

 for(int i=0;i<lun;i++){double t=vet[i];n2=n2+t*t;}

 n2=sqrt(n2);

 return n2;

}

void costra(){

int c[NNZTIL];int z[NNZTIL];int t[NNZTIL];



cout<<"  p=? "<<endl;

cin>>p;

n=pow(2,p);

cout<<"  n= "<<n<<' ';

cout<<"  tests(=0 ,=1(examples in Von Matt-Stewart) )=?  ";

cin>>tests;

if(!tests){

cout<<" d=? "<<endl;

cin>>d;

cout<<"  d= "<<d<<' ';

nn=n*d; bouv=n*(n-1)*d*d/2;

cout<<"  nn= "<<nn<<endl<<endl;

cout<<"  test(=0  random,=1 special)=?  ";

cin>>test;



if(!test){  // random

	cout<<"  spn < "<<bouv<<", "<<(NNZTIL-n*d*d-(n-1)*d*d)<<"   ";

   cout<<" if spn= "<<bouv<<" SA full"<<endl;

	cin>>spn;

   if(spn!=bouv){

	initrand(0);

	for( int i=0;i<spn;i++){

   	c[i]=irand(bouv);

   	}

      quicksort(c,0,spn-1);

	int spnn=spn-contdop(c,spn);

   quicksort(c,0,spnn-1);

	cout<<" # SA= "<<spnn<<endl;

	spn=spnn;

	m=spn+n*d*d+(n-1)*d*d;

	quicksort(c,0,spn-1);

   }else{for( int i=0;i<spn;i++)c[i]=i;m=spn+n*d*d+(n-1)*d*d;}

}

else{ // special

   for(int k=0;k<spn;k++)c[k]=0;

   int q=0; spn=(n-1)*d*d;

   for(int pp=0;pp<d;pp++)for(int k=0;k<nn-d;k++){

   	c[q]=q;

      q++;}

   cout<<" # SA= "<<spn<<endl;

	m=spn+n*d*d+(n-1)*d*d;

	}

crcc(c,z,t,spn);

crs(z,t,spn);

}

else {

	d=2;

	cout<<"  d= "<<d<<' ';

	nn=n*d;

	cout<<"  nn= "<<nn<<endl<<endl;

	spn=n-1;

	for(int k=0;k<spn;k++)c[k]=0;

   	int q=nn-2;

   	for(int k=0;k<spn;k++){

   		c[k]=q;

      	q=q+2*(nn-(2*(k+1)+1));

      	}

	cout<<"  test(=1 ,=2 ,=3  )=?  ";

	cin>>test;

	m=spn+n*d*d+(n-1)*d*d;

	cout<<" #SA= "<<spn<<endl;

	crccs(c,z,t,spn);

	crss(z,t,test,spn);

	}

}



void crcc(int vet[],int vetc[],int vetr[],int p){

	for( int i=0;i<p;i++) vetc[i]=0;

	for( int i=0;i<p;i++) {

   	int lb=0, u=n-1,j=0,ofm=d,l=0,uv=0;

   	while(vet[i]>=u*d*d){uv=u;lb=u;u=u+n-2-l;j=j+d;ofm=ofm+d;l++;}

   	int le=lb*d*d, ue=uv*d*d+(n-l-1)*d,k=n-l-1;

   	while(vet[i]>=ue){le=ue;ue=ue+k*d;j=j++;

   	}

	vetc[i]=j;vetr[i]=vet[i]-le+ofm;

	}

}





void crs(int vetc[],int vetr[],int p){

	int k;int j;

	s.r[0]=0;int jk=0;k=0;

	for(int j=0;j<d;j++){

		s.r[j]=jk;

   	for(int h=0;h<d;h++){s.a[jk]=0.;//-rrand();

   		s.b[jk]=h;jk++;}

		int alf=vetc[k];

		while((alf==j)&&(k<p)){s.a[jk]=-1.;//-rrand();

   	s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

   	}

 	int pp;

 	for(int jj=1;jj<n;jj++)for(int q=0;q<d;q++){pp=(jj-1)*d;j=d*jj+q;

		s.r[j]=jk;



			for(int h=0;h<d;h++){

				s.a[jk]=-1.;//-rrand();

         	s.b[jk]=h+pp;jk++;}



   		for(int h=0;h<d;h++){

   			s.a[jk]=0.;//-rrand();

         	s.b[jk]=h+pp+d;jk++;}

   	int alf=vetc[k];

		while((alf==j)&&(k<p)){s.a[jk]=-1;//-rrand();

   	s.b[jk]=vetr[k];jk++;k++;alf=vetc[k];}

		}

	for(int k=0;k<d;k++){

		double sumn=0.;

		for(int i=0;i<k;i++)sumn=sumn+ s.a[s.r[k]+i];

		for(int i=k+1;i<d;i++)sumn=sumn+ s.a[s.r[k]+i];

		for(int h=s.r[k]+d;h<s.r[k+1];h++)

   		sumn=sumn+ s.a[h];

   	s.a[s.r[k]+k]= -sumn+1.;

		}

	for(int j=1;j<n;j++){for(int l=0;l<d;l++){

		int k=j*d+l;double sumn=0.;

		for(int i=0;i<d;i++)sumn=sumn+ s.a[s.r[k]+i];

		for(int i=d;i< d+l  ;i++)sumn=sumn+ s.a[s.r[k]+i];

		for(int i=  d+l +1;i<2*d;i++)sumn=sumn+ s.a[s.r[k]+i];

   	for(int h=s.r[k]+2*d;h<s.r[k+1];h++)sumn=sumn+ s.a[h];

		s.a[s.r[k]+d+l]= -sumn+1.;

		}

	}

}





void scambia(int& a,int& b){

	int temp=a;

	a=b;

	b=temp;

}



void quicksort(int vet[],int l, int r){

int pp,li,ri;

int i;

li=l; ri=r; pp=vet[(l+r)/2];



for(i=l; i<=r; i++)



do{

	while(vet[ri]>pp) ri--;

	while(vet[li]<pp) li++;

	if(li<=ri){

		scambia(vet[li],vet[ri]);



      li++;

		ri--;



		}

	}

while(ri>=li);

if(l<ri) quicksort(vet,l,ri);

if(li<r) quicksort(vet,li,r);

}



int contdop(int vet[],int p){

int cont=0;int k,q,r;int sd=p;

	k=0;

   while(k<sd){

	q=k+1;

   	while(q<sd){

 		r=q+1;

         if(vet[k]==vet[q]){

         	cont++;

            for(int h=q;h<sd-1;h++)vet[h]=vet[h+1];

            sd--;

            r--;

            }

         q=r;

         }

   k++;}

return cont;

 }





void initrand(int v){

seed=v;}



double rrand(){

seed=(seed*15937+33503)%65521;

return double(seed)/65521;

}



int irand(int n){

return int(rrand()*n);

}





void stspa_pie(ofstream& out){

	int i,j,count,lc;

   out.setf(ios::scientific);out.precision(1);

	for( i=0;i<nn;i++){

		out<<endl;

		for (j=0;j<(nn-1);j++){

      	count=s.r[j]; lc=s.r[j+1]-1;

      	while ((s.b[count]<i)&&(count<lc))count++;



         if (s.b[count]==i)out<<s.a[count]<<' ';

         else out<<"       "<<' ';

         }

     j=nn-1;

     count=s.r[j]; lc=m-1;

     while ((s.b[count]<i)&&(count<lc))count++;

     if (s.b[count]==i){

      out<<s.a[count]<<' ';

      }

     else out<<"    "<<' ';

     }

   out<<endl;

   }