#include<stdio.h>
#include <math.h>

main()
{
double  I,b,w,k[4][4],D,H1,H2,H3;
int i,j,m,n,x,e,f,g,p,q,V,Z;
static double Q[4][4][4],c[4][4][4][4];
double d[4];

I=0;
double GaussPoints2[]={-0.57735026918963,
    0.57735026918963};
double GaussPtCeof2[]={1,
    1};
double GaussPoints3[]={-0.77459666924148,
    0,
    0.77459666924148};
double GaussPtCeof3[]={0.55555555555556,
    0.88888888888889,
    0.55555555555556};
double GaussPoints4[]={-0.86113631159405,
    -0.33998104358486,
    0.33998104358486,
    0.86113631159405};
double GaussPtCeof4[]={0.34785484513745,
    0.65214515486255,
    0.65214515486255,
    0.34785484513745};
double GaussPoints5[]={-0.90617984593866,
    -0.53846931010568,
    0,
    0.53846931010568,
    0.90617984593866};
double GaussPtCeof5[]={0.23692688505619,
    0.47862867049937,
    0.56888888888889,
    0.47862867049937,
    0.23692688505619};;
c[1][1][1][1]=39;
c[2][2][2][2]=39;
c[1][1][2][2]=30;
c[2][2][1][1]=30;
c[3][3][1][1]=27.6;
c[1][1][3][3]=27.6;
c[2][2][3][3]=27.6;
c[3][3][2][2]=27.6;
c[3][3][3][3]=28;
c[2][3][2][3]=51;
c[2][3][3][2]=51;
c[3][2][2][3]=51;
c[3][2][3][2]=51;
c[3][1][3][1]=51;
c[3][1][1][3]=51;
c[1][3][3][1]=51;
c[1][3][1][3]=51;
c[1][2][1][2]=49;
c[1][2][2][1]=49;
c[2][1][2][1]=49;
c[2][1][1][2]=49;
Q[1][2][3]=1;
Q[2][3][1]=1;
Q[3][1][2]=1;
Q[3][2][1]=-1;
Q[2][1][3]=-1;
Q[1][3][2]=-1;
I=0;
printf("input 4 subscript:\n");
scanf("%d%d%d%d",&i,&j,&m,&n);
for(Z=0;Z<1000;Z++)
   {for(V=0;V<1000;V++)
    {

b=GaussPoints4[V];
w=GaussPtCeof4[Z];
d[1]=((sqrt(1-(b*b)))*(cos(3.141593+3.141593*w)))/3;
d[2]=((sqrt(1-(b*b)))*(sin(3.141593+3.141593*w)))/3;
d[3]=b/11;


k[1][1]=((c[1][1][1][1])*d[1]*d[1]+(c[1][1][1][2])*d[1]*d[2]+(c[1][1][1][3])*d[1]*d[3]+(c[1][2][1][1])*d[2]*d[1]+(c[1][2][1][2])*d[2]*d[2]+(c[1][2][1][3])*d[2]*d[3]+(c[1][3][1][1])*d[3]*d[1]+(c[1][3][1][2])*d[3]*d[2]+(c[1][3][1][3])*d[3]*d[3]);

k[1][2]=((c[1][1][2][1])*d[1]*d[1]+(c[1][1][2][2])*d[1]*d[2]+(c[1][1][2][3])*d[1]*d[3]+(c[1][2][2][1])*d[2]*d[1]+(c[1][2][2][2])*d[2]*d[2]+(c[1][2][2][3])*d[2]*d[3]+(c[1][3][2][1])*d[3]*d[1]+(c[1][3][2][2])*d[3]*d[2]+(c[1][3][2][3])*d[3]*d[3]);

k[1][3]=((c[1][1][3][1])*d[1]*d[1]+(c[1][1][3][2])*d[1]*d[2]+(c[1][1][3][3])*d[1]*d[3]+(c[1][2][3][1])*d[2]*d[1]+(c[1][2][3][2])*d[2]*d[2]+(c[1][2][3][3])*d[2]*d[3]+(c[1][3][3][1])*d[3]*d[1]+(c[1][3][3][2])*d[3]*d[2]+(c[1][3][3][3])*d[3]*d[3]);

k[2][1]=((c[2][1][1][1])*d[1]*d[1]+(c[2][1][1][2])*d[1]*d[2]+(c[2][1][1][3])*d[1]*d[3]+(c[2][2][1][1])*d[2]*d[1]+(c[2][2][1][2])*d[2]*d[2]+(c[2][2][1][3])*d[2]*d[3]+(c[2][3][1][1])*d[3]*d[1]+(c[2][3][1][2])*d[3]*d[2]+(c[2][3][1][3])*d[3]*d[3]);

k[2][2]=((c[2][1][2][1])*d[1]*d[1]+(c[2][1][2][2])*d[1]*d[2]+(c[2][1][2][3])*d[1]*d[3]+(c[2][2][2][1])*d[2]*d[1]+(c[2][2][2][2])*d[2]*d[2]+(c[2][2][2][3])*d[2]*d[3]+(c[2][3][2][1])*d[3]*d[1]+(c[2][3][2][2])*d[3]*d[2]+(c[2][3][2][3])*d[3]*d[3]);

k[2][3]=((c[2][1][3][1])*d[1]*d[1]+(c[2][1][3][2])*d[1]*d[2]+(c[2][1][3][3])*d[1]*d[3]+(c[2][2][3][1])*d[2]*d[1]+(c[2][2][3][2])*d[2]*d[2]+(c[2][2][3][3])*d[2]*d[3]+(c[2][3][3][1])*d[3]*d[1]+(c[2][3][3][2])*d[3]*d[2]+(c[2][3][3][3])*d[3]*d[3]);

k[3][1]=((c[3][1][1][1])*d[1]*d[1]+(c[3][1][1][2])*d[1]*d[2]+(c[3][1][1][3])*d[1]*d[3]+(c[3][2][1][1])*d[2]*d[1]+(c[3][2][1][2])*d[2]*d[2]+(c[3][2][1][3])*d[2]*d[3]+(c[3][3][1][1])*d[3]*d[1]+(c[3][3][1][2])*d[3]*d[2]+(c[3][3][1][3])*d[3]*d[3]);

k[3][2]=((c[3][1][2][1])*d[1]*d[1]+(c[3][1][2][2])*d[1]*d[2]+(c[3][1][2][3])*d[1]*d[3]+(c[3][2][2][1])*d[2]*d[1]+(c[3][2][2][2])*d[2]*d[2]+(c[3][2][2][3])*d[2]*d[3]+(c[3][3][2][1])*d[3]*d[1]+(c[3][3][2][2])*d[3]*d[2]+(c[3][3][2][3])*d[3]*d[3]);

k[3][3]=((c[3][1][3][1])*d[1]*d[1]+(c[3][1][3][2])*d[1]*d[2]+(c[3][1][3][3])*d[1]*d[3]+(c[3][2][3][1])*d[2]*d[1]+(c[3][2][3][2])*d[2]*d[2]+(c[3][2][3][3])*d[2]*d[3]+(c[3][3][3][1])*d[3]*d[1]+(c[3][3][3][2])*d[3]*d[2]+(c[3][3][3][3])*d[3]*d[3]);


D=k[1][1]*k[2][2]*k[3][3]+k[1][2]*k[2][3]*k[3][1]+k[1][3]*k[2][1]*k[3][2]-k[1][1]*k[2][3]*k[3][2]-k[1][2]*k[2][1]*k[3][3]-k[1][3]*k[2][2]*k[3][1];
H1=0;
H2=0;
for(p=1;p<4;p++)
{   for(q=1;q<4;q++)
    {   for(x=1;x<4;x++)
          {   for(e=1;e<4;e++)
               {   for(f=1;f<4;f++)
                      { for(g=1;g<4;g++)
                           {H1+=((0.5*c[p][q][m][n]*d[j]*d[q]*Q[i][x][e]*Q[p][f][g]*k[x][f]*k[e][g])/(8*D));
                            H2+=((0.5*c[p][q][m][n]*d[i]*d[q]*Q[j][x][e]*Q[p][f][g]*k[x][f]*k[e][g])/(8*D));}

}}}}}

H3=(H1+H2);

I+=GaussPtCeof4[V]*GaussPtCeof4[Z]*(H3);}
}

printf("I=%f",I);
return 0;
}

搞这么复杂.

这是要干嘛啊！看着我都怕！

这个用matlab 会好点