e02race.c

/* nag_fit_pade_app (e02rac) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
Integerexit_status,i,k,m,ia,ib,ic,itmax,maxkm;
Nag_Root_Polishpolish;
NagErrorfail;
/* Arrays */
double*aa=0,*bb=0,*berr=0,*cc=0,*cond=0,*dd=0;
Complex*z=0;
Integer*conv=0;
INIT_FAIL(fail);
exit_status=0;
printf("nag_fit_pade_app (e02rac) Example Program Results\n");
k=4;
m=4;
ia=k+1;
ib=m+1;
ic=ia+ib-1;
maxkm=k>m?k:m;
/* Allocate memory */
if(!(aa=NAG_ALLOC(ia,double))||!(bb=NAG_ALLOC(ib,double))||
!(cc=NAG_ALLOC(ic,double))||!(dd=NAG_ALLOC(ia+ib,double))||
!(berr=NAG_ALLOC(maxkm,double))||
!(cond=NAG_ALLOC(maxkm,double))||
!(conv=NAG_ALLOC(maxkm,Integer))||!(z=NAG_ALLOC(maxkm,Complex))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* Power series coefficients in cc */
cc[0]=1.0;
for(i=1;i<=8;++i)
cc[i]=cc[i-1]/(double)i;
printf("\n");
printf("The given series coefficients are\n");
for(i=1;i<=ic;++i){
printf("%13.4e",cc[i-1]);
printf(i%5==0||i==ic?"\n":" ");
}
/* nag_fit_pade_app (e02rac).
 * Pade-approximants
 */
nag_fit_pade_app(ia,ib,cc,aa,bb,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_fit_pade_app (e02rac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
printf("\n");
printf("Numerator coefficients\n");
for(i=1;i<=ia;++i){
printf("%13.4e",aa[i-1]);
printf(i%5==0||i==ia?"\n":" ");
}
printf("\n");
printf("Denominator coefficients\n");
for(i=1;i<=ib;++i){
printf("%13.4e",bb[i-1]);
printf(i%5==0||i==ib?"\n":" ");
}
/* Calculate zeros of the approximant using
 nag_zeros_poly_real_fpml (c02abc) */
/* First need to reverse order of coefficients */
for(i=1;i<=ia;++i)
dd[ia-i]=aa[i-1];
/* nag_zeros_poly_real_fpml (c02abc).
 * Zeros of a polynomial with real coefficients
 */
itmax=30;
polish=Nag_Root_Polish_Simple;
nag_zeros_poly_real_fpml(dd,k,itmax,polish,z,berr,cond,conv,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_zeros_poly_real_fmpl (c02abc), 1st call.\n%s\n",
fail.message);
exit_status=1;
gotoEND;
}
printf("\n");
printf("Zeros of approximant are at\n");
printf(" Real part Imag part\n");
for(i=1;i<=k;++i)
printf("%13.4e%13.4e\n",z[i-1].re,z[i-1].im);
/* Calculate poles of the approximant using
 nag_zeros_poly_real_fpml (c02abc) */
/* Reverse order of coefficients */
for(i=1;i<=ib;++i)
dd[ib-i]=bb[i-1];
/* nag_zeros_poly_real_fpml (c02abc), see above. */
nag_zeros_poly_real_fpml(dd,m,itmax,polish,z,berr,cond,conv,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_zeros_poly_real_fpml (c02abc), 2nd call.\n%s\n",
fail.message);
exit_status=1;
gotoEND;
}
printf("\n");
printf("Poles of approximant are at\n");
printf(" Real part Imag part\n");
for(i=1;i<=m;++i)
printf("%13.4e%13.4e\n",z[i-1].re,z[i-1].im);
END:
NAG_FREE(aa);
NAG_FREE(bb);
NAG_FREE(cc);
NAG_FREE(dd);
NAG_FREE(berr);
NAG_FREE(cond);
NAG_FREE(conv);
NAG_FREE(z);
returnexit_status;
}