1 文本格式
using System;
namespace Legalsoft.Truffer
{
public class Ratfn
{
private double[] cofs { get; set; }
private int nn { get; set; }
private int dd { get; set; }
public Ratfn(double[] num, double[] den)
{
this.cofs = new double[num.Length + den.Length - 1];
this.nn = num.Length;
this.dd = den.Length;
for (int j = 0; j < nn; j++)
{
cofs[j] = num[j] / den[0];
}
for (int j = 1; j < dd; j++)
{
cofs[j + nn - 1] = den[j] / den[0];
}
}
public Ratfn(double[] coffs, int n, int d)
{
this.cofs = Globals.CopyFrom(coffs);
this.nn = n;
this.dd = d;
}
public double get(double x)
{
double sumn = 0.0;
double sumd = 0.0;
for (int j = nn - 1; j >= 0; j--)
{
sumn = sumn * x + cofs[j];
}
for (int j = nn + dd - 2; j >= nn; j--)
{
sumd = sumd * x + cofs[j];
}
return sumn / (1.0 + x * sumd);
}
public static Ratfn pade(double[] cof)
{
int n = (cof.Length - 1) / 2;
double[,] q = new double[n, n];
double[,] qlu = new double[n, n];
double[] x = new double[n];
double[] y = new double[n];
for (int j = 0; j < n; j++)
{
y[j] = cof[n + j + 1];
for (int k = 0; k < n; k++)
{
q[j, k] = cof[j - k + n];
}
}
LUdcmp lu = new LUdcmp(q);
lu.solve( y, x);
for (int j = 0; j < 4; j++)
{
lu.mprove(y, x);
}
for (int k = 0; k < n; k++)
{
double sum = cof[k + 1];
for (int j = 0; j <= k; j++)
{
sum -= x[j] * cof[k - j];
}
y[k] = sum;
}
double[] num = new double[n + 1];
double[] denom = new double[n + 1];
num[0] = cof[0];
denom[0] = 1.0;
for (int j = 0; j < n; j++)
{
num[j + 1] = y[j];
denom[j + 1] = -x[j];
}
return new Ratfn(num, denom);
}
public static Ratfn ratlsq(UniVarRealValueFun fn, double a, double b, int mm, int kk, ref double dev)
{
const int NPFAC = 8;
const int MAXIT = 5;
const double BIG = 1.0e99;
const double PIO2 = 1.570796326794896619;
int ncof = mm + kk + 1;
int npt = NPFAC * ncof;
double[] bb = new double[npt];
double[] coff = new double[ncof];
double[] ee = new double[npt];
double[] fs = new double[npt];
double[] wt = new double[npt];
double[] xs = new double[npt];
double[,] u = new double[npt, ncof];
Ratfn ratbest = new Ratfn(coff, mm + 1, kk + 1);
dev = BIG;
for (int i = 0; i < npt; i++)
{
if (i < (npt / 2) - 1)
{
double hth = PIO2 * i / (npt - 1.0);
xs[i] = a + (b - a) * Globals.SQR(Math.Sin(hth));
}
else
{
double hth = PIO2 * (npt - i) / (npt - 1.0);
xs[i] = b - (b - a) * Globals.SQR(Math.Sin(hth));
}
fs[i] = fn.funk(xs[i]);
wt[i] = 1.0;
ee[i] = 1.0;
}
double e = 0.0;
for (int it = 0; it < MAXIT; it++)
{
for (int i = 0; i < npt; i++)
{
double power = wt[i];
bb[i] = power * (fs[i] + Globals.SIGN(e, ee[i]));
for (int j = 0; j < mm + 1; j++)
{
u[i, j] = power;
power *= xs[i];
}
power = -bb[i];
for (int j = mm + 1; j < ncof; j++)
{
power *= xs[i];
u[i, j] = power;
}
}
SVD svd = new SVD(u);
svd.solve(bb, coff);
double devmax = 0.0;
double sum = 0.0;
Ratfn rat = new Ratfn(coff, mm + 1, kk + 1);
for (int j = 0; j < npt; j++)
{
ee[j] = rat.get(xs[j]) - fs[j];
wt[j] = Math.Abs(ee[j]);
sum += wt[j];
if (wt[j] > devmax)
{
devmax = wt[j];
}
}
e = sum / npt;
if (devmax <= dev)
{
ratbest = rat;
//ratbest.CopyFrom(rat);
dev = devmax;
}
//Console.Write(" ratlsq iteration= ");
//Console.Write(it);
//Console.Write(" max error= ");
//Console.Write("{0,10}", devmax);
//Console.Write("{0}", "\n");
}
return ratbest;
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
public class Ratfn
{
private double[] cofs { get; set; }
private int nn { get; set; }
private int dd { get; set; }
public Ratfn(double[] num, double[] den)
{
this.cofs = new double[num.Length + den.Length - 1];
this.nn = num.Length;
this.dd = den.Length;
for (int j = 0; j < nn; j++)
{
cofs[j] = num[j] / den[0];
}
for (int j = 1; j < dd; j++)
{
cofs[j + nn - 1] = den[j] / den[0];
}
}
public Ratfn(double[] coffs, int n, int d)
{
this.cofs = Globals.CopyFrom(coffs);
this.nn = n;
this.dd = d;
}
public double get(double x)
{
double sumn = 0.0;
double sumd = 0.0;
for (int j = nn - 1; j >= 0; j--)
{
sumn = sumn * x + cofs[j];
}
for (int j = nn + dd - 2; j >= nn; j--)
{
sumd = sumd * x + cofs[j];
}
return sumn / (1.0 + x * sumd);
}
public static Ratfn pade(double[] cof)
{
int n = (cof.Length - 1) / 2;
double[,] q = new double[n, n];
double[,] qlu = new double[n, n];
double[] x = new double[n];
double[] y = new double[n];
for (int j = 0; j < n; j++)
{
y[j] = cof[n + j + 1];
for (int k = 0; k < n; k++)
{
q[j, k] = cof[j - k + n];
}
}
LUdcmp lu = new LUdcmp(q);
lu.solve( y, x);
for (int j = 0; j < 4; j++)
{
lu.mprove(y, x);
}
for (int k = 0; k < n; k++)
{
double sum = cof[k + 1];
for (int j = 0; j <= k; j++)
{
sum -= x[j] * cof[k - j];
}
y[k] = sum;
}
double[] num = new double[n + 1];
double[] denom = new double[n + 1];
num[0] = cof[0];
denom[0] = 1.0;
for (int j = 0; j < n; j++)
{
num[j + 1] = y[j];
denom[j + 1] = -x[j];
}
return new Ratfn(num, denom);
}
public static Ratfn ratlsq(UniVarRealValueFun fn, double a, double b, int mm, int kk, ref double dev)
{
const int NPFAC = 8;
const int MAXIT = 5;
const double BIG = 1.0e99;
const double PIO2 = 1.570796326794896619;
int ncof = mm + kk + 1;
int npt = NPFAC * ncof;
double[] bb = new double[npt];
double[] coff = new double[ncof];
double[] ee = new double[npt];
double[] fs = new double[npt];
double[] wt = new double[npt];
double[] xs = new double[npt];
double[,] u = new double[npt, ncof];
Ratfn ratbest = new Ratfn(coff, mm + 1, kk + 1);
dev = BIG;
for (int i = 0; i < npt; i++)
{
if (i < (npt / 2) - 1)
{
double hth = PIO2 * i / (npt - 1.0);
xs[i] = a + (b - a) * Globals.SQR(Math.Sin(hth));
}
else
{
double hth = PIO2 * (npt - i) / (npt - 1.0);
xs[i] = b - (b - a) * Globals.SQR(Math.Sin(hth));
}
fs[i] = fn.funk(xs[i]);
wt[i] = 1.0;
ee[i] = 1.0;
}
double e = 0.0;
for (int it = 0; it < MAXIT; it++)
{
for (int i = 0; i < npt; i++)
{
double power = wt[i];
bb[i] = power * (fs[i] + Globals.SIGN(e, ee[i]));
for (int j = 0; j < mm + 1; j++)
{
u[i, j] = power;
power *= xs[i];
}
power = -bb[i];
for (int j = mm + 1; j < ncof; j++)
{
power *= xs[i];
u[i, j] = power;
}
}
SVD svd = new SVD(u);
svd.solve(bb, coff);
double devmax = 0.0;
double sum = 0.0;
Ratfn rat = new Ratfn(coff, mm + 1, kk + 1);
for (int j = 0; j < npt; j++)
{
ee[j] = rat.get(xs[j]) - fs[j];
wt[j] = Math.Abs(ee[j]);
sum += wt[j];
if (wt[j] > devmax)
{
devmax = wt[j];
}
}
e = sum / npt;
if (devmax <= dev)
{
ratbest = rat;
//ratbest.CopyFrom(rat);
dev = devmax;
}
//Console.Write(" ratlsq iteration= ");
//Console.Write(it);
//Console.Write(" max error= ");
//Console.Write("{0,10}", devmax);
//Console.Write("{0}", "\n");
}
return ratbest;
}
}
}
