C言語 プログラミング DK法

質問

至急回答お願いしたいです！
プログラミング初心者です．
C言語のプログラミングについての質問です．
代数方程式
P_８ (z)=z^8-36z^7+546z^6-4536z^5+22449z^4-67284z^3+118124z^2-109584z+40320
の解をDK法（Durand-Kerner法）で求めるプログラムをC言語で書いていただきたいです．
できれば初期値の与え方は，すべての解を内に含むような半径Rの円周上に配置する方法でお願いしたいです．
よろしくお願いします．

回答No.2

f272

ちょっと検索すると例えば
http://ysrken.blog.fc2.com/blog-entry-21.html
が見つかって，ここにあったDKA法.cppをC言語で書き直せば，こんな感じ。

#include <complex.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
typedef double _Complex dcomplex;
const dcomplex im=0.0+1.0*_Complex_I;

int main() {
FILE *fpi=fopen("input.txt","r");
FILE *fpo=fopen("output.txt","w");
int n,i, j, k, L;
double maxerror,err=DBL_MAX,err_,real,imag,R,R_;
dcomplex *a,*z,*z0,f,fp,z0i,ci,zc,*b;
fscanf(fpi,"%d %lf",&n,&maxerror);
a=(dcomplex*)malloc((n+1)*sizeof(dcomplex));
for(i=0;i<=n;i++) {
fscanf(fpi,"%lf %lf",&real,&imag);
a[i]=real+imag*_Complex_I;
}
z=(dcomplex*)malloc((n)*sizeof(dcomplex));
if(n==1) {
z[0]=-a[1]/a[0];
}
else {
R=0.0;
ci=2.0*M_PI/n*im;
zc=-a[1]/((double)n*a[0]);
b=(dcomplex*)malloc((n+1)*sizeof(dcomplex));
for(k=0;k<=n;k++)
b[k]=a[k];
for(L=0;L<n;L++) {
for(k=1;k<=n-L;k++)
b[k]=b[k-1]*zc+b[k];
R_=n*cpow(cabs(b[n-L]/a[0]),(1/(n-L)));
if(R_>R)
R=R_;
}
for(i=0;i<n;i++) {
z[i]=zc+R*cexp(ci*(i+1-0.75));
}
z0=(dcomplex*)malloc((n)*sizeof(dcomplex));
while(err>maxerror) {
for(i=0;i<n;i++)
z0[i]=z[i];
err=0.0;
for(i=0;i<n;i++) {
z0i=z0[i];
f=a[0];
for(j=1;j<=n;j++)
f=f*z0i+a[j];
fp=a[0];
for(j=0;j<i;j++)
fp*=z0i-z0[j];
for(j=i+1;j<n;j++)
fp*=z0i-z0[j];
z[i]=z0i-f/fp;
err_=cabs(z[i]-z0i);
if(err<err_)
err=err_;
}
}
}
for(i=0;i<n;i++) {
fprintf(fpo,"x[%d]=",real);
real=floor(creal(z[i])/maxerror+0.5)*maxerror;
imag=floor(cimag(z[i])/maxerror+0.5)*maxerror;
if(real!=0.0)
fprintf(fpo,"%f",real);
if((real!=0.0)&&(imag>0.0))
fprintf(fpo,"+");
if((imag!=0.0)&&(imag!=1.0))
fprintf(fpo,"%f",imag);
if(imag!=0.0)
fprintf(fpo,"i");
if((real==0.0)&&(imag==0.0))
fprintf(fpo,"%f",0);
fprintf(fpo,"\n");
}
free(a);
free(z);
return 0;
}

input.txtはこれ。
8 0.00001
1 0
-36 0
546 0
-4536 0
22449 0
-67284 0
118124 0
-109584 0
40320 0

output.txtはこうなる。
x[0]=8.000000
x[0]=6.000000
x[1]=4.000000
x[0]=2.000000
x[0]=1.000000
x[0]=3.000000
x[1]=5.000000
x[0]=7.000000

gcc 4.8.2を使って計算してみた。

回答No.1

f272

ネットで検索すれば，色々と見つかると思うよ。

補足

回答ありがとうございます．
私もネットで色々と検索してみたのですが，なかなかコンパイルできずに困っています．
説明不足で申し訳ありません．