Category : Files from Magazines
Archive   : CUJ9306.ZIP
Filename : 1106032A

 
Output of file : 1106032A contained in archive : CUJ9306.ZIP
/******************************************************
* RPFT.C -- Curve fitting with optional extrapolation.
* Fits either polynomial or rational function. Based
* on algoritms defined in Press, et al., "Numerical
* Recipes", 2nd ed., Cambridge University Press, 1992
* Developed using the Borland C compiler.
*
* Lowell Smith
* 3368 Crestwood Drive
* Salt Lake City, UT 84109-3202
*****************************************************/
#include
#include
#include
#include
#include
#include
#include

#define NPFAC 8
#define PI 3.141592653589793
#define PIO2 (PI/2.0)
#define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
void ratlsq(double xs[],double fs[],int npt,int mm,
int kk, double cof[], double *dev);
double *dvector(long nl, long nh, const char *ner);
double **dmatrix(long nrl,long nrh,long ncl,long nch,
const char *ner);
void free_dvector(double *v, long ncl);
void free_dmatrix(double **m, long nrl, long ncl);

/*****************************************************/
void main(int argc, char **argv)
{ double a,b,*xs,*fs,cof[42],dev=0.,sumd,x,yy;
int i,j,nn,nd,ncof,dig,npt=0,xln,yln;
char nar[90];
char filenam[80];
void commd(int argc, char**argv, double *a,double *b,
int *nn, int *nd, int *dig, int *xln,
int *yln, char *filenam);
void inpt(int nn,int npt,double a, double b, int xln,
int yln, double *xs,double *fs,char filenam[]);
void pcshft(double a,double b,double d[],int n);
void chebpc(double c[],double d[],int n);

commd(argc,argv,&a,&b,&nn,&nd,&dig,&xln,&yln,
filenam);
if (xln) { a=log(a); b=log(b); }
if (nd) /* Fit rational function */
{ ncof=nn+nd+1; npt=NPFAC*ncof;
xs=dvector(1L,(long)npt,"xs in main");
fs=dvector(1L,(long)npt,"fs in main");
inpt(0, npt, a, b, xln, yln, xs, fs, filenam);
ratlsq(xs,fs,npt,nn,nd,&cof[0],&dev);
fprintf(stderr,"Est. max error = %.7G\n",dev);
free_dvector(xs,1L);
free_dvector(fs,1L);
}
else /* Fit polynomial */
{ ncof=nn+5; xs=dvector(0L,(long)ncof,"xs in main");
fs=dvector(0L,(long)ncof,"fs in main");
inpt(ncof, npt, a, b, xln, yln, xs, fs, filenam);
chebpc(fs,cof,nn+1);
pcshft(a,b,cof,nn+1);
free_dvector(xs,0L); free_dvector(fs,0L);
}
for (i=0;i<=nn+nd;++i)
{ sprintf(nar,"%.*G ",dig,cof[i]);
sscanf(nar,"%lf",&cof[i]);
}
if (nd) printf("\n Rational function coefficien"
"ts\n Numerator Denominator\n\n");
else printf("\n Polynomial\n Coefficients\n\n");
for (i=0;i { if (i<=nn) printf("%-24.16G",cof[i]);
else printf("%*s",24," ");
if (i else printf("\n");
}
for(;;) /* Calculate trial values */
{ i=-1;
fprintf(stderr,
"Enter E to quit or a trial X value: ");
i=scanf("%lf",&x);
if (i<1) exit(1); if (xln) x=log(x);
yy=cof[nn]; for (j=nn-1;j>=0;j--) yy=yy*x+cof[j];
if (nd)
{ for (sumd=0.,j=nn+nd;j>=nn+1;j--)
sumd=(sumd+cof[j])*x;
yy=yy/(1.0+sumd);
}
if (yln) yy=exp(yy); if (xln) x=exp(x);
fprintf(stderr,"For X=%.8G Y=%.8G\n",x,yy);
}
}
/******************************************************
* COMMD - Parses the command line
*****************************************************/
void commd(int argc, char**argv, double *a,double *b,
int *nn, int *nd, int *dig, int *xln,
int *yln, char *filenam)
{ struct ffblk ffblk;
int i,j,k,l;
void help(char *msg);

*a=1.11e300; *b=-1.11e300; *nn=-32768; *nd=-32768;
*dig=7, *xln=0, *yln=0;
if (argc<2) help(""); filenam[0]=0;
for (i=1,k=0;i { for (j=0; j<(l=strlen(argv[i])); ++j)
{ argv[i][j]=toupper(argv[i][j]);
if (argv[i][j]=='=') ++k;
}
if (k)
{ if (k>1 || (!isdigit(argv[i][l-1])
&& argv[i][l-1] !='.'))
{ fprintf(stderr,"Invalid command line parameter"
" %s\n",argv[i]);
help("");
}
if (!strncmp(argv[i],"LL=",3))
{ k=sscanf(&argv[i][3],"%lf",a);
if (k<1) help("Invalid value for command line"
" parameter LL\n"); continue;
}
if (!strncmp(argv[i],"UL=",3))
{ k=sscanf(&argv[i][3],"%lf",b);
if (k<1) help("Invalid value for command line"
" parameter UL\n"); continue;
}
if (!strncmp(argv[i],"ND=",3))
{ k=sscanf(&argv[i][3],"%d",nd);
if (k<1) help("Invalid value for command line"
" parameter ND\n"); continue;
}
if (!strncmp(argv[i],"NN=",3))
{ k=sscanf(&argv[i][3],"%d",nn);
if (k<1) help("Invalid value for command line"
" parameter NN\n"); continue;
}
if (!strncmp(&argv[i][0],"-DIG=",5))
{ k=sscanf(&argv[i][6],"%d",dig);
if (k<1) help("Invalid value for optional"
" command line parameter DIG\n");
continue;
}
fprintf(stderr,"Invalid command line parameter "
"%s\n",argv[i]);
help("");
}
else
{ if (!strncmp(&argv[i][0],"-XLN",4))
{ *xln=1; continue; }
if (!strncmp(&argv[i][0],"-YLN",4))
{ *yln=1; continue; }
if (!filenam[0] && argv[i][0] != '-')
{ strcpy(filenam,argv[i]);
if (findfirst(filenam,&ffblk,0))
{ fprintf(stderr,"Input file %s doesn't "
"exist\n",filenam); help("");
}
}
else
{ fprintf(stderr,"Invalid command line parameter"
" %s\n",argv[i]); help("");
}
}
}
k=0;
if (*a>1.1e300)
{ fprintf(stderr,"Parameter LL undefined\n"); ++k; }
if (*b<-1.1e300)
{ fprintf(stderr,"Parameter UL undefined\n"); ++k; }
if (*nn==-32768)
{ fprintf(stderr,"Parameter NN undefined\n"); ++k;
*nn=5;
}
if (*nd==-32768)
{ fprintf(stderr,"Parameter ND undefined\n"); ++k;
*nd=5;
}
if (filenam[0]==0)
{ fprintf(stderr,"Input file is not defined\n");
++k;
}
if (*nn<1 || (*nd==0 && *nn>20) || (*nd && *nn>10))
{ fprintf(stderr,"Invalid value %d for "
"command line parameter NN\n",*nn); ++k;
}
if (*nd<0 || *nd>10)
{ fprintf(stderr,"Invalid value %d for "
"command line parameter ND\n",*nd); ++k;
}
if (*a>*b)
{ fprintf(stderr,"Lower limit > upper limit\n");
++k;
}
if (*dig<1 || *dig>20)
{ fprintf(stderr,"Invalid value %d for "
"command line option DIG\n",*dig); ++k;
}
if (*xln && *a<=0.)
{ fprintf(stderr,"Lower limit is <=0 with "
"logarithmic transform of X values\n"); ++k;
}
if (k) help("");
}
/******************************************************
* HELP - Prints the help screen
*****************************************************/
void help(char *msg)
{ char *hlp[] =
{ "USAGE: rpft LL=a UL=b NN=n ND=m [-DIG=n -XLN -YLN"
"] file",""," LL=a a is the lower limit of the"
"region for fit"," UL=b b is the upper limit o"
"f the region for fit",""," ND=m | m>0: rational"
" function, denominator degree m,"," NN=n | "
" numerator degree n. 1<=m<=10 1<=n<=10"," "
" | m=0: Polynomial degree n. 1<=n<=20",""," -DI"
"G=n (optional) n = number of significant digits"
," for coefficients on output. Defaul"
"t=7."," -XLN (optional) Transform X values to"
" ln(X)"," -YLN (optional) Transform Y values "
"to ln(Y)",""," file File with input X Y data "
"points, 1 per line","","All command line paramete"
"rs except DIG, XLN and YLN are ","required. They "
"may be in any order, upper or lower case."
};
int i,n=sizeof(hlp)/sizeof(char *)-1;

fprintf(stderr,"%s\n",msg);
for (i=0; i fprintf(stderr,"%s",hlp[n]); exit(1);
}
/******************************************************
* INPT - Read input data and calculate data for the
* Chebyshev polynomial or rational polynomial
* curve fit.
*****************************************************/
void inpt(int nn, int npt, double a, double b, int xln,
int yln,double *xs, double *fs, char filenam[])
{ double bma,bpa,hth,yy,*xa,*ya,*c,*d;
int i,j,n=0;
char nar[82];
FILE *infile;
void ratint(double xa[], double ya[], double *c,
double *d, int n, double x, double *y);

if ((infile=fopen(filenam,"r"))==NULL)
{ fprintf(stderr,"Can't open file\n"); exit(1); }
while (!feof(infile) && fgets(nar,80,infile))
{ i=sscanf(nar,"%lf%lf",&hth,&yy);
if (i<1) break;
++n;
if (i<2)
{ fprintf(stderr,"Need 2 numbers on line %d\n",n);
exit(1);
}
if (xln && hth<=0.)
{ fprintf(stderr,"Nonpositive X = %.8G on line %d "
"with ln(X) transformation\n",hth,n);
exit(1);
}
if (yln && yy<=0.)
{ fprintf(stderr,"Nonpositive Y = %.8G on line %d "
"with ln(Y) transformation\n",yy,n);
exit(1);
}
}
if (n<3)
{ fprintf(stderr,
"You provided only %d data points\n",n);
exit(1);
}
fseek(infile,0L,0);
xa=dvector(1L,(long)n,"xa in inpt");
ya=dvector(1L,(long)n,"ya in inpt");
c=dvector(1L,(long)n,"c in inpt");
d=dvector(1L,(long)n,"d in inpt");
for (i=1; i<=n; ++i)
{ (void)fgets(nar,80,infile);
sscanf(nar,"%lf%lf",&xa[i],&ya[i]);
if (xln) xa[i]=log(xa[i]);
if (yln) ya[i]=log(ya[i]);
}
fclose(infile);
if (!nn)
/* Calculate arrays of abcissa and function values
* for rational function evaluation in ratlsq.
* Return in xs and fs
*/
{ for (i=1;i<=npt;i++)
{ if (i < npt/2)
{ hth=PIO2*(i-1)/(npt-1.0);
yy=sin(hth); xs[i]=a+(b-a)*yy*yy;
}
else
{ hth=PIO2*(npt-i)/(npt-1.0);
yy=sin(hth); xs[i]=b-(b-a)*yy*yy;
}
ratint(xa, ya, c, d, n, xs[i], &yy);
fs[i]=yy;
}
}
else
/* Calculate Chebyshev coefficients and return in
* the array fs
*/
{ bma=0.5*(b-a); bpa=0.5*(b+a);
for (i=0;i { ratint(xa,ya,c,d,n,
cos(PI*(i+0.5)/nn)*bma+bpa,&xs[i]); }
for (i=0;i { yy=0.; for (j=0;j yy += xs[j]*cos(PI*i*(j+0.5)/nn);
fs[i]=2.0*yy/nn;
}
}
free_dvector(xa,1L); free_dvector(ya,1L);
free_dvector(c,1L); free_dvector(d,1L); return;
}
/******************************************************
* RATLSQ - Returns in cof[0..mm+kk] the coefficients
* of a rational function approximation to the X,Y data
* points xs[1..npt],fs[1..npt]. The deviation of the
* approximation (insofar as is known) returned as dev.
*****************************************************/
void ratlsq(double xs[],double fs[],int npt, int mm,
int kk, double cof[], double *dev)
{ void dsvbksb(double **u, double w[], double **v,
int m, int n, double b[], double x[]);
void dsvdcmp(double **a, int m, int n, double w[],
double **v);
int i,it,j,ncof;
double devmax,e=0.0,power,sum,sumn,sumd,*bb,*coff,
*ee,**u,**v,*w,*wt;

ncof=mm+kk+1;
bb=dvector(1L,(long)npt,"bb in ratlsq");
ee=dvector(1L,(long)npt,"ee in ratlsq");
coff=dvector(0L,(long)(ncof-1),"coff in ratlsq");
u=dmatrix(1L,(long)npt,1L,(long)ncof,"u in ratlsq");
v=dmatrix(1L,(long)ncof,1L,(long)ncof,"v in ratlsq");
w=dvector(1L,(long)ncof,"w in ratlsq");
wt=dvector(1L,(long)npt,"wt in ratlsq");
*dev=1.e50;
for (i=1;i<=npt;++i) { wt[i]=1.0; ee[i]=1.0; }
for (it=1;it<=5;it++)
{ for (i=1;i<=npt;i++)
{ power=wt[i]; bb[i]=power*(fs[i]+SIGN(e,ee[i]));
for (j=1;j<=mm+1;j++)
{ u[i][j]=power; power *= xs[i]; }
power = -bb[i];
for (j=mm+2;j<=ncof;j++)
{ power *= xs[i]; u[i][j]=power; }
}
dsvdcmp(u,npt,ncof,w,v);
dsvbksb(u,w,v,npt,ncof,bb,coff-1);
devmax=sum=0.0;
for (j=1;j<=npt;j++)
{ e=xs[j]; for (sumn=coff[mm],i=mm-1;i>=0;i--)
sumn=sumn*e+coff[i];
for (sumd=0.0,i=mm+kk;i>=mm+1;i--)
sumd=(sumd+coff[i])*e;
ee[j]=sumn/(1.0+sumd)-fs[j];
wt[j]=fabs(ee[j]); sum += wt[j];
if (wt[j] > devmax) devmax=wt[j];
}
e=sum/npt;
if (devmax <= *dev)
{ for(j=0;j fprintf(stderr,"For ratlsq iteration %d max error="
" %.5G\n",it,devmax);
}
free_dvector(wt,1L); free_dvector(w,1L);
free_dmatrix(v,1L,1L); free_dmatrix(u,1L,1L);
free_dvector(ee,1L); free_dvector(coff,0L);
free_dvector(bb,1L);
}
/******************************************************
* DSVDCMP - Computes singular value decomposition of
* the matrix a[1..m][1..n].
*****************************************************/
void dsvdcmp(double **a, int m, int n, double w[],
double **v)
{ double dpythag(double a, double b);
int flag,i,its,j,jj,k,l=0,nm=0;
double anorm,c,f,g,h,s,scale,x,y,z,*rv1;

rv1=dvector(1L,(long)n,"rv1 in dsvdcmp");
g=scale=anorm=0.0;
for (i=1;i<=n;i++)
{ l=i+1; rv1[i]=scale*g; g=s=scale=0.0;
if (i <= m)
{ for (k=i;k<=m;k++) scale += fabs(a[k][i]);
if (scale)
{ for (k=i;k<=m;k++)
{ a[k][i] /= scale; s += a[k][i]*a[k][i]; }
f=a[i][i]; g = -SIGN(sqrt(s),f); h=f*g-s;
a[i][i]=f-g;
for (j=l;j<=n;j++)
{ for(s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
f=s/h;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for (k=i;k<=m;k++) a[k][i] *= scale;
}
}
w[i]=scale*g; g=s=scale=0.0;
if (i <= m && i != n)
{ for (k=l;k<=n;k++) scale += fabs(a[i][k]);
if (scale)
{ for (k=l;k<=n;k++)
{ a[i][k] /= scale; s += a[i][k]*a[i][k]; }
f=a[i][l]; g = -SIGN(sqrt(s),f); h=f*g-s;
a[i][l]=f-g;
for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
for (j=l;j<=m;j++)
{ for(s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
}
for (k=l;k<=n;k++) a[i][k] *= scale;
}
}
x=fabs(w[i])+fabs(rv1[i]); if (x>anorm) anorm=x;
}
for (i=n;i>=1;i--)
{ if (i < n)
{ if (g)
{ for (j=l;j<=n;j++)
v[j][i]=(a[i][j]/a[i][l])/g;
for (j=l;j<=n;j++)
{ for(s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
}
}
for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
}
v[i][i]=1.0; g=rv1[i]; l=i;
}
for (i=min(m,n);i>=1;i--)
{ l=i+1; g=w[i];
for (j=l;j<=n;j++) a[i][j]=0.0;
if (g)
{ g=1.0/g;
for (j=l;j<=n;j++)
{ for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
f=(s/a[i][i])*g;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for (j=i;j<=m;j++) a[j][i] *= g;
}
else for (j=i;j<=m;j++) a[j][i]=0.0;
++a[i][i];
}
for (k=n;k>=1;k--)
{ for (its=1;its<=30;its++)
{ flag=1;
for (l=k;l>=1;l--)
{ nm=l-1;
if ((fabs(rv1[l])+anorm) == anorm)
{ flag=0; break; }
if ((fabs(w[nm])+anorm) == anorm) break;
}
if (flag)
{ c=0.0; s=1.0;
for (i=l;i<=k;i++)
{ f=s*rv1[i]; rv1[i]=c*rv1[i];
if ((fabs(f)+anorm) == anorm) break;
g=w[i]; h=dpythag(f,g); w[i]=h; h=1.0/h;
c=g*h; s = -f*h;
for (j=1;i<=m;j++)
{ y=a[j][nm]; z=a[j][i]; a[j][nm]=y*c+z*s;
a[j][i]=z*c-y*s;
}
}
}
z=w[k];
if (l == k)
{ if (z < 0.0)
{ w[k]=-z;
for (j=1;j<=n;j++) v[j][k] = -v[j][k];
}
break;
}
if (its == 30)
{ fprintf(stderr,"No convergence in 30 "
"dsvdcmp iterations\n");
exit(1);
}
x=w[l]; nm=k-1; y=w[nm]; g=rv1[nm]; h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=dpythag(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
c=s=1.0;
for (j=l;j<=nm;j++)
{ i=j+1; g=rv1[i]; y=w[i]; h=s*g; g=c*g;
z=dpythag(f,h); rv1[j]=z; c=f/z; s=h/z;
f=x*c+g*s; g=g*c-x*s; h=y*s; y *= c;
for (jj=1;jj<=n;jj++)
{ x=v[jj][j]; z=v[jj][i]; v[jj][j]=x*c+z*s;
v[jj][i]=z*c-x*s;
}
z=dpythag(f,h); w[j]=z;
if (z) { z=1.0/z; c=f*z; s=h*z; }
f=c*g+s*y; x=c*y-s*g;
for (jj=1;jj<=m;jj++)
{ y=a[jj][j]; z=a[jj][i]; a[jj][j]=y*c+z*s;
a[jj][i]=z*c-y*s;
}
}
rv1[l]=0.0; rv1[k]=f; w[k]=x;
}
}
free_dvector(rv1,1L);
}
/******************************************************
* DPYTHAG - Computes sqrt(a*a + b*b) without
* destructive underflow or overflow.
*****************************************************/
double dpythag(double a, double b)
{ double aba,abb,xx,yy;
aba=fabs(a); abb=fabs(b); xx=abb/aba; yy=aba/abb;
if (aba > abb) return aba*sqrt(1.0+xx*xx);
else return (abb == 0.0 ? 0.0 : abb*sqrt(1.0+yy*yy));
}
/******************************************************
* DSVBKSB - Solves A.X = B for a vector X, where A is
* specified by the arrays u[1..m][1..n], w[1..n],
* v[1..n][1..n] as returned by dsvdcmp.
*****************************************************/
void dsvbksb(double **u, double w[], double **v, int m,
int n, double b[], double x[])
{ int jj,j,i;
double s,*tmp;

tmp=dvector(1L,(long)n,"tmp in dsvbksb");
for (j=1;j<=n;j++)
{ s=0.0;
if (w[j])
{ for (i=1;i<=m;i++) s += u[i][j]*b[i]; s /= w[j];}
tmp[j]=s;
}
for (j=1;j<=n;j++)
{ s=0.0;
for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
x[j]=s;
}
free_dvector(tmp,1L);
}
/******************************************************
* CHEBPC - Transformation of Chebyshev coefficients
* generated in INPT to the interval -1 to 1
*****************************************************/
void chebpc(double c[],double d[],int n)
{ int k,j;
double sv,*dd;

dd=dvector(0L,(long)n-1,"dd in chebpc");
for (j=0;j d[0]=c[n-1];
for (j=n-2;j>=1;j--)
{ for (k=n-j;k>=1;k--)
{ sv=d[k]; d[k]=2.0*d[k-1]-dd[k]; dd[k]=sv; }
sv=d[0]; d[0] = -dd[0]+c[j]; dd[0]=sv;
}
for (j=n-1;j>=1;j--) d[j]=d[j-1]-dd[j];
d[0] = -dd[0]+0.5*c[0];
free_dvector(dd,0L);
}
/******************************************************
* PCSHFT - Generate the polynomial coefficients
* using the Chebyshev coefficients from CHEBPC
*****************************************************/
void pcshft(double a,double b,double d[],int n)
{ int k,j;
double fac,cnst;

cnst=2.0/(b-a); fac=cnst;
for (j=1;j cnst=0.5*(a+b);
for (j=0;j<=n-2;j++)
for (k=n-2;k>=j;k--) d[k] -= cnst*d[k+1];
}
/******************************************************
* RATINT - Diagonal rational function interpolation in
* the arrays xa[1..n] and ya[1..n].
*****************************************************/
void ratint(double xa[], double ya[], double *c,
double *d, int n, double x, double *y)
{ int m,i,ns=1;
double w,t,hh,h,dd;

hh=fabs(x-xa[1]);
for (i=1;i<=n;i++)
{ h=fabs(x-xa[i]);
if (h == 0.0) { *y=ya[i]; return; }
else if (h < hh) { ns=i; hh=h; }
c[i]=ya[i]; d[i]=ya[i]+1.e-50;
}
*y=ya[ns--];
for (m=1;m { for (i=1;i<=n-m;i++)
{ w=c[i+1]-d[i] ; h=xa[i+m]-x; t=(xa[i]-x)*d[i]/h;
dd=t-c[i+1];
if (dd == 0.0)
{ fprintf(stderr,"Error in routine ratint\n");
exit(1); /* The function may have a pole */
} /* at the requested value of x. */
dd=w/dd; d[i] =c [i+1]*dd; c[i]=t*dd;
}
*y += (2*ns < (n-m) ? c[ns+1] : d[ns--]);
}
return;
}
/******************************************************
* Utility routines based on Numerical Recipes
*****************************************************/
double *dvector(long nl, long nh, const char *ner)
{ double *v;

v=(double *)malloc((size_t) ((nh-nl+1+1)*
sizeof(double)));
if (!v)
{ fprintf(stderr,"Allocation failure for vector %s\n"
,ner); exit(1);
}
return v-nl+1;
}
double **dmatrix(long nrl,long nrh,long ncl,long nch,
const char *ner)
{ long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;

m=(double **) malloc((size_t)((nrow+1)*
sizeof(double*)));
if (m)
{ m += 1; m -= nrl;
m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*
sizeof(double)));
}
if (m[nrl])
{ m[nrl] += 1; m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
}
if (!m || !m[nrl])
{ fprintf(stderr,"Allocation failure for matrix %s\n"
,ner); exit(1);
}
return m;
}
void free_dvector(double *v, long nl)
{ free(v+nl-1); }

void free_dmatrix(double **m, long nrl, long ncl)
{ free(m[nrl]+ncl-1); free(m+nrl-1); }
/************* End of file RPFT.C *******************/


  3 Responses to “Category : Files from Magazines
Archive   : CUJ9306.ZIP
Filename : 1106032A

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/