/* 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 = '
'; form += '
'; form += ''; form += '
 Beispiele   '; if (opt == 1) form += ''; else form += ''; if (opt == 2) form += ''; else form += ''; if (opt == 3) form += ''; else form += ''; form += ''; form += '
'; form += '. Konvergenzstufe   '; form += '. Ableitung   '; form += '. Integral   '; form += ' Exponent   '; 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();