/* FFT form of Taylor and Fourier series and roots of polynomials 05.10.2022 */
/* global complex */
var canvas, centre = 314, colour, dif = 0, num = 32, opt = 0, twice = centre << 1;
function get_colour() {
// Get colour from cookie
// var ca = document.cookie.split(";"), ch, i, name = "bh_style=";
// for (i = 0; i < ca.length; i++) {
// ch = ca[i];
// while (ch.charAt(0) === " ") {
// ch = ch.substring(1);
// }
// if (ch.indexOf(name) === 0) {
// return (ch.substring(name.length, ch.length) === "bh_inblackandwhite") ? "#000000" : "#FFFFFF";
// }
// }
return "#000000"; // default colour
}
function form()
{
// Build and write form
var form = '';
document.writeln(form);
canvas = new Raphael("canvas", twice, twice + 1);
colour = get_colour();
}
function initialise_roots(size, root) {
// determine powers of roots
root[0] = new Complex(Decimal(1), Decimal(0));
root[1] = new Complex(Decimal(0), Decimal(-2).mul(Decimal.acos(-1)).div(Decimal(size))).cexp();
for (var k = 1; k < size / 2; k *= 2) {
for (var m = 1; m <= k; m++) root[k + m] = root[k].cmul(root[m]);
}
return;
}
function initialise_tft(signal, pbsize, realft, xvalue, method) {
var s = new Complex(Decimal(1), Decimal(0)), t = new Complex((method == 1) ? Decimal(0.5) : Decimal(1/16384), Decimal(0));
var u = new Complex(Decimal(0), Decimal.acos(-1).mul(Decimal(2).div(Decimal(pbsize)))).cexp(), x = xvalue.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])";
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);
x = xvalue.cadd(t.cmul(s));
signal[k] = eval(realft);
}
}
}
else
alert("Ausdruck enthält ungültige(s) Zeichen!");
return realft;
}
function initialise_fft(signal, pbsize, realft) {
var z = new Complex(Decimal.acos(-1).neg(), Decimal(0)), j, k, x, y;
var u = Decimal.acos(-1).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])";
realft = realft.replace(/\[/g, "Complex(Decimal(").replace(/\]/g, "),Decimal(0))");
try {
x = z;
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++) {
x = z.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, out, roo) {
var n = ins.length;
// base case
if (n == 1) out[0] = ins[0]; else {
// compute FFT of even terms
var even = [], evenFFT = [];
for (var k = 0; k < n / 2; k++) {
even[k] = ins[2 * k];
}
cfft(even, evenFFT, roo);
// compute FFT of odd terms
var odd = even, oddFFT = [];
for (var k = 0; k < n / 2; k++) {
odd[k] = ins[2 * k + 1];
}
cfft(odd, oddFFT, roo);
// combine
var l = n / 2, m = num / l;
for (var k = 0; k < l; k++) {
var uni = roo[k * m].cmul(oddFFT[k]);
out[k] = evenFFT[k].cadd(uni);
out[k + l] = evenFFT[k].csub(uni);
}
}
return;
}
function compute_taylor(result, pbsize, factor) {
var c = new Complex((pbsize >= 32) ? Decimal(0.5) : Decimal(1/16384), Decimal(0)), f = Decimal(1), s;
var h = new Complex(Decimal(pbsize), Decimal(0)), str = "Auf sechzehn Dezimalstellen gerundete Koeffizienten der Taylorreihe:\n";
if (factor < 0) for (var k = 1; k <= -factor; k++) f = f.div(Decimal(k));
for (var k = 0; k < 32 + ((factor > 0) ? factor : 0); k++) {
result[k] = result[k].cdiv(h);
h = h.cmul(c);
// integrate or differentiate
if (factor < 0) {
if (k > 0) f = f.mul(Decimal(k)).div(Decimal(k - factor));
}
else if (factor > 0) {
if (k > 0) f = f.mul(Decimal(k));
if (k > factor) f = f.div(Decimal(k - factor));
}
result[k].re = result[k].re.mul(Decimal(f.mul((factor > k) ? Decimal(0) : Decimal(1))));
s = (k < factor + 10) ? "a[0" : "a[";
if (k >= factor) {
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 */
}
}
return str;
}
function compute_fourier(result, pbsize, differ, factor) {
var f = new Complex(Decimal(differ), Decimal(0)), h = new Complex(Decimal(pbsize), Decimal(0)), n = 16, s;
var str = "Auf sechzehn Dezimalstellen gerundete Koeffizienten der Fourierreihe:\n";
// differentiate if necessary
if (differ > 0) result[0] = new Complex(Decimal(0), Decimal(0));
for (var j = -n; j < n; j++) {
k = (j < 0) ? j + pbsize : j;
result[k] = result[k].cdiv(h);
if (differ > 0 && k > 0) result[k] = result[k].cmul(Complex(Decimal(0), Decimal(k)).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, xvalue, conver, realft) {
var c0 = new Complex(Decimal( 0), Decimal(0)), c1 = new Complex(Decimal(1), Decimal(0 ));
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, z, r = Decimal(1e16 );
var str = "Auf sechzehn Dezimalstellen gerundete Iterationen:\n", e = xvalue, x = xvalue;
// 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 = z = c0;
var k = pbsize;
do {
f = f.cmul(v).cadd( fa[k]);
z = z.cmul(v).cadd(result[k]);
} while (--k > 0);
z = z.cmul(v).cadd(result[0]);
}
// precise and fast, quadratic, cubic, quartic, quintic, sextic, septic, and octic convergence
switch(conver) {
case 0:
y = x;
x = y.cadd(c9);
z = eval(realft);
x = y.csub(c9);
f = eval(realft);
x = c1.csub( z.cdiv(f));
x = c2.cdiv(x).csub(c1);
v = v.cadd( x.cmul(c9));
break;
case 1:
y = x;
x = y.cadd(c9);
z = eval(realft);
x = y.csub(c9);
f = eval(realft);
x = c1.csub( z.cdiv(f));
x = c2.cdiv(x).csub(c1);
w = x.cmul(c9);
z = z.cadd( f).cdiv(c2);
f = z.csub( f).cdiv(c9);
// see case 8
x = y.cadd(w);
g = eval(realft);
y = w.cmul(g);
w = z.csub(g.cmul(c2));
x = x.cadd(y.cdiv(w ));
h = eval(realft);
w = z.csub(g).cdiv( w);
y = h.cdiv( g.csub(h));
w = w.cmul(w).cadd( y);
y = h.cmul(c4);
y = y.cdiv( z.cadd(h));
y = y.cadd( w);
y = y.cmul( h.cdiv(f));
y = y.cadd(xvalue);
v = x.csub( y);
break;
case 2:
v = v.csub( z.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(z));
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(z);
x = f.cmul(f);
y = x.csub(g).cmul(c2);
g = g.cadd(y);
f = f.cmul( y.cdiv(g));
y = z.cdiv(x).cdiv(c6);
y = y.cmul(h);
y = f.cdiv(z).cadd(y );
v = v.csub(c1.cdiv(y));
break;
case 5:
x = v.csub( z.cdiv(f));
k = pbsize;
do {
g = g.cmul(x).cadd(result[k]);
h = h.cmul(x).cadd( fa[k]);
} while (--k > 0);
g = g.cmul(x).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));
z = f.cadd(c7.cmul(h));
y = y.cdiv(z).cmul(g );
v = x.csub(y );
break;
case 6:
x = v.csub( z.cdiv(f));
k = pbsize;
do {
g = g.cmul(x).cadd( fa[k]);
} while (--k > 0);
x = z.cdiv( f.cadd(g));
x = v.csub(c2.cmul(x));
k = pbsize;
do {
h = h.cmul(x).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 = x.csub(y);
break;
case 7:
y = z.cdiv(f);
x = v.csub(y);
k = pbsize;
do {
g = g.cmul(x).cadd(result[k]);
} while (k-- > 0);
y = y.cmul(g);
w = z.csub(g.cmul(c2));
x = x.csub(y.cdiv( w));
k = pbsize;
do {
h = h.cmul(x).cadd(result[k]);
} while (k-- > 0);
z = z.csub(g).cdiv( w);
y = h.cdiv( g.csub(h));
z = z.cmul(z).cadd( y);
z = z.cdiv(f).cmul( h);
v = x.csub(z);
break;
case 8:
y = z.cdiv(f);
x = v.csub(y);
k = pbsize;
do {
g = g.cmul(x).cadd(result[k]);
} while (k-- > 0);
y = y.cmul(g);
w = z.csub(g.cmul(c2));
x = x.csub(y.cdiv(w ));
k = pbsize;
do {
h = h.cmul(x).cadd(result[k]);
} while (k-- > 0);
w = z.csub(g).cdiv( w);
y = h.cdiv( g.csub(h));
w = w.cmul(w).cadd( y);
y = h.cmul(c4);
y = y.cdiv( z.cadd(h));
y = y.cadd( w);
y = y.cmul( h.cdiv(f));
v = x.csub( y);
break;
default :
}
// output if interesting
p = (m < 10) ? " " + m : m;
x = v.cadd(xvalue);
if (x.re.isNaN() || x.im.isNaN() || e.csub(x).cabs().re.lt(Decimal(2e-80))) {
m = 17;
str = str + "Berechnung vorzeitig beendet: Bitte eventuell andere Vorgabe wählen.";
}
else {
s = ((x.re.gte(Decimal(0))) ? "+" : "") + x.re.toExponential(16, 4);
str = str + p + ": Re(x) = " + s.padEnd(24) + " (" + ((x.re.gte(Decimal(0))) ? "+" : "") + x.re.toExponential(80, 4) + ")\n";
s = ((x.im.gte(Decimal(0))) ? "+" : "") + x.im.toExponential(16, 4);
str = str + "und Im(x) = " + s.padEnd(24) + " (" + ((x.im.gte(Decimal(0))) ? "+" : "") + x.im.toExponential(80, 4) + ")\n";
}
e = x;
}
return str;
}
function compute_coeffs() {
var aft, con, fac, fct, img, inp = [], n = Decimal(0), z = Decimal(0), oup = [], roo = [], index = document.fft.b.selectedIndex, rea, str = "", val, y = Decimal(1);
// Insert examples
if (index > 0)
{
switch (index)
{
case 1 :
document.fft.pbsize.value = "x.add([1]).log()";
break;
case 2 :
document.fft.pbsize.value = "x.atan()";
break;
case 3 :
document.fft.pbsize.value = "[1].div([1].sub(x))";
break;
case 4 :
document.fft.pbsize.value = "[200e-2].mul([I]).mul(x.exp())";
break;
case 5 :
document.fft.pbsize.value = "[PI].add([E].mul(x.mul(x).mul(x)))";
break;
case 6 :
document.fft.pbsize.value = "x.mul([2]).cos()";
break;
case 7 :
document.fft.pbsize.value = "x.abs()";
break;
case 8 :
document.fft.pbsize.value = "x.div(x.abs())";
break;
case 9 :
document.fft.pbsize.value = "x.sin().abs()";
break;
case 10 :
document.fft.pbsize.value = "[1].div(x).sin()";
break;
}
}
// determine Taylor or Fourier series, or compute zeros
fct = document.fft.pbsize.value.toString();
rea = (isNaN(document.fft.real.value)) ? 0 : Number(document.fft.real.value);
img = (isNaN(document.fft.imag.value)) ? 0 : Number(document.fft.imag.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.num.value) > 0) num = Math.pow(2, Number(document.fft.num.value));
if (opt == 1) {
val = new Complex(Decimal(rea), Decimal(img));
initialise_tft(inp, 2 * num, fct, val, opt);
if (inp[0] !== null) {
initialise_roots(2 * num, roo);
cfft(inp, oup, roo);
document.fft.result.value = compute_taylor(oup, num * 2, dif);
}
}
else if (opt == 2) {
num /= 2;
fac = initialise_fft(inp, num, fct);
if (inp[0] !== null) {
initialise_roots(2 * num, roo);
cfft(inp, oup, roo);
document.fft.result.value = compute_fourier(oup, num * 2, dif, fac);
}
num *= 2;
}
else if (opt == 3) {
num /= 2;
val = new Complex(Decimal(rea), Decimal(img));
con = Number(document.fft.conv.value);
fct = (con > 1) ? initialise_tft(inp, num, fct, val, opt) : initialise_tft(inp, 1, fct, val, opt);
if (inp[0] !== null) {
if (con > 1) {
initialise_roots(2 * num, roo);
cfft(inp, oup, roo);
}
document.fft.result.value = compute_simpson(oup, num, val, con, fct);
}
num *= 2;
}
show_fft(num, fct);
return;
}
function show_fft(asize, realf) {
var bnd = 1e5, fra = new Complex(Decimal(Math.PI / num), Decimal(0)), graph = [], inc, len, max = -bnd, min = bnd, k = -twice, r = asize - 0.5, quo = 2 * r / twice, tst;
var end_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])";
realf = realf.replace(/\[/g, "Complex(Decimal(").replace(/\]/g, "),Decimal(0))");
x = new Complex(Decimal(-Math.PI), Decimal(0));
// Determine minimum and maximum
for (var m = 0; m < 2 * num; m++) {
graph[m] = eval(realf);
tst = graph[m].re.toNumber();
if (tst > max) max = tst;
if (tst < min) min = tst;
x = x.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;
// Paint 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, "Re(f(x))").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});
// Paint function
inc = quo * k / 2 + r;
bnd = centre / bnd;
while (++k < twice) {
if (k % 2 == 0) {
start_y = centre - Math.round(bnd * graph[Math.round(inc)].re.toNumber());
inc += quo;
end_y = centre - Math.round(bnd * graph[Math.round(inc)].re.toNumber());
horizontal_axis = canvas.path("M" + (start_x - 1) + " " + start_y + "L" + ++start_x + " " + end_y);
horizontal_axis.attr({"stroke": colour, "stroke-width": 1 });
}
}
return;
}
form();
document.fft.result.value = "Beispiele:\n1) x.add([1]).log() für ln(x+1)\n2) x.atan() für arctan(x)\n3) [1].div([1].sub(x)) für 1/(1-x)\n4) [200e-2].mul([I]).mul(x.exp()) für 2ieˣ \n5) [PI].add([E].mul(x.mul(x).mul(x))) für π+ex³\n6) x.mul([2]).cos() für cos(2x)\n7) x.abs() für |x|\n8) x.div(x.abs()) für x/|x|\n9) x.sin().abs() für |sin(x)|\n10) [1].div(x).sin() für sin(1/x)";
compute_coeffs();