/* RESOLUCIO D'UN SISTEMA D'EQUACIONS LINEALS AMB PIVOTATGE PARCIAL*/
#define n 31
#include "math.h"
char m,p[n]; double a=-1,bb=1,h;
main()
{double U[n][n], b[n], x[n],f(double); char i,j;
 clrscr();
 printf(
 "PROGRAMA PER CALCULAR ELS PESOS DE NEWTON-COTES AMB a=-1, b=1, m<31\n");
 printf(
 "(c) Rafael Pla López, Dpt.Matemàtica Aplicada, Universitat de València");
 printf("\nSoftware de Domini Públic: ");
 printf("de lliure distribució sempre que es faça gratuitament\n");
 do {printf("m="); scanf("%d",&m);} while(m>=31); h=(bb-a)/m;
 for(i=0;i<=m;i++) p[i]=i;
 triangul(U,b);
 for(i=m;i>=0;i--)
	{x[i]=b[p[i]]; for(j=i+1;j<=m;j++) x[i]-=U[p[i]][j]*x[j];
	 printf("W%2d=%9.6lg  ",i,x[i]/=U[p[i]][i]);
	 if(i%5==0) printf("\n");
	}
 getch();
}

triangul(A,b)
double A[n][n],b[n];
{char k;
 introduir(A,b);
 for(k=0;k<m;k++) pivot(A,b,k);
}

pivot(A,b,k)
double A[n][n],b[n]; char k;
{char i,j; double l[n];
 for(i=m-1;i>=k;i--) if(fabs(A[p[i]][k])<fabs(A[p[i+1]][k]))
			{j=p[i+1]; p[i+1]=p[i]; p[i]=j;}
 if(A[p[k]][k]==0) {printf("\nKO:A%d,%d=%.3e\n",p[k],k,A[p[k]][k]); exit(1);}
 for(i=k+1;i<=m;i++) 
	{l[p[i]]=-A[p[i]][k]/A[p[k]][k];
	 for(j=0;j<=m;j++) A[p[i]][j]+=A[p[k]][j]*l[p[i]];
	 b[p[i]]+=b[p[k]]*l[p[i]];
	}
}

introduir(P,q)
double P[n][n],q[n];
{char i,j,k; double p1,p2;
 for(i=0;i<=m;i++) 
	{for(j=0;j<=m;j++) {P[i][j]=1; for(k=0;k<i;k++) P[i][j]*=a+j*h;}
	 p1=p2=1; for(k=0;k<=i;k++) {p1*=bb; p2*=a;}
	 q[i]=(p1-p2)/(i+1);}
}