/* FFT form of Taylor and Fourier series, Padé approximant, spline approximation and roots of polynomials 17.02.2024 */
/* global complex */
const CENTRE = 314, COLOUR = "#FFFFFF", LZ = Decimal.acos(-1), ZER = Complex(Decimal(0), Decimal(0)), ONE = Complex(Decimal(1), Decimal(0));
var canvas, com = 1, dif = 0, exp = 32, opt = 0, twice = CENTRE << 1;
function form()
{
// build and write form
var form = '';
document.writeln(form);
canvas = new Raphael("canvas", twice, twice + 1);
}
function solve_ls(mat) {
// solve linear equation of rank d with coefficients in m using Gauss Jordan algorithm
// m: array[d][d+1]
// returns solution in m[i][d+1]
// returns false if no solution found
//
// usage:
//
// EQ 1: a * x + b * y = c
// EQ 2: e * x + f * y = g
//
// var m = [ [ a, b, c ], [ e, f, g ] ];
// if (solve_ls(m)) {
// x = m[0][2];
// y = m[1][2];
// } else {
// no solution
// }
function find_max(rv) {
var ma = d, mi = Decimal(0);
for (var c = rv + 1; c < d; c++) {
if (mat[c][rv].cabs().re.gt(mi)) {
mi = mat[c][rv].cabs().re;
ma = c;
}
}
return ma;
}
function swap_rows(r1, r2) {
for (var c = 0; c <= d; c++) {
[mat[r1][c], mat[r2][c]] = [mat[r2][c], mat[r1][c]];
}
}
function norm_row(rp) {
// require m[r][r] != 0
mat[r][r] = ONE;
for (var c = r + 1; c <= d; c++) {
mat[r][c] = mat[r][c].cmul(rp);
}
}
function pivo_cell(za) {
mat[zr][r] = ZER;
for (var c = r + 1; c <= d; c++) {
mat[zr][c] = mat[zr][c].cadd(za.cmul(mat[r][c]));
}
}
var d = mat.length;
for (var r = 0; r < d; r++) {
if (mat[r][r].re.eq(Decimal(0)) && mat[r][r].im.eq(Decimal(0))) {
var nzr = find_max(r);
if (nzr == d) return false;
swap_rows(r, nzr);
}
// assert m[r][r] != 0
if (!(mat[r][r].re.eq(Decimal(0)) && mat[r][r].im.eq(Decimal(0)))) norm_row(ONE.cdiv(mat[r][r]));
for (var zr = 0; zr < d; zr++) {
if (zr != r && !(mat[zr][r].re.eq(Decimal(0)) && mat[zr][r].im.eq(Decimal(0)))) {
pivo_cell(mat[zr][r].cneg());
}
}
}
return true;
}
function initialise_tft(signal, pbsize, realft, zvalue, prefac, pshift, method, epsilo) {
var s = ONE, t = new Complex((method == 1) ? Decimal(epsilo) : Decimal(1/16384), Decimal(0));
var u = new Complex(Decimal(0), LZ.mul(Decimal(2).div(Decimal(pbsize)))).cexp(), z = zvalue.cadd(t);
// check input and compute coefficients
if (realft.search(/[^\.\(\)\w\-\[\]]/g) == -1) { /* check real function */
realft = realft.replace(/\./g, ".c").replace(/\[E\]/g, "Complex(Decimal(1).exp(),Decimal(0))");
realft = realft.replace(/\[I\]/g, "Complex(Decimal(0),Decimal(1))");
realft = realft.replace(/\[PI\]/g, "Complex(Decimal.acos(-1),Decimal(0))") + ".cadd([0])";
if (!(pshift.re.eq(Decimal(0)) && pshift.im.eq(Decimal(0)))) realft = realft.replace(/z/g, "(z.csub(Complex(Decimal(" + pshift.re.toNumber() + "), Decimal(" + + pshift.im.toNumber() + "))))");
if (!(prefac.re.eq(Decimal(1)) && prefac.im.eq(Decimal(0)))) realft = realft.replace(/z/g, "(z.cmul(Complex(Decimal(" + prefac.re.toNumber() + "), Decimal(" + + prefac.im.toNumber() + "))))");
realft = realft.replace(/\[/g, "Complex(Decimal(").replace(/\]/g, "),Decimal(0))");
try {
signal[0] = eval(realft);
} catch (exception) {
alert("Ausdruck ist nicht auswertbar!");
}
if (signal[0] !== null) {
for (var k = 1; k < pbsize; k++) {
s = s.cmul(u);
z = zvalue.cadd(t.cmul(s));
signal[k] = eval(realft);
}
}
}
else
alert("Ausdruck enthält ungültige(s) Zeichen!");
return realft;
}
function initialise_fft(signal, pbsize, prefac, pshift, realft) {
var x = new Complex(LZ.neg(), Decimal(0)), j, k, y, z;
var u = LZ.div(Decimal(pbsize)), e = Decimal(1e-17), f = 0;
// check input and compute coefficients
if (realft.search(/[^\.\(\)\w\-\[\]]/g) == -1) { /* check real function */
realft = realft.replace(/\./g, ".c").replace(/\[E\]/g, "Complex(Decimal(1).exp(),Decimal(0))");
realft = realft.replace(/\[I\]/g, "Complex(Decimal(0),Decimal(1))");
realft = realft.replace(/\[PI\]/g, "Complex(Decimal.acos(-1),Decimal(0))") + ".cadd([0])";
if (!(pshift.re.eq(Decimal(0)) && pshift.im.eq(Decimal(0)))) realft = realft.replace(/z/g, "(z.csub(Complex(Decimal(" + pshift.re.toNumber() + "), Decimal(" + + pshift.im.toNumber() + "))))");
if (!(prefac.re.eq(Decimal(1)) && prefac.im.eq(Decimal(0)))) realft = realft.replace(/z/g, "(z.cmul(Complex(Decimal(" + prefac.re.toNumber() + "), Decimal(" + + prefac.im.toNumber() + "))))");
realft = realft.replace(/\[/g, "Complex(Decimal(").replace(/\]/g, "),Decimal(0))");
try {
z = x;
signal[0] = eval(realft);
if (signal[0].re.abs().lt(e) || signal[0].re.isNaN()) signal[0].re = Decimal(0);
if (signal[0].im.abs().lt(e) || signal[0].im.isNaN()) signal[0].im = Decimal(0);
} catch (exception) {
alert("Ausdruck ist nicht auswertbar!");
}
if (signal[0] !== null) {
for (k = 1; k < 2 * pbsize; k++) {
z = x.cadd(Complex(u.mul(Decimal(k)), Decimal(0)));
signal[k] = eval(realft);
if (signal[k].re.abs().lt(e) || signal[k].re.isNaN()) signal[k].re = Decimal(0);
if (signal[k].im.abs().lt(e) || signal[k].im.isNaN()) signal[k].im = Decimal(0);
}
}
}
else
alert("Ausdruck enthält ungültige(s) Zeichen!");
return f;
}
function cfft(ins) {
var max = ins.length, i, j = 0, k, l = 1, m = max >> 1, n;
var one = Decimal(1), two = Decimal(2), c1 = ONE, cn = c1.cneg(), t, u;
// reverse bits
for (i = 0; i < max - 1; i++) {
if (i < j) [ins[i], ins[j]] = [ins[j], ins[i]];
k = m;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
// compute FFT
for (k = 0; k < Math.log2(max); k++) {
u = c1;
m = l;
l <<= 1;
for (j = 0; j < m; j++) {
for (i = j; i < max; i += l) {
n = i + m;
t = u.cmul(ins[n]);
ins[n] = ins[i].csub(t);
ins[i] = ins[i].cadd(t);
}
u = u.cmul(cn);
}
cn.im = one.sub(cn.re).div(two).sqrt().neg(); // do not negate for ifft, but normalise
cn.re = one.add(cn.re).div(two).sqrt();
}
}
function output_pade(result, factor, elldeg, emmdeg, numera, denomi) {
var d, s, str = "Auf sechzehn Dezimalstellen gerundete Koeffizienten der Taylorpolynome von Zähler p und Nenner q:\n";
var z = Decimal(0), matrix = [], l = elldeg, m = emmdeg, n = l - m;
if (factor > 0) for (var k = 0; k < 32; k++) result[k] = result[k + factor];
// build matrix
for (var i = 0; i < m; i++) {
matrix[i] = [];
matrix[i][m] = result[++n + m].cneg();
for (var j = 0; j < m; j++) matrix[i][j] = (n + j < 0) ? ZER : result[n + j];
}
// solve linear system and determine coefficients of quotient
if (solve_ls(matrix)) {
denomi[0] = ONE;
for (var j = 1; j <= m; j++) denomi[j] = matrix[m - j][m];
for (var i = 0; i <= l; i++) {
d = result[i];
for (var j = 1; j <= Math.min(i, m); j++) d = d.cadd(denomi[j].cmul(result[i - j]));
numera[i] = d;
}
// output quotient
for (var k = 0; k <= l; k++) {
s = (k < 10) ? "p[0" : "p[";
str = str + s + k + "] = ";
s = ((numera[k].re.gte(Decimal(0))) ? "+" : "") + numera[k].re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17); /* real coefficient */
s = ((numera[k].im.gte(Decimal(0))) ? " +" : " ") + numera[k].im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n"; /* complex coefficient */
}
for (var k = 0; k <= m; k++) {
s = (k < 10) ? "q[0" : "q[";
str = str + s + k + "] = ";
s = ((denomi[k].re.gte(Decimal(0))) ? "+" : "") + denomi[k].re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17); /* real coefficient */
s = ((denomi[k].im.gte(Decimal(0))) ? " +" : " ") + denomi[k].im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n"; /* complex coefficient */
}
}
else {
str = str + "Das lineare Gleichungssystem ist nicht lösbar: Bitte eventuell andere Vorgabe wählen.";
}
return str;
}
function output_taylor(result, factor, intnum, ulimit, llimit) {
var b = [], d, l = [], p = 32, s, str = "Auf sechzehn Dezimalstellen gerundete Koeffizienten der Taylorreihe:\n", v = ZER, w = ZER, u = [];
if (factor > 0) p += factor - 1;
b[0] = llimit;
d = ulimit.csub(llimit).cdiv(new Complex(Decimal(intnum), Decimal(0)));
for (var j = 0; j < intnum; j++) {
l[j] = u[j] = result[p];
b[j + 1] = b[j].cadd(d);
}
for (var j = 0; j < intnum; j++) {
for (var k = p - 1; k >= Math.max(0, factor); k--) {
u[j] = u[j].cmul(b[j + 1]).cadd(result[k]);
l[j] = l[j].cmul(b[j] ).cadd(result[k]);
}
}
for (var k = Math.max(0, factor); k <= p; k++) {
s = (k < factor + 10) ? "a[0" : "a[";
str = str + s + (k - factor) + "] = ";
s = ((result[k].re.gte(Decimal(0))) ? "+" : "") + result[k].re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17); /* real coefficient */
s = ((result[k].im.gte(Decimal(0))) ? " +" : " ") + result[k].im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n"; /* complex coefficient */
}
if (factor < 0) {
d = new Complex(Decimal(-factor), Decimal(0));
for (var j = 0; j < intnum; j++) {
u[j] = u[j].cmul(b[j + 1].cpow(d));
l[j] = l[j].cmul(b[j ].cpow(d));
}
}
for (var j = 0; j < intnum; j++) {
v = v.cadd(u[j]);
w = w.cadd(l[j]);
}
str = str + "Die Summe der oberen Werte beträgt:\n";
s = ((v.re.gte(Decimal(0))) ? "+" : "") + v.re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17);
s = ((v.im.gte(Decimal(0))) ? " +" : " ") + v.im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n";
str = str + "Die Summe der unteren Werte beträgt:\n";
s = ((w.re.gte(Decimal(0))) ? "+" : "") + w.re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17);
s = ((w.im.gte(Decimal(0))) ? " +" : " ") + w.im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n";
d = v.csub(w);
str = str + "Die Differenz der Werte beträgt:\n";
s = ((d.re.gte(Decimal(0))) ? "+" : "") + d.re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17);
s = ((d.im.gte(Decimal(0))) ? " +" : " ") + d.im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n";
return str;
}
function output_spline(spline, values, degree, knotno) {
var s, str = "\nAuf sechzehn Dezimalstellen gerundete Koeffizienten der Splines:";
for (var i = 0; i < knotno; i++) {
str = str + "\nKnoten " + (i + 1) + " bei ";
s = ((values[i].re.gte(Decimal(0))) ? "+" : "") + values[i].re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17).trimEnd() + ":\n"; /* real coefficient */
for (var j = 0; j <= degree; j++) {
str = str + "s[" + (i + 1) + "][" + j + "] = ";
s = ((spline[i][j].re.gte(Decimal(0))) ? "+" : "") + spline[i][j].re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17); /* real coefficient */
s = ((spline[i][j].im.gte(Decimal(0))) ? " +" : " ") + spline[i][j].im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n"; /* complex coefficient */
}
}
return str;
}
function compute_taylor(result, pbsize, factor, epsilo) {
var z = Decimal(0), p = 32, c = new Complex((pbsize >= p) ? Decimal(epsilo) : Decimal(1/16384), Decimal(0));
var f = Decimal(1), h = new Complex(Decimal(pbsize), Decimal(0));
if (factor < 0) for (var k = 1; k <= -factor; k++) f = f.div(Decimal(k)); else p += factor;
result[0] = result[0].cdiv(h);
h = h.cmul(c);
result[0] = result[0].cmul(Complex(f, z).cmul((factor > 0) ? ZER : ONE));
for (var k = 1; k < p; k++) {
result[k] = result[k].cdiv(h);
h = h.cmul(c);
// differentiate or integrate
if (factor < 0) {
f = f.mul(Decimal(k)).div(Decimal(k - factor));
}
else if (factor > 0) {
f = f.mul(Decimal(k));
if (k > factor) f = f.div(Decimal(k - factor));
}
result[k] = result[k].cmul(Complex(f, z).cmul((k < factor) ? ZER : ONE));
}
}
function compute_fourier(result, pbsize, differ) {
var f = new Complex(Decimal(differ), Decimal(0)), h = new Complex(Decimal(pbsize), Decimal(0)), n = pbsize / 2, s;
var str = "Auf sechzehn Dezimalstellen gerundete Koeffizienten der Fourierreihe:\n";
// differentiate if necessary
for (var j = -n; j < n; j++) {
k = (j < 0) ? j + pbsize : j;
result[k] = result[k].cdiv(h).cmul(Complex(Decimal(0), Decimal(Math.round(j))).cpow(f));
h = h.cneg();
s = (j < 0) ? "c[" : "c[ ";
if (Math.abs(j) < 10) s = s + " ";
str = str + s + j + "] = ";
s = ((result[k].re.gte(Decimal(0))) ? "+" : "") + result[k].re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17); /* real coefficient */
s = ((result[k].im.gte(Decimal(0))) ? " +" : " ") + result[k].im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n"; /* complex coefficient */
}
return str;
}
function compute_simpson(result, pbsize, zvalue, conver, realft) {
var c0 = ZER, c1 = ONE;
var c2 = new Complex(Decimal( 2), Decimal(0)), c3 = new Complex(Decimal(3), Decimal(0 ));
var c4 = new Complex(Decimal( 4), Decimal(0)), c5 = new Complex(Decimal(5), Decimal(0 ));
var c6 = new Complex(Decimal( 6), Decimal(0)), c7 = new Complex(Decimal(7), Decimal(0 ));
var c8 = new Complex(Decimal(16384), Decimal(0)), c9 = new Complex(Decimal(0), Decimal(1e-40));
var fa = [], ga = [], ha = [], height = [], f, g, h, p, v = c0, w, y, x, r = Decimal(1e16 );
var str = "Auf sechzehn Dezimalstellen gerundete Iterationen:\n", e = zvalue, z = zvalue;
// determine coefficients by differentiating Taylor series
if (conver > 1) {
height[0] = new Complex(Decimal(1 / pbsize), Decimal(0));
result[0] = result[0].cmul(height[0]);
for (var k = 1; k < pbsize; k++) {
height[k] = height[k - 1].cmul(c8);
result[k] = result[k].cmul(height[k]);
fa[k] = result[k].cmul(Complex(Decimal(k), Decimal(0)));
ga[k] = fa[k].cmul(Complex(Decimal(k - 1), Decimal(0)));
ha[k] = ga[k].cmul(Complex(Decimal(k - 2), Decimal(0)));
}
pbsize--;
}
// determine derivatives by using Horner's scheme
for (var m = 1; m < 17; m++) {
if (conver > 1) {
f = g = h = x = c0;
var k = pbsize;
do {
f = f.cmul(v).cadd( fa[k]);
x = x.cmul(v).cadd(result[k]);
} while (--k > 0);
x = x.cmul(v).cadd(result[0]);
}
// precise and fast, quadratic, cubic, quartic, quintic, sextic, septic, and octic convergence
switch(conver) {
case 0:
y = z;
z = y.cadd(c9);
x = eval(realft);
z = y.csub(c9);
f = eval(realft);
z = c1.csub( x.cdiv(f));
z = c2.cdiv(z).csub(c1);
v = v.cadd( z.cmul(c9));
break;
case 1:
y = z;
z = y.cadd(c9);
x = eval(realft);
z = y.csub(c9);
f = eval(realft);
z = c1.csub( x.cdiv(f));
z = c2.cdiv(z).csub(c1);
w = z.cmul(c9);
x = x.cadd( f).cdiv(c2);
f = x.csub( f).cdiv(c9);
// see case 8
z = y.cadd(w);
g = eval(realft);
y = w.cmul(g);
w = x.csub(g.cmul(c2));
z = z.cadd(y.cdiv(w ));
h = eval(realft);
w = x.csub(g).cdiv( w);
y = h.cdiv( g.csub(h));
w = w.cmul(w).cadd( y);
y = h.cmul(c4);
y = y.cdiv( x.cadd(h));
y = y.cadd( w);
y = y.cmul( h.cdiv(f));
y = y.cadd(zvalue);
v = z.csub( y);
break;
case 2:
v = v.csub( x.cdiv(f));
break;
case 3:
k = pbsize;
do {
g = g.cmul(v).cadd( ga[k]);
} while (--k > 1);
g = g.cdiv(c2.cmul(f));
g = g.csub( f.cdiv(x));
v = v.cadd(c1.cdiv(g));
break;
case 4:
k = pbsize;
do {
g = g.cmul(v).cadd( ga[k]);
h = h.cmul(v).cadd( ha[k]);
} while (--k > 2);
g = g.cmul(v).cadd( ga[2]);
g = g.cmul(x);
z = f.cmul(f);
y = z.csub(g).cmul(c2);
g = g.cadd(y);
f = f.cmul( y.cdiv(g));
y = x.cdiv(z).cdiv(c6);
y = y.cmul(h);
y = f.cdiv(x).cadd(y );
v = v.csub(c1.cdiv(y));
break;
case 5:
z = v.csub( x.cdiv(f));
k = pbsize;
do {
g = g.cmul(z).cadd(result[k]);
h = h.cmul(z).cadd( fa[k]);
} while (--k > 0);
g = g.cmul(z).cadd(result[0]);
h = h.cmul(h );
g = g.cdiv(f );
f = f.cmul(f );
y = f.cmul(c5);
y = y.cadd(c3.cmul(h));
x = f.cadd(c7.cmul(h));
y = y.cdiv(x).cmul(g );
v = z.csub(y );
break;
case 6:
z = v.csub( x.cdiv(f));
k = pbsize;
do {
g = g.cmul(z).cadd( fa[k]);
} while (--k > 0);
z = x.cdiv( f.cadd(g));
z = v.csub(c2.cmul(z));
k = pbsize;
do {
h = h.cmul(z).cadd(result[k]);
} while (k-- > 0);
y = c3.cmul(g).csub(f);
y = f.cadd(g).cdiv(y);
y = h.cdiv(f).cmul(y);
v = z.csub(y);
break;
case 7:
y = x.cdiv(f);
z = v.csub(y);
k = pbsize;
do {
g = g.cmul(z).cadd(result[k]);
} while (k-- > 0);
y = y.cmul(g);
w = x.csub(g.cmul(c2));
z = z.csub(y.cdiv( w));
k = pbsize;
do {
h = h.cmul(z).cadd(result[k]);
} while (k-- > 0);
x = x.csub(g).cdiv( w);
y = h.cdiv( g.csub(h));
x = x.cmul(x).cadd( y);
x = x.cdiv(f).cmul( h);
v = z.csub(x);
break;
case 8:
y = x.cdiv(f);
z = v.csub(y);
k = pbsize;
do {
g = g.cmul(z).cadd(result[k]);
} while (k-- > 0);
y = y.cmul(g);
w = x.csub(g.cmul(c2));
z = z.csub(y.cdiv(w ));
k = pbsize;
do {
h = h.cmul(z).cadd(result[k]);
} while (k-- > 0);
w = x.csub(g).cdiv( w);
y = h.cdiv( g.csub(h));
w = w.cmul(w).cadd( y);
y = h.cmul(c4);
y = y.cdiv( x.cadd(h));
y = y.cadd( w);
y = y.cmul( h.cdiv(f));
v = z.csub( y);
break;
default :
}
// output if interesting
p = (m < 10) ? " " + m : m;
z = v.cadd(zvalue);
if (z.re.isNaN() || z.im.isNaN() || e.csub(z).cabs().re.lt(Decimal(2e-80))) {
m = 17;
str = str + "Berechnung vorzeitig beendet: Bitte eventuell andere Vorgabe wählen.";
}
else {
s = ((z.re.gte(Decimal(0))) ? "+" : "") + z.re.toExponential(16, 4);
str = str + p + ": Re(z) = " + s.padEnd(24) + " (" + ((z.re.gte(Decimal(0))) ? "+" : "") + z.re.toExponential(80, 4) + ")\n";
s = ((z.im.gte(Decimal(0))) ? "+" : "") + z.im.toExponential(16, 4);
str = str + "und Im(z) = " + s.padEnd(24) + " (" + ((z.im.gte(Decimal(0))) ? "+" : "") + z.im.toExponential(80, 4) + ")\n";
}
e = z;
}
return str + "|" + e.re.toNumber();
}
function compute_coeffs() {
var deg, dif, eps, fac, fct, inp = [], inq = [], inr = [], kno = [], pee = [], pol = [], quu = [], lev, n = Decimal(0), num, x = n, index = document.fft.b.selectedIndex;
var lli, llr, lol, pfi, pfr, poi, por, shi, shr, str = "", uli, ulr, upl, vpf, vpo, vsh, y = Decimal(1);
// insert examples
if (index > 0)
{
switch (index)
{
case 1 :
document.fft.pbsize.value = "z.add([1]).log()";
break;
case 2 :
document.fft.pbsize.value = "z.atan()";
break;
case 3 :
document.fft.pbsize.value = "[1].div([1].sub(z))";
break;
case 4 :
document.fft.pbsize.value = "[200e-2].mul([I]).mul(z.exp())";
break;
case 5 :
document.fft.pbsize.value = "[PI].add([E].mul(z.mul(z).mul(z)))";
break;
case 6 :
document.fft.pbsize.value = "z.mul([2]).cos()";
break;
case 7 :
document.fft.pbsize.value = "z.abs()";
break;
case 8 :
document.fft.pbsize.value = "z.div(z.abs())";
break;
case 9 :
document.fft.pbsize.value = "z.sin().abs()";
break;
case 10 :
document.fft.pbsize.value = "[1].div(z).sin()";
break;
}
}
// determine Taylor or Fourier series, compute Padé approximant or spline approximation and zeros resp.
fct = document.fft.pbsize.value.toString();
eps = 1 / Math.pow(2, Number(document.fft.eps.value));
por = (isNaN(document.fft.pore.value)) ? 0 : Number(document.fft.pore.value);
poi = (isNaN(document.fft.poim.value)) ? 0 : Number(document.fft.poim.value);
ulr = (isNaN(document.fft.ulre.value)) ? 0 : Number(document.fft.ulre.value);
uli = (isNaN(document.fft.ulim.value)) ? 0 : Number(document.fft.ulim.value);
llr = (isNaN(document.fft.llre.value)) ? 0 : Number(document.fft.llre.value);
lli = (isNaN(document.fft.llim.value)) ? 0 : Number(document.fft.llim.value);
pfr = (isNaN(document.fft.pfre.value)) ? 0 : Number(document.fft.pfre.value);
pfi = (isNaN(document.fft.pfim.value)) ? 0 : Number(document.fft.pfim.value);
shr = (isNaN(document.fft.shre.value)) ? 0 : Number(document.fft.shre.value);
shi = (isNaN(document.fft.shim.value)) ? 0 : Number(document.fft.shim.value);
vpo = new Complex(Decimal(por), Decimal(poi));
upl = new Complex(Decimal(ulr), Decimal(uli));
lol = new Complex(Decimal(llr), Decimal(lli));
vpf = new Complex(Decimal(pfr), Decimal(pfi));
vsh = new Complex(Decimal(shr), Decimal(shi));
if (Number(document.fft.inn.value) > 0) inn = Number(document.fft.inn.value);
if (Number(document.fft.com.value) > 0) com = Number(document.fft.com.value);
if (Number(document.fft.opt.value) > 0) opt = Number(document.fft.opt.value);
if (Number(document.fft.dif.value) >= 0) dif = Number(document.fft.dif.value);
if (Number(document.fft.ite.value) >= 0 && dif == 0 && opt != 2) {
dif = -Number(document.fft.ite.value);
document.fft.dif.value = 0;
}
else {
document.fft.ite.value = 0;
}
if (Number(document.fft.exp.value) > 0) exp = Math.pow(2, Number(document.fft.exp.value));
// process options
switch (opt)
{
case 1 :
initialise_tft(inp, exp * 2, fct, vpo, vpf, vsh, opt, eps);
if (inp[0] !== null) {
cfft(inp);
compute_taylor(inp, exp * 2, dif, eps);
document.fft.result.value = output_taylor(inp, dif, inn, upl, lol);
pee = inp;
ell = exp;
quu[0] = ONE;
emm = 0;
}
break;
case 2 :
exp /= 2;
initialise_fft(inp, exp, vpf, vsh, fct);
if (inp[0] !== null) {
cfft(inp);
document.fft.result.value = compute_fourier(inp, exp * 2, 0);
pee = inp;
ell = exp;
quu[0] = ONE;
emm = 0;
}
exp /= 2;
break;
case 3 :
exp /= 2;
lev = Number(document.fft.lev.value);
fac = (lev > 1) ? initialise_tft(inp, exp, fct, vpo, vpf, vsh, opt, eps) : initialise_tft(inp, 1, fct, vpo, vpf, vsh, opt, eps);
if (inp[0] !== null) {
if (lev > 1) cfft(inp);
quu[0] = compute_simpson(inp, exp, vpo, lev, fac);
pee = quu[0].split("|");
document.fft.result.value = pee[0];
quu[0] = pee[1];
pee = inp;
ell = exp;
emm = 0;
}
exp *= 2;
break;
case 4 :
initialise_tft(inp, exp * 2, fct, vpo, vpf, vsh, 1, eps);
if (inp[0] !== null) {
cfft(inp);
compute_taylor(inp, exp * 2, 0, eps);
ell = Number(document.fft.ell.value);
emm = Number(document.fft.emm.value);
document.fft.result.value = output_pade(inp, 0, ell, emm, pee, quu);
}
break;
case 5 :
initialise_tft(inp, exp * 2, fct, vpo, vpf, vsh, 1, eps);
if (inp[0] !== null) {
deg = Number(document.fft.deg.value);
num = Number(document.fft.num.value);
cfft(inp);
compute_taylor(inp, exp * 2, 0, eps);
kno[0] = new Complex(LZ.neg(), Decimal(0));
lev = new Complex(LZ.mul(Decimal(2).div(Decimal(num - 1))), Decimal(0));
for (var i = 0; i < num; i++) {
pol[i] = [];
initialise_tft(pol[i], exp * 2, fct, kno[i], vpf, vsh, 1, eps);
cfft(pol[i]);
compute_taylor(pol[i], exp * 2, 0, eps);
kno[i+1] = kno[i].cadd(lev);
}
document.fft.result.value = output_taylor(inp, 0, inn, upl, lol) + output_spline(pol, kno, deg, num);
pee = inp;
quu[0] = ONE;
ell = num;
emm = deg;
}
break;
}
if (inp[0] === null) {
document.fft.result.value = str + "Auswertung ist nicht möglich.";
}
// draw graph
show_fft(exp, fct, com, opt, ell, emm, pee, quu, vpo, vpf, vsh, dif, kno, pol);
}
function show_fft(asize, realf, compl, optio, eldeg, emdeg, numer, denom, value, prefa, pshif, diffe, knots, polys) {
var bnd = 1e5, fac = new Complex(Decimal(-diffe), Decimal(0)), fra = new Complex(LZ.div(Decimal(asize)), Decimal(0)), graph = [], padeg = [], den, hlp, inc, len, max = -bnd, min = bnd, meths =[" ", "Taylor", "Fourier", "Nullstelle", "Padé", "Spline"];
var k = -twice, m = 2 * asize, r = asize - 0.5, z = Complex(LZ, Decimal(0)).cneg(), quo = 2 * r / twice, str, tst, zend_x, end_y, horizontal_axis, vertical_axis, start_x, start_y;
realf = realf.replace(/\./g, ".c").replace(/\[E\]/g, "Complex(Decimal(1).exp(),Decimal(0))");
realf = realf.replace(/\[I\]/g, "Complex(Decimal(0),Decimal(1))");
realf = realf.replace(/\[PI\]/g, "Complex(Decimal.acos(-1),Decimal(0))") + ".cadd([0])";
if (!(pshif.re.eq(Decimal(0)) && pshif.im.eq(Decimal(0)))) realf = realf.replace(/z/g, "(z.csub(Complex(Decimal(" + pshif.re.toNumber() + "), Decimal(" + + pshif.im.toNumber() + "))))");
if (!(prefa.re.eq(Decimal(1)) && prefa.im.eq(Decimal(0)))) realf = realf.replace(/z/g, "(z.cmul(Complex(Decimal(" + prefa.re.toNumber() + "), Decimal(" + + prefa.im.toNumber() + "))))");
realf = realf.replace(/\[/g, "Complex(Decimal(").replace(/\]/g, "),Decimal(0))");
switch (optio)
{
// compute Taylor series if requested
case 1:
z = z.csub(value);
inc = (diffe > 0) ? diffe : 0;
for (var n = 0; n < 2 * asize; n++) {
hlp = numer[eldeg];
for (var l = eldeg + inc - 1; l >= inc; l--) hlp = hlp.cmul(z).cadd(numer[l]);
padeg[n] = (diffe < 0) ? hlp.cmul(z.cpow(fac)) : hlp;
z = z.cadd(fra);
}
break;
// compute Fourier series if requested
case 2:
for (var n = 0; n < 2 * asize; n++) {
hlp = numer[0];
for (var l = 1 - m; l < m; l++) hlp = hlp.cadd(z.cmul(Complex(Decimal(0), Decimal(-l))).cexp().cmul(numer[m - l]));
padeg[n] = hlp;
z = z.cadd(fra);
}
break;
// compute Padé approximant if requested
case 4:
z = z.csub(value);
for (var n = 0; n < 2 * asize; n++) {
hlp = numer[eldeg];
for (var l = eldeg - 1; l >= 0; l--) hlp = hlp.cmul(z).cadd(numer[l]);
den = denom[emdeg];
for (var m = emdeg - 1; m >= 0; m--) den = den.cmul(z).cadd(denom[m]);
padeg[n] = hlp.cdiv(den);
z = z.cadd(fra);
}
break;
// compute spline approximation if requested
case 5:
fac = new Complex(LZ.mul(Decimal(2).div(Decimal(eldeg - 1))), Decimal(0));
eldeg = 0;
for (var n = 0; n < 2 * asize; n++) {
hlp = polys[eldeg][emdeg];
z = z.csub(knots[eldeg]);
for (var m = emdeg - 1; m >= 0; m--) hlp = hlp.cmul(z).cadd(polys[eldeg][m]);
z = z.cadd(knots[eldeg]);
padeg[n] = hlp.cmul(knots[++eldeg].csub(z));
hlp = polys[eldeg][emdeg];
z = z.csub(knots[eldeg]);
for (var m = emdeg - 1; m >= 0; m--) hlp = hlp.cmul(z).cadd(polys[eldeg][m]);
z = z.cadd(knots[eldeg]);
padeg[n] = padeg[n].cadd(hlp.cmul(z.csub(knots[--eldeg]))).cdiv(fac);
z = z.cadd(fra);
if (z.re.toNumber() > knots[eldeg+1].re.toNumber()) eldeg++;
}
break;
}
// determine minimum and maximum
z = new Complex(LZ.neg(), Decimal(0));
for (var n = 0; n < 2 * asize; n++) {
graph[n] = eval(realf);
tst = (compl == 1) ? graph[n].re.toNumber() : graph[n].im.toNumber();
if (tst > max) max = tst;
if (tst < min) min = tst;
z = z.cadd(fra);
}
if (max > bnd) max = bnd;
if (min < -bnd) min = -bnd;
bnd = Math.round(Math.max(Math.abs(max), Math.abs(min)) + 0.499);
len = bnd.toString().length;
// draw axes
canvas.clear();
start_x = 0;
start_y = 0;
end_x = twice;
end_y = twice;
horizontal_axis = canvas.path("M" + start_x + " " + CENTRE + "L" + end_x + " " + CENTRE);
horizontal_axis.attr({"stroke": COLOUR, "stroke-width": 1 });
vertical_axis = canvas.path("M" + CENTRE + " " + start_y + "L" + CENTRE + " " + end_y);
vertical_axis.attr({ "stroke": COLOUR, "stroke-width": 1 });
// add texts
canvas.text(20, 6, ((compl == 1) ? "Re" : "Im") + "(f(z))").attr({fill: COLOUR});
canvas.text(6, CENTRE + 6, "-π").attr({fill: COLOUR});
canvas.text(twice - 6, CENTRE + 6, "π").attr({fill: COLOUR});
canvas.text(CENTRE - 4 * len++, 6, bnd.toString()).attr({fill: COLOUR});
canvas.text(CENTRE - 4 * len, twice - 6, "-" + bnd.toString()).attr({fill: COLOUR});
canvas.text((optio == 3) ? 30 : 20, 20, ((optio == 3) ? "X " : "--- ") + meths[optio]).attr({fill: "#FF0000"});
// draw function
inc = quo * k / 2 + r;
bnd = CENTRE / bnd;
while (++k < twice) {
if (k % 2 == 0) {
start_y = (compl == 1) ? CENTRE - Math.round(bnd * graph[Math.round(inc)].re.toNumber()) : CENTRE - Math.round(bnd * graph[Math.round(inc)].im.toNumber());
inc += quo;
end_y = (compl == 1) ? CENTRE - Math.round(bnd * graph[Math.round(inc)].re.toNumber()) : CENTRE - Math.round(bnd * graph[Math.round(inc)].im.toNumber());
horizontal_axis = canvas.path("M" + (start_x - 1) + " " + start_y + "L" + ++start_x + " " + end_y);
horizontal_axis.attr({"stroke": COLOUR, "stroke-width": 1 });
}
}
// draw second result in red if available
if (optio == 3) {
canvas.text(CENTRE + Math.round(CENTRE * denom[0] / Math.PI), CENTRE, "X").attr({fill: "#FF0000"});
}
else {
start_x = 0;
k = -twice;
inc = quo * k / 2 + r;
while (++k < twice) {
if (k % 2 == 0) {
start_y = (compl == 1) ? CENTRE - Math.round(bnd * padeg[Math.round(inc)].re.toNumber()) : CENTRE - Math.round(bnd * padeg[Math.round(inc)].im.toNumber());
inc += quo;
end_y = (compl == 1) ? CENTRE - Math.round(bnd * padeg[Math.round(inc)].re.toNumber()) : CENTRE - Math.round(bnd * padeg[Math.round(inc)].im.toNumber());
horizontal_axis = canvas.path("M" + (start_x - 1) + " " + start_y + "L" + ++start_x + " " + end_y);
horizontal_axis.attr({"stroke": "#FF0000", "stroke-width": 1 });
}
}
}
}
// main programme
form();
document.fft.result.value = "Beispiele:\n1) z.add([1]).log() für ln(z+1)\n2) z.atan() für arctan(z)\n3) [1].div([1].sub(z)) für 1/(1-z)\n4) [200e-2].mul([I]).mul(z.exp()) für 2ieᶻ \n5) [PI].add([E].mul(z.mul(z).mul(z))) für π+ez³\n6) z.mul([2]).cos() für cos(2z)\n7) z.abs() für |z|\n8) z.div(z.abs()) für z/|z|\n9) z.sin().abs() für |sin(z)|\n10) [1].div(z).sin() für sin(1/z)";
compute_coeffs();