求助:关与C语言的计算程序,解卡尔丹(x**3+p*X+q=0)一元三次的实根,结果存入xa[],出不来正确结果。
/**************晚上**************/
#include<stdio.h>
#include<math.h>
#include<float.h>
#define Real double
#define PI 3.14159265354
#define Real double
#define Real_MIN DBL_MIN
/*double PI2(double xx){return 2*PI;}*/
/*解卡尔丹(x**3+p*X+q=0)一元三次的实根,结果存入xa[],实
根个数为*n*/
void getslovp3l(real p,real q,real xa,int n)
{real d,r,s;
p*=(1./.3);
q*=0.5;
d=q*q+p*p*p;
if(d>0.)
{*n=1;d=sqrt(d);
xa[0]=(((d-q)>=0.?pow(d-q,1./3.):-pow(q-d,1./3.))+((-q-d)>=0.?pow(-q-d,1./3.):-pow(q+d,1./3./)));
}
else if(d==0)
{*n=3;
if(p==q==0.)xa[0]=xa[1]=xa[2]=0.;
else {r=(q.0?-pow(q,1./3.):pow(-q,1./3.));
xa[0]=r+r;
xa[1]=xa[2]=-r; }
}
else
{*n=3;
r=pow(q*q-d,1./6.);
s=(1./3.)*atan2(sqrt(-d),-q);
xa[0]=r*cos(s);xa[0]+=xa[0];
xa[1]=r*cos(s+PI2/3.);xa[1]+=xa[1];
xa[2]=r*cos(s-PI2/3.);xa[2]+=xa[2];
}
}
/*解一元二次方程的实根,结果存入x[],实根个数为*n*/
void getslov2l(Real a,Real b,Real c,Real *x,int *n)
{Real d;
d=b*b-4.*a*c;
if(d>0.)
{*n=2;
if(b>=0)x[0]=(-b-sqrt(d))/(2.*a);
else x[0]=(-b+sqrt(d))/(2.*a);
if(fabs(x[0])>Real_MINx[1]=c/(a*x[0]);
else x[1]=-b/a;
}
else if(d==0){*n=2;x[0]=x[1]=-b/(2.*a);}
else {*n=0;return;}
}
/*解任意一元三次方程ax**3+b*X**2+c*x+d=0(a!=0)方程的实
根,其中*n取值为0-3,为实根个数*/
void getslovp3l(Real a,Real b,Real c,Real d,real *x,int *n)
{Real h;
b/=a;c/=a;d/=a;
h=b*(1./3.);
getslovp3l(c-b*b*(1./3.),2.*b*b*b*(1./27)-b*c*(1./3)+d,x,n);
x[0]-=h;x[1]-=h;x[2]-=h;
}
*解ax**3+b*X**2+c*x+d=0一元三次方程的实根,(其中包括a=0退化为低次方程的情况)
其中*n取值为0-3时为实根个数,取值为4时表示所有实数都是该方程的解*/
void getslovp3lall(Real a,Real b,Real c,Real d,real *x,int *n)
{Real h;
if(fabs(a)>Real_MIN)
{
b/=a;c/=a;d/=a;
h=b*(1./3.);
getslop3l(c-b*b*(1./3.),b*b*b*(2./27.)-b*c*(1./3.)+d,x,n);
x[0]-=h;x[1]-=h;x[2]-=h;
}
else if(fabs(b)>Real_MIN)
{
getslov2l(b,c,d,x,n);
}
else if(fabs(c)>Real_MIN)
{
*n=1;x[0]=d/c;
}
else if(d==0.)*n=4;
else *n=0;
}
mian[]
{double a,b,c,d,x[3];
char str1[10],str2[10];
float aa;
int i,n,tn,tmin,lmin.lmax,tmax;
FILE *fp;
/*scanf("%f%f%f%f,&a,&b,&c,&d);*/
/*scanf("%f,&a);*/
/*scanf("%f,&b);*/
float m0,s0,e0,a0,fpp,fmax,detp,d0,l,lm,af;
float c0,f1,f6,f,g,t,v,det;
float t1,t4,t6,v1,v4,v6;
Real g1,g4,g6;
if((fp=fopen("jcsj.txt","r"))==NULL{puts("cannot openfile jg.txt\n");}
fscanf(fp,"%f%f%f%f%f%f%f%f%f%d%d',&m0,&s0,&d0,&fmax,&e0,&a0,&af,&k,&lm,&tim,&tmax);
fclose(fp);
/*mo 导线计算单位质量(kg/km)*/
/*s0 导线计算截面积(mm~2)*/
/*d0 导线计算直径(mm)*/
/*fmax 导线拉断力(N)*/
/*e0 导线弹性系数(kg/mm~2)*/
/*a0 导线线膨胀系数(1/℃)*/
/*af 风速不均匀系数,绝缘导线af=1*/
/*lmin 导线需计算的最小档距(m)*/
/*lmax 导线需计算的最大档距(m)*/
/*tmin* 导线需计算的最低温度(℃)*/
/*tmax 导线需计算的最高温度(℃)*/
detp=0.;
t1=-5.;t4=20.;t6=20.;v1=0.;v4=25.;v6=0;/*气象条件*/
tn=0;
detp=(detp==0)?fmax/s0:detp;
fmax=(fmax==0)?detp*s0:fmax;
g1=9.80665*m0/(1000.*s0);
c0=(d0>=17)?1.1:1.2;
/*不考虑风压高度变化系数情况下*/
g4=0.6128*v4*v4*c0*d0*af/(s0*1000.);
g6=sqrt(g1*g1+g4*g4);
a=1.;
c=0.;
if((fp=fopen("jp.txt","w"))==NULL)(puts("cannot openfile jg.txt\n");}
fprintf(fp,"\n导线单位质量(kg/km)m=%g",m0);fprintf(fp,"\n 导线截面(mm~2) s=%g",s0);
fprintf(fp,"\n导线直径(mm)d=%g",d0);fprintf(fp,"\n 导线弹性系数(N/mm~2) E=%g",e0);
fprintf(fp,"\n导线线膨胀系数(1/℃)α=%g",a0);fpritf(fp,"\n导线破坏应力&p=%g",dept);
fprintf(fp,"\n最大使用应力detp=%g",detp/k);fprint(fp,"\n安全系数k=%g",k);
fprintf(fp,"\n自重比载个g1=%g",g1);fprint(fp,"\n风比载g4=%g,g4);
fprint(fp,"\n无冰综合比载g6=%g,g6);
for(tn=tmin;tn<=tman;tn=tn+5)
{printf("\n tn=%d",tn);
fprintf(fp,"\n温度tn=%d",tn);
{
f1=e0*g1*g1*lm*lm*(24*(detp/k)*(detp/k))-detp/k+0*e0*t1);
f4=e0*g6*g6*lm*lm*(24*(detp/k)*(detp/k))-detp/k+0*e0*t4);
f6=e0*g1*g1*lm*lm*(24*(detp/k)*(detp/k))-detp/4+0*e0*t6);
f=f1;t=t1;v=v1;g=g1;det=detp/k;
if(f1<f4){t=t4;v=v4;f=f4;g=g6;}
if(f<f6){t=t6;v=v6;f=f6;g=g1;det=detp/4;}
b=e0*g*g*lm*lm/(24.*det*det)-det+a0*e0)*(tn-t);
d=-e0*g1*lm*lm/24.;
/* forintf(fp,"\n a= %d",b);
fprintf(fp," b= %d",d)*/
getslov3lall(a,b,c,d,x,&n);
if(n!=4)for(i=0;i<n;++i)
{fpp=g1*lm*lm/(8*x[i]);
if(x[i]>0){
printf("\n1=%5.1f",lm);
printf(" det=%7.3f",x[i]);
printf(" f=%6.3f",fpp);
fprintf(fp,"\n=档距l=%5.1f",lm);
fprintf(fp," 应力&=%7.3f",lm);
fprintf(fp," 弧垂f=%6.3f",fpp);
}
else if(n==4)print("\n all Real numbers 解为所有实数");
}
}
fclose(fp);
getch0;
}