/* 28.05.2026: FFT form of Fourier and (reverse) Taylor series, Padé approximant, spline approximation and roots of polynomials by using 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;
// ------------------------------------------------------------
// Parser: "gewohnte" Eingabe -> Complex-Auswertung ohne umständliche Syntax
// Unterstützt: + - * / ^ (auch **), Klammern, implizite Multiplikation,
// Funktionen: sin, cos, tan, exp, log/ln, sqrt, abs, atan, ...
// Konstanten: pi, e, i
// Variable: z
// ------------------------------------------------------------
const COMPLEX_CONST = {
pi: "Complex(Decimal.acos(-1), Decimal(0))",
e: "Complex(Decimal(1).exp(), Decimal(0))",
i: "Complex(Decimal(0), Decimal(1))"
};
// Funktionsname (User) -> Complex-Prototyp-Methode (bh_bigcomplex.js)
const FUNC_MAP = {
sin: "csin", cos: "ccos", tan: "ctan",
cot: "ccot", sec: "csec", csc: "ccsc",
asin: "casin", acos: "cacos", atan: "catan",
acot: "cacot", asec: "casec", acsc: "cacsc",
sinh: "csinh", cosh: "ccosh", tanh: "ctanh",
coth: "ccoth", sech: "csech", csch: "ccsch",
asinh: "casinh", acosh: "cacosh", atanh: "catanh",
acoth: "cacoth", asech: "casech", acsch: "cacsch",
exp: "cexp",
log: "clog",
ln: "clog",
sqrt: "csqrt",
abs: "cabs"
};
// Operatoren: Präzedenz/Assoziativität
const OP_INFO = {
"+": { prec: 10, assoc: "L", arity: 2 },
"-": { prec: 10, assoc: "L", arity: 2 },
"*": { prec: 20, assoc: "L", arity: 2 },
"/": { prec: 20, assoc: "L", arity: 2 },
"^": { prec: 30, assoc: "R", arity: 2 },
"u+": { prec: 25, assoc: "R", arity: 1 },
"u-": { prec: 25, assoc: "R", arity: 1 }
};
function compileMathExprToJS(input) {
const tokens0 = tokenize(input);
const tokens1 = insertImplicitMultiplication(tokens0);
const rpn = toRPN(tokens1);
return compileRPNToJS(rpn); // liefert JS-Expr als String, Variable: zz (Complex)
}
function tokenize(s) {
const tokens = [];
let i = 0;
const isWS = c => /\s/.test(c);
const isDigit = c => /[0-9]/.test(c);
const isIdentStart = c => /[A-Za-z_]/.test(c);
const isIdent = c => /[A-Za-z0-9_]/.test(c);
while (i < s.length) {
const c = s[i];
if (isWS(c)) { i++; continue; }
if ("+-*/^(),".includes(c)) {
if (c === "*" && s[i + 1] === "*") { tokens.push({ type: "op", value: "^" }); i += 2; continue; }
if (c === "(") tokens.push({ type: "lparen", value: "(" });
else if (c === ")") tokens.push({ type: "rparen", value: ")" });
else if (c === ",") tokens.push({ type: "comma", value: "," });
else tokens.push({ type: "op", value: c });
i++;
continue;
}
// Zahl: 12, 12.3, .5, 1e-3, 1.2E+4
if (isDigit(c) || (c === "." && isDigit(s[i + 1] || ""))) {
let j = i;
while (isDigit(s[j] || "")) j++;
if (s[j] === ".") {
j++;
while (isDigit(s[j] || "")) j++;
}
if (s[j] === "e" || s[j] === "E") {
let k = j + 1;
if (s[k] === "+" || s[k] === "-") k++;
if (!isDigit(s[k] || "")) throw new Error("Ungültiger Exponent in Zahl bei Position " + j);
while (isDigit(s[k] || "")) k++;
j = k;
}
tokens.push({ type: "number", value: s.slice(i, j) });
i = j;
continue;
}
// Identifier
if (isIdentStart(c)) {
let j = i + 1;
while (isIdent(s[j] || "")) j++;
const name = s.slice(i, j);
tokens.push({ type: "ident", value: name });
i = j;
continue;
}
throw new Error("Ungültiges Zeichen '" + c + "' bei Position " + i);
}
return tokens;
}
function isFunctionName(name) {
return Object.prototype.hasOwnProperty.call(FUNC_MAP, name.toLowerCase());
}
function insertImplicitMultiplication(tokens) {
const out = [];
const isValue = (t) => t.type === "number" || t.type === "ident" || t.type === "rparen";
const isValueStart = (t, prev) => {
if (!t) return false;
// ident '(' ist Funktionsaufruf => keine implizite Multiplikation
if (t.type === "lparen" && prev && prev.type === "ident" && isFunctionName(prev.value)) return false;
return t.type === "number" || t.type === "ident" || t.type === "lparen";
};
for (let i = 0; i < tokens.length; i++) {
const t = tokens[i];
const prev = out[out.length - 1];
if (prev && isValue(prev) && isValueStart(t, prev)) out.push({ type: "op", value: "*" });
out.push(t);
}
return out;
}
function toRPN(tokens) {
const output = [];
const stack = [];
const argCountStack = [];
const peek = () => stack[stack.length - 1];
let prev = null;
for (let i = 0; i < tokens.length; i++) {
const t = tokens[i];
if (t.type === "number") { output.push(t); prev = t; continue; }
if (t.type === "ident") {
const name = t.value.toLowerCase();
const next = tokens[i + 1];
if (next && next.type === "lparen" && isFunctionName(name)) {
stack.push({ type: "func", value: name });
prev = t;
} else {
output.push({ type: "ident", value: name });
prev = t;
}
continue;
}
if (t.type === "comma") {
while (stack.length && peek().type !== "lparen") output.push(stack.pop());
if (!stack.length) throw new Error("Komma ohne passende Klammer.");
if (!argCountStack.length) throw new Error("Komma außerhalb eines Funktionsaufrufs.");
argCountStack[argCountStack.length - 1] += 1;
prev = t;
continue;
}
if (t.type === "lparen") {
stack.push(t);
const before = stack[stack.length - 2];
if (before && before.type === "func") argCountStack.push(0);
prev = t;
continue;
}
if (t.type === "rparen") {
while (stack.length && peek().type !== "lparen") output.push(stack.pop());
if (!stack.length) throw new Error("Schließende Klammer ohne öffnende Klammer.");
stack.pop(); // '('
if (stack.length && peek().type === "func") {
const fn = stack.pop();
const commas = argCountStack.pop() ?? 0;
output.push({ type: "func", value: fn.value, argc: commas + 1 });
}
prev = t;
continue;
}
if (t.type === "op") {
let op = t.value;
const prevIsOpOrStart = !prev || prev.type === "op" || prev.type === "lparen" || prev.type === "comma";
if ((op === "+" || op === "-") && prevIsOpOrStart) op = (op === "+") ? "u+" : "u-";
const o1 = OP_INFO[op];
if (!o1) throw new Error("Unbekannter Operator: " + op);
while (stack.length) {
const top = peek();
if (top.type !== "op") break;
const o2 = OP_INFO[top.value];
if (!o2) break;
// Prefix unary operators must not pop a preceding binary operator.
// This keeps 2^-3 valid, while -z^2 becomes -(z^2).
if (o1.arity === 1) break;
const cond =
(o1.assoc === "L" && o1.prec <= o2.prec) ||
(o1.assoc === "R" && o1.prec < o2.prec);
if (!cond) break;
output.push(stack.pop());
}
stack.push({ type: "op", value: op });
prev = t;
continue;
}
throw new Error("Interner Parserfehler bei Token: " + JSON.stringify(t));
}
while (stack.length) {
const t = stack.pop();
if (t.type === "lparen" || t.type === "rparen") throw new Error("Nicht geschlossene Klammer.");
output.push(t);
}
return output;
}
function compileRPNToJS(rpn) {
const st = [];
const numToJS = (numStr) => `Complex(Decimal("${numStr}"), Decimal(0))`;
const identToJS = (name) => {
if (name === "z") return "zz";
if (COMPLEX_CONST[name]) return COMPLEX_CONST[name];
throw new Error("Unbekannte Variable/Konstante: " + name + " (erlaubt: z, pi, e, i)");
};
for (const t of rpn) {
if (t.type === "number") { st.push(numToJS(t.value)); continue; }
if (t.type === "ident") { st.push(identToJS(t.value)); continue; }
if (t.type === "op") {
const info = OP_INFO[t.value];
if (info.arity === 1) {
const a = st.pop();
if (a === undefined) throw new Error("Operator " + t.value + ": Stack unterlaufen.");
if (t.value === "u+") st.push(`(${a})`);
else st.push(`(${a}).cneg()`);
} else {
const b = st.pop(), a = st.pop();
if (a === undefined || b === undefined) throw new Error("Operator " + t.value + ": Stack unterlaufen.");
switch (t.value) {
case "+": st.push(`(${a}).cadd(${b})`); break;
case "-": st.push(`(${a}).csub(${b})`); break;
case "*": st.push(`(${a}).cmul(${b})`); break;
case "/": st.push(`(${a}).cdiv(${b})`); break;
case "^": st.push(`(${a}).cpow(${b})`); break;
default: throw new Error("Nicht implementierter Operator: " + t.value);
}
}
continue;
}
if (t.type === "func") {
const fn = t.value.toLowerCase();
const method = FUNC_MAP[fn];
if (!method) throw new Error("Unbekannte Funktion: " + fn);
if (t.argc !== 1) throw new Error("Funktion " + fn + " erwartet genau ein Argument, erhalten: " + t.argc);
const a = st.pop();
if (a === undefined) throw new Error("Funktion " + fn + ": Stack unterlaufen.");
st.push(`(${a}).${method}()`);
continue;
}
throw new Error("Interner Compilerfehler bei Token: " + JSON.stringify(t));
}
if (st.length !== 1) throw new Error("Ausdruck ist nicht eindeutig (Rest-Stackgröße = " + st.length + ").");
return st[0];
}
// Legacy: bisherige "Kettensyntax" weiter unterstützen (optional)
function legacyExprToJS(expr) {
// entspricht dem bisherigen check_function ohne affine Einbettung und ohne Test-eval
let s = expr;
s = s.replace(/\./g, ".c")
.replace(/\[E\]/g, "Complex(Decimal(1).exp(),Decimal(0))")
.replace(/\[I\]/g, "Complex(Decimal(0),Decimal(1))")
.replace(/\[PI\]/g, "Complex(Decimal.acos(-1),Decimal(0))")
+ ".cadd([0])";
s = s.replace(/\[/g, "Complex(Decimal(").replace(/\]/g, "),Decimal(0))");
return s;
}
function isLegacySyntaxSafe(src) {
// Prevent arbitrary JavaScript in legacy mode.
// This still allows the old syntax:
// z.add([1]).log(), [1].div([1].sub(z)), [PI], [E], [I], ...
return /^[A-Za-z0-9_\.\(\)\[\]\+\-\*\/,Ee\s]+$/.test(src);
}
// Compiles the user function once to F(z)->Complex, including prefactor/shift.
// Accepted input forms:
// (A) conventional syntax: log(z+1), 1/(1-z), 2*i*exp(z), pi+e*z^3, ...
// (B) legacy syntax: z.add([1]).log(), [1].div([1].sub(z)), ...
function compileUserFunction(expr, prefac, pshift) {
const src = (expr ?? "").toString().trim();
if (src.length === 0) return null;
const looksLegacy = /[\[\]]|\b[zZ]\./.test(src);
let core;
if (looksLegacy) {
if (!isLegacySyntaxSafe(src)) {
throw new Error("Alte Syntax enthält nicht erlaubte Zeichen.");
}
const jsExpr = legacyExprToJS(src);
core = new Function(
"zz",
"Complex",
"Decimal",
"const z = zz; return (" + jsExpr + ");"
);
} else {
const jsExpr = compileMathExprToJS(src);
core = new Function(
"zz",
"Complex",
"Decimal",
"return (" + jsExpr + ");"
);
}
// Affine embedding as before:
// f(z) -> f(prefac*z - pshift)
//
// Important:
// Do not test at z=0 here. Functions such as sin(1/z), 1/z,
// log(z), z/abs(z) may be valid for the chosen sampling path.
return function (z) {
const zz = z.cmul(prefac).csub(pshift);
return core(zz, Complex, Decimal);
};
}
function decIsFinite(x) {
if (!x) return false;
if (typeof x.isFinite === "function") return x.isFinite();
if (typeof x.isNaN === "function" && x.isNaN()) return false;
var n = x.toNumber();
return Number.isFinite(n);
}
function complexIsFinite(z) {
return z &&
decIsFinite(z.re) &&
decIsFinite(z.im);
}
function safeEvalComplex(F, z) {
try {
var v = F(z);
if (!complexIsFinite(v)) return null;
return v;
} catch (e) {
return null;
}
}
function complexMagnitude1(z) {
return z.re.abs().add(z.im.abs());
}
function makeReciprocalFunction(F) {
// G(z) = 1/F(z).
// A pole of F should become a zero of G.
//
// If F(z) blows up or cannot be evaluated exactly at the pole,
// we interpret this as G(z) = 0 for the limiting root test.
return function(z) {
var v;
try {
v = F(z);
} catch (e) {
return ZER;
}
if (!complexIsFinite(v)) {
return ZER;
}
// If F(z) is nearly zero, then 1/F is singular.
// This is not a pole of F, but a zero of F.
if (complexMagnitude1(v).lt(Decimal("1e-80"))) {
return Complex(Decimal("1e80"), Decimal(0));
}
return ONE.cdiv(v);
};
}
function verifyPoleCandidate(F, z) {
// Heuristic check:
// For a genuine pole, |F(z+h)| should grow as h -> 0.
// This helps to avoid calling essential singularities or bad starts "poles".
var h1 = Complex(Decimal("1e-3"), Decimal(0));
var h2 = Complex(Decimal("1e-4"), Decimal(0));
var h3 = Complex(Decimal("1e-5"), Decimal(0));
var v1 = safeEvalComplex(F, z.cadd(h1));
var v2 = safeEvalComplex(F, z.cadd(h2));
var v3 = safeEvalComplex(F, z.cadd(h3));
if (v1 === null || v2 === null || v3 === null) {
return "Lokaler Pol-Check: nicht eindeutig, da die Funktionswerte in der Nähe nicht ausgewertet werden konnten.\n";
}
var m1 = complexMagnitude1(v1);
var m2 = complexMagnitude1(v2);
var m3 = complexMagnitude1(v3);
if (m2.gt(m1.mul(Decimal(3))) && m3.gt(m2.mul(Decimal(3)))) {
return "Lokaler Pol-Check: die Werte in der Nähe steigen rapide an; dies ist mit einem Pol vereinbar.\n";
}
return "Lokaler Pol-Check: es wurde kein eindeutiger Ausreißer festgestellt. Der Punkt könnte ein schwacher Pol,\nein schlechter Startpunkt, eine hebbare Singularität oder eine essentielle Singularität sein.\n";
}
function build_form() {
// build and write form
var form = '';
document.writeln(form);
canvas = new Raphael("canvas", twice, twice + 1);
}
function updateThumb(slider, thumbId) {
if (!slider) return;
const thumb = document.getElementById(thumbId);
if (!thumb) return;
const val = parseInt(slider.value, 10);
const min = parseInt(slider.min, 10);
const max = parseInt(slider.max, 10);
if (!isFinite(val) || !isFinite(min) || !isFinite(max) || max === min) {
thumb.textContent = slider.value;
return;
}
const percent = (val - min) / (max - min);
const sliderWidth = slider.offsetWidth || 1;
const thumbWidth = 40;
const offset = percent * (sliderWidth - thumbWidth);
thumb.style.left = Math.max(0, offset) + "px";
thumb.textContent = val;
}
// Initialisierung beim Laden, damit die Blase nicht bei 0 startet
function initSliders() {
const s1 = document.getElementById('ell');
const s2 = document.getElementById('emm');
if(s1) updateThumb(s1, 'thumb-ell');
if(s2) updateThumb(s2, 'thumb-emm');
}
// Aufruf nach dem Laden der Seite
window.addEventListener('load', initSliders);
// Und falls die Form dynamisch neu gebaut wird, auch nach build_form() aufrufen!
// Initialisierung nach dem Laden der Seite
window.addEventListener('DOMContentLoaded', (event) => {
setTimeout(() => {
updateThumb(document.getElementById('ell'), 'thumb-ell');
updateThumb(document.getElementById('emm'), 'thumb-emm');
}, 200);
});
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
// }
const d = mat.length;
for (var r = 0; r < d; r++) {
if (mat[r][r].cisZero()) {
var nzr = find_max(r);
if (nzr == d) return false;
for (var c = 0; c <= d; c++) [mat[r][c], mat[nzr][c]] = [mat[nzr][c], mat[r][c]]; // swap row
} else {
var piv = ONE.cdiv(mat[r][r]); // assert m[r][r] != 0
mat[r][r] = ONE;
for (var c = r + 1; c <= d; c++) mat[r][c] = mat[r][c].cmul(piv); // norm row
}
for (var zr = 0; zr < d; zr++) if (zr != r && !mat[zr][r].cisZero()) pivot_cell(mat[zr][r].cneg());
}
return true;
function find_max(rv) {
// find maximum in column
for (var c = rv + 1, mi = Decimal(0), ma = d; c < d; c++) if (mat[c][rv].cabs().re.gt(mi)) [mi, ma] = [mat[c][rv].cabs().re, c];
return ma;
}
function pivot_cell(za) {
// determine pivot
mat[zr][r] = ZER;
for (var c = r + 1; c <= d; c++) mat[zr][c] = mat[zr][c].cadd(za.cmul(mat[r][c]));
}
}
function check_function(realft, z, prefac, pshift) {
// Deprecated: Die Auswertung erfolgt nun über compileUserFunction / Parser.
// Diese Funktion bleibt nur aus Kompatibilitätsgründen erhalten.
try {
var F = compileUserFunction(realft, prefac, pshift);
F(z);
} catch (exception) {
return alert("Ausdruck ist nicht auswertbar!");
}
return realft;
}
function compute_cpower(result, pbsize, factor) {
result[0] = ONE;
result[1] = factor;
for (var k = 1; k <= pbsize / 2; k *= 2) for (var l = 1; l <= k; l++) result[k + l] = result[k].cmul(result[l]);
}
function compute_rpower(result, pbsize, factor) {
result[0] = Decimal(1);
result[1] = factor;
for (var k = 1; k <= pbsize / 2; k *= 2) for (var l = 1; l <= k; l++) result[k + l] = result[k].mul( result[l]);
}
function initialise_tft(signal, pbsize, F, zvalue, method, epsilo) {
var o = Decimal(0), t = new Complex((method <= 2) ? Decimal(epsilo) : Decimal(1/16384), o);
var s = [], u = new Complex(o, LZ.mul(Decimal(2).div(Decimal(pbsize)))).cexp(), z = zvalue.cadd(t);
compute_cpower(s, pbsize, u);
for (var k = 0; k < pbsize; k++, z = zvalue.cadd(t.cmul(s[k]))) {
try { signal[k] = F(z); } catch (exception) { signal[0] = null; return alert("Ausdruck ist nicht auswertbar!"); }
}
}
function initialise_fft(signal, pbsize, F) {
var e = Decimal(1e-17), o = Decimal(0), u = LZ.div(Decimal(pbsize)), x = new Complex(LZ.neg(), o), z = x;
for (var k = 0; k < 2 * pbsize;) {
try {
signal[k] = F(z);
} catch (exception) {
signal[0] = null;
return alert("Ausdruck ist nicht auswertbar!");
}
if (signal[k].re.abs().lt(e) || signal[k].re.isNaN()) signal[k].re = o;
if (signal[k].im.abs().lt(e) || signal[k].im.isNaN()) signal[k].im = o;
z = x.cadd(Complex(u.mul(Decimal(++k)), o));
}
}
function compute_fft(ins) {
const max = ins.length, one = Decimal(1), two = Decimal(2);
// reverse bits
for (var i = 0, j = 0, k, m = max >> 1; i < max - 1; i++, j += k) {
if (i < j) [ins[i], ins[j]] = [ins[j], ins[i]];
for (k = m; k <= j; k >>= 1) j -= k;
}
// compute FFT
for (var k = 0, l = 2, m = 1, c = ONE.cneg(), u = []; k < Math.log2(max); k++, m = l, l <<= 1) {
compute_cpower(u, m, c);
c = Complex(one.add(c.re).div(two).sqrt(), one.sub(c.re).div(two).sqrt().neg()); // do not negate for ifft, but normalise
for (var j = 0; j < m; j++) {
for (var i = j, n = j + m, t; i < max; i += l, n = i + m) {
t = u[j].cmul(ins[n]);
[ins[n], ins[i]] = [ins[i].csub(t), ins[i].cadd(t)];
}
}
}
}
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)) return str = "Das zugehörige lineare Gleichungssystem ist nicht lösbar: Bitte eventuell andere Vorgabe wählen.";
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(z)) ? "+" : "") + 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(z)) ? " +" : " ") + 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(z)) ? "+" : "") + 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(z)) ? " +" : " ") + denomi[k].im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n"; /* complex coefficient */
}
return str;
}
function output_taylor(result, realft, factor, intnum, ulimit, llimit, nummer, F, rootReport) {
var o = Decimal(0), v = ZER, d = ulimit.csub(llimit).cdiv(new Complex(Decimal(intnum), o)), p = (factor > 0) ? factor + 32 : 32, s;
var b = [], l = [], u = [], w = ZER, z = nummer, str = "Auf sechzehn Dezimalstellen gerundete Koeffizienten der Taylorreihe:\n";
rootReport = rootReport || "";
b[0] = llimit;
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]);
}
}
// output coefficients
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(o)) ? "+" : "") + 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(o)) ? " +" : " ") + 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), o);
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));
}
}
if (rootReport !== "") {
str = str + rootReport + "\n";
}
// output further information (e. g. for integration)
for (var j = 0; j < intnum; j++) [v, w] = [v.cadd(u[j]), w.cadd(l[j])];
str = str + "Die Summe der oberen Werte beträgt:\n";
s = ((v.re.gte(o)) ? "+" : "") + 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(o)) ? " +" : " ") + 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(o)) ? "+" : "") + 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(o)) ? " +" : " ") + 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(o)) ? "+" : "") + 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(o)) ? " +" : " ") + d.im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n";
// output zero of function and the value of the latter there
if (isNaN(nummer)) {
str = str + "Die berechnete Nullstelle liegt bei:\n";
s = ((nummer.re.gte(o)) ? "+" : "") + nummer.re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17);
s = ((nummer.im.gte(o)) ? " +" : " ") + nummer.im.toExponential(50, 4).replace("e", ")e");
str = str + s.substring(0, 17) + "(" + s.substring(17) + "im\n";
z = F(z);
str = str + "Der Funktionswert an der berechneten Nullstelle ist:\n";
s = ((z.re.gte(o)) ? "+" : "") + z.re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
str = str + s.substring(0, 17) + "(" + s.substring(17);
s = ((z.im.gte(o)) ? " +" : " ") + z.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:", z = Decimal(0);
for (var i = 0; i < knotno; i++) {
str = str + "\nKnoten " + (i + 1) + " bei ";
s = ((values[i].re.gte(z)) ? "+" : "") + 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(z)) ? "+" : "") + 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(z)) ? " +" : " ") + 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), n = 32, m = (factor > 0) ? n + factor : n, e = Decimal(epsilo), p = [], q = [], o = Decimal(1), f = o, h = new Complex(Decimal(pbsize), z);
compute_rpower(p, n, e);
for (var k = 1; k <= -factor; k++) f = f.mul(Decimal(k));
if (factor < 0) f = o.div(f);
result[0] = (factor > 0) ? ZER : result[0].cdiv(h).cmul(Complex(f, z));
for (var k = 1; k < n; k++) result[k] = result[k].cdiv(h.cmul(Complex(Decimal(p[k]), z)));
if (factor == 0) return;
// differentiate or integrate
if (factor > 0) compute_rpower(q, Math.pow(2, Math.ceil(Math.log2(factor))), e);
for (var k = n; k < m; k++) result[k] = result[k].cdiv(h.cmul(Complex(p[n] * q[k - n], z)));
for (var k = 1; k < m; k++) {
f = f.mul(Decimal(k));
if (k > factor) f = f.div(Decimal(k - factor));
result[k] = (k < factor) ? ZER : result[k].cmul(Complex(f, z));
}
}
function compute_reverse(result, pbsize, vpoint, epsilo) {
var z = Decimal(0), r = new Complex(z, LZ.mul(Decimal(2).div(Decimal(pbsize)))).cexp(), x;
var j = 1, e = Decimal(epsilo), f, h = new Complex(e, z), p = [], q = [], t = [], u = [];
if (result[1].cabs().re.lt(Decimal(1E-25))) alert("Die Umkehrung existiert nicht wegen a[1] = 0!");
// build quotients by Horner scheme
compute_rpower(p, pbsize, e);
compute_cpower(u, pbsize, ONE.cdiv(r));
for (var k = 0; k < pbsize; k++) {
x = h.cdiv(u[k]);
f = result[pbsize - 1];
for (var l = pbsize - 2; l > 0; l--) f = f.cmul(x).cadd(result[l]);
t[k] = [];
compute_cpower(t[k], pbsize, ONE.cdiv(f));
}
// revert Taylor series
result[0] = vpoint;
f = ZER;
for (var k = 0; k < pbsize; k++) f = f.cadd(t[k][1]);
result[1] = f.cdiv(Complex(Decimal(pbsize), z));
for (var l = 2; l < pbsize; l++) {
compute_cpower(q, pbsize, u[j]);
f = t[0][l];
for (var k = 1; k < pbsize; k++) f = f.cadd(t[k][l].cmul(q[k]));
result[l] = f.cdiv(Complex(Decimal(pbsize * l).mul(p[j++]), z));
}
}
function compute_fourier(result, pbsize, differ) {
var str = "Auf sechzehn Dezimalstellen gerundete Koeffizienten der Fourierreihe:\n", o = Decimal(0);
var f = new Complex(Decimal(differ), o), h = new Complex(Decimal(pbsize), o), n = pbsize / 2, s;
// differentiate if necessary
for (var j = -n; j < n; j++) {
k = (j < 0) ? j + pbsize : j;
result[k] = result[k].cdiv(h).cmul(Complex(o, 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(o)) ? "+" : "") + 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(o)) ? " +" : " ") + 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, F) {
var c0 = ZER, c1 = ONE, str = "Auf sechzehn Dezimalstellen gerundete Iterationen:\n", b = 16384, d = 1e-40;
var fa = [], ga = [], ha = [], height = [], e, f, g, h, m, o = Decimal(0), p, v = c0, w, x, y, z = zvalue;
var c2 = new Complex(Decimal(2), o), c3 = new Complex(Decimal(3), o);
var c4 = new Complex(Decimal(4), o), c5 = new Complex(Decimal(5), o);
var c6 = new Complex(Decimal(6), o), c7 = new Complex(Decimal(7), o);
var c8 = new Complex(Decimal(b), o), c9 = new Complex(o, Decimal(d));
// determine coefficients by differentiating Taylor series
height[0] = new Complex(Decimal(1 / pbsize), o);
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), o));
ga[k] = fa[k].cmul(Complex(Decimal(k - 1), o));
ha[k] = ga[k].cmul(Complex(Decimal(k - 2), o));
}
pbsize--;
// determine derivatives by using Horner's scheme
for (m = 1; m < 17; m++) {
e = z;
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 = F(z);
z = y.csub(c9);
f = F(z);
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 = F(z);
z = y.csub(c9);
f = F(z);
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 = F(z);
y = w.cmul(g);
w = x.csub(g.cmul(c2));
z = z.cadd(y.cdiv(w ));
h = F(z);
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 :
}
p = (m < 10) ? " " + m : m;
z = v.cadd(zvalue);
// output if interesting
s = ((z.re.gte(o)) ? "+" : "") + z.re.toExponential(16, 4);
str = str + p + ": Re(z) = " + s.padEnd(24) + " (" + ((z.re.gte(o)) ? "+" : "") + z.re.toExponential(80, 4) + ")\n";
s = ((z.im.gte(o)) ? "+" : "") + z.im.toExponential(16, 4);
str = str + "und Im(z) = " + s.padEnd(24) + " (" + ((z.im.gte(o)) ? "+" : "") + z.im.toExponential(80, 4) + ")\n";
if (z.re.isNaN() || z.im.isNaN() || conver < 2) m = 17;
}
str = (m == 18) ? "Die Berechnung wurde vorzeitig abgebrochen; bitte versuchen Sie es mit einer anderen Eingabe." : str + "|" + z.re.toString() + "|" + z.im.toString();
return str;
}
// ============================================================================
// START: Komplexe Nullstellen- und Polextraktion (Toleranzschlauch)
// ============================================================================
/**
* Scannt einen Streifen um die reelle Achse [-π, π] nach komplexen Nullstellen.
* Identifiziert lokale Minima des Betrags, um verborgene Singularitäten aufzuspüren.
*/
function findRootsInStrip(F, isReciprocal, steps = 1000) {
let roots = [];
let start = LZ.neg(), end = LZ;
let stepSize = end.sub(start).div(Decimal(steps));
let mags = [], xs = [];
// Vorab-Scan der reellen Achse, um lokale Minima des Betrags zu lokalisieren
for (let i = 0; i <= steps; i++) {
let x = start.add(stepSize.mul(Decimal(i)));
xs.push(x);
let val = F(Complex(x, Decimal(0)));
if (isReciprocal) val = (val.re.abs().lt(1e-100) && val.im.abs().lt(1e-100)) ? Complex(Decimal(1e100), Decimal(0)) : ONE.cdiv(val);
mags.push(val.re.abs().add(val.im.abs())); // 1-Norm
}
for (let i = 1; i < steps; i++) {
let mag = mags[i];
// Volltreffer oder lokales Minimum, das auf eine nahegelegene komplexe Nullstelle hinweist
let isDirectHit = mag.lt(Decimal("1e-8"));
let isLocalMin = (mag.lt(mags[i - 1]) && mag.lt(mags[i + 1]) && mag.lt(Decimal("0.5")));
if (isDirectHit || isLocalMin) {
// Leichte Verschiebung in die komplexe Ebene, um f'(x)=0 Sattelpunkte zu umgehen
let z0 = Complex(xs[i], Decimal("1e-5"));
let rootZ = refineComplexRoot(F, isReciprocal, z0);
// Filtere nach Nullstellen innerhalb des festgelegten Toleranzschlauchs (|Im| < 1.0)
if (rootZ !== null && rootZ.im.abs().lt(Decimal("1.0"))) {
let isNew = true;
for (let r of roots) {
if (r.csub(rootZ).re.abs().add(r.csub(rootZ).im.abs()).lt(Decimal("1e-6"))) isNew = false;
}
if (isNew) roots.push(rootZ);
// Prüfe auf die konjugiert komplexe Nullstelle
let conjZ = Complex(rootZ.re, rootZ.im.neg());
let conjRoot = refineComplexRoot(F, isReciprocal, conjZ);
if (conjRoot !== null && conjRoot.im.abs().lt(Decimal("1.0"))) {
let isNewConj = true;
for (let r of roots) {
if (r.csub(conjRoot).re.abs().add(r.csub(conjRoot).im.abs()).lt(Decimal("1e-6"))) isNewConj = false;
}
if (isNewConj) roots.push(conjRoot);
}
}
}
}
return roots;
}
/**
* Verfeinerung durch komplexes Newton-Verfahren
*/
function refineComplexRoot(F, isReciprocal, z0) {
let z = z0, h = Complex(Decimal("1e-7"), Decimal(0)), tol = Decimal("1e-12");
for (let i = 0; i < 50; i++) {
let val = F(z);
if (isReciprocal) val = (val.re.abs().lt(1e-100) && val.im.abs().lt(1e-100)) ? Complex(Decimal(1e100), Decimal(0)) : ONE.cdiv(val);
let magnitude = val.re.abs().add(val.im.abs());
if (magnitude.lt(tol)) return z;
let valH = F(z.cadd(h));
if (isReciprocal) valH = (valH.re.abs().lt(1e-100) && valH.im.abs().lt(1e-100)) ? Complex(Decimal(1e100), Decimal(0)) : ONE.cdiv(valH);
let df = valH.csub(val).cdiv(h);
// Abbruch, falls die Ableitung verschwindet, um Division durch Null zu verhindern
if (df.re.abs().add(df.im.abs()).lt(Decimal("1e-18"))) return null;
z = z.csub(val.cdiv(df));
}
return null;
}
// ============================================================================
// ENDE: Komplexe Nullstellen- und Polextraktion
// ============================================================================
/**
* Generiert die reelle, deflationierte Funktion durch direkte Produktdivision.
* Die Funktion übernimmt ein reelles x (als Decimal) und liefert den
* deflatierten Funktionswert für die Kronreihenentwicklung.
* @param {Function} F - Ursprungsfunktion
* @param {Array} rootsF - Nullstellen von f(x)
* @param {Array} rootsG - Nullstellen von 1/f(x)
* @returns {Function}
*/
function createDeflatedRealFunctionProduct(F, rootsF, rootsG) {
return function(x) { // x muss zwingend ein Decimal-Objekt sein
let tol = Decimal("1e-11");
let h = Decimal("1e-7");
// Prüfen, ob x kritisch nah an einer Singularität liegt
let nearSingularity = false;
for (let r of rootsF) if (x.sub(r).abs().lt(tol)) nearSingularity = true;
for (let p of rootsG) if (x.sub(p).abs().lt(tol)) nearSingularity = true;
// Implizite L'Hôpital-Grenzwertbildung durch minimalen Shift
let evalX = nearSingularity ? x.add(h) : x;
let val = F(Complex(evalX, Decimal(0))).re;
// Zähler-Dynamik (Nullstellen von f)
let prodF = Decimal(1);
for (let r of rootsF) prodF = prodF.mul(evalX.sub(r));
// Nenner-Dynamik (Polstellen von f / Nullstellen von g)
let prodG = Decimal(1);
for (let p of rootsG) prodG = prodG.mul(evalX.sub(p));
// Rationale Bereinigung in einem Schritt
return val.mul(prodG).div(prodF);
};
}
// ============================================================================
// ENDE: Reelle Null- und Polstellen-Abspaltung
// ============================================================================
// ============================================================================
// START: Nullstellen des Taylorpolynoms per Kron-DFT-Aberth
// ============================================================================
function complex_abs_dec(z) {
if (!z) return Decimal(0);
return z.cabs().re;
}
function complex_is_tiny(z, tol) {
return complex_abs_dec(z).lt(tol);
}
function fmt_dec16(x) {
return ((x.gte(Decimal(0))) ? "+" : "") + x.toExponential(16, 4);
}
function fmt_complex16(z) {
return fmt_dec16(z.re) + " " + fmt_dec16(z.im) + "im";
}
function next_power_two_int_gt(n) {
var k = 1;
while (k <= n) k <<= 1;
return Math.max(8, k);
}
function next_power_two_decimal(x) {
var xn = x.toNumber();
if (!Number.isFinite(xn) || xn <= 1.0) return Decimal(1);
var e = Math.ceil(Math.log2(xn));
if (e < 0) e = 0;
if (e > 4096) e = 4096;
var y = Decimal(1);
var two = Decimal(2);
for (var k = 0; k < e; k++) y = y.mul(two);
return y;
}
function extract_taylor_coeffs_for_roots(result, factor) {
// Verwendet die in output_taylor tatsächlich ausgegebenen Koeffizienten.
// Für die gewöhnliche Taylorreihe ist factor = 0.
var first = Math.max(0, factor);
var coeffs = [];
for (var k = 0; k < 32; k++) {
coeffs[k] = result[first + k] || ZER;
}
return coeffs;
}
function effective_poly_degree(coeffs) {
var maxabs = Decimal(0);
for (var k = 0; k < coeffs.length; k++) {
var a = complex_abs_dec(coeffs[k]);
if (a.gt(maxabs)) maxabs = a;
}
if (maxabs.lte(Decimal(0))) return 0;
var tol = maxabs.mul(Decimal("1e-26"));
if (tol.lt(Decimal("1e-60"))) tol = Decimal("1e-60");
for (var k = coeffs.length - 1; k >= 1; k--) {
if (complex_abs_dec(coeffs[k]).gt(tol)) return k;
}
return 0;
}
function cauchy_radius_for_local_coeffs(coeffs, N) {
if (N < 1) return Decimal(1);
var lead = complex_abs_dec(coeffs[N]);
if (lead.lte(Decimal(0))) return Decimal(1);
var m = Decimal(0);
for (var j = 0; j < N; j++) {
var a = complex_abs_dec(coeffs[j]);
if (a.gt(m)) m = a;
}
return Decimal(1).add(m.div(lead));
}
function poly_eval_coeffs(coeffs, z, N) {
var p = coeffs[N] || ZER;
for (var k = N - 1; k >= 0; k--) {
p = p.cmul(z).cadd(coeffs[k] || ZER);
}
return p;
}
function poly_eval_deriv_coeffs(coeffs, z, N) {
var p = coeffs[N] || ZER;
var dp = ZER;
for (var k = N - 1; k >= 0; k--) {
dp = dp.cmul(z).cadd(p);
p = p.cmul(z).cadd(coeffs[k] || ZER);
}
return { p: p, dp: dp };
}
function scaled_poly_residual(coeffs, z, N) {
var p = poly_eval_coeffs(coeffs, z, N);
var az = complex_abs_dec(z);
var pow = Decimal(1);
var den = Decimal(0);
for (var j = 0; j <= N; j++) {
den = den.add(complex_abs_dec(coeffs[j] || ZER).mul(pow));
pow = pow.mul(az);
}
if (den.lte(Decimal(0))) den = Decimal(1);
return complex_abs_dec(p).div(den);
}
function kron_dft_coeffs_from_taylor(coeffs, N, rho, K) {
// Lokales Taylorpolynom:
// T(h) = a_0 + ... + a_N h^N
//
// Kronform:
// h = rho*q, |q| = 1
// P(q) = T(rho*q) = d_0 + ... + d_N q^N
//
// Die d_j werden bewusst per direkter DFT aus den Kronwerten rekonstruiert.
// Damit entspricht die Erweiterung dem Kron-DFT-Aberth-Verfahren.
var samples = [];
var d = [];
var upow = [];
var zero = Decimal(0);
var rhoC = Complex(rho, zero);
var u = Complex(
zero,
LZ.mul(Decimal(2).div(Decimal(K)))
).cexp();
compute_cpower(upow, K, u);
for (var k = 0; k < K; k++) {
var h = rhoC.cmul(upow[k]);
samples[k] = poly_eval_coeffs(coeffs, h, N);
}
var invK = ONE.cdiv(Complex(Decimal(K), zero));
for (var j = 0; j <= N; j++) {
var s = ZER;
for (var k = 0; k < K; k++) {
var idx = (K - ((j * k) % K)) % K;
s = s.cadd(samples[k].cmul(upow[idx]));
}
d[j] = s.cmul(invK);
}
return d;
}
function aberth_roots_coeffs(coeffs, N) {
var roots = [];
var zero = Decimal(0);
var oneD = Decimal(1);
var tiny = Decimal("1e-60");
var RA = cauchy_radius_for_local_coeffs(coeffs, N);
if (RA.lt(oneD)) RA = oneD;
// Halber Winkelschritt gegen einfache Symmetrie-Sonderfälle.
for (var i = 0; i < N; i++) {
var ang = 2.0 * Math.PI * (i + 0.5) / N;
roots[i] = Complex(
RA.mul(Decimal(Math.cos(ang))),
RA.mul(Decimal(Math.sin(ang)))
);
}
var tolStep = Decimal("1e-20");
var tolRes = Decimal("1e-20");
var maxIter = 100;
var minIter = 5;
var maxStep = Decimal(0);
var maxRes = Decimal(0);
for (var it = 0; it < maxIter; it++) {
var next = [];
maxStep = Decimal(0);
maxRes = Decimal(0);
for (var i = 0; i < N; i++) {
var ev = poly_eval_deriv_coeffs(coeffs, roots[i], N);
var p = ev.p;
var dp = ev.dp;
var res = scaled_poly_residual(coeffs, roots[i], N);
if (res.gt(maxRes)) maxRes = res;
if (complex_abs_dec(dp).lt(tiny)) {
// Sehr flache Stelle: minimal stören statt Division durch fast Null.
next[i] = roots[i].cadd(Complex(tiny, tiny));
continue;
}
var delta = p.cdiv(dp);
var rep = ZER;
for (var j = 0; j < N; j++) {
if (j === i) continue;
var diff = roots[i].csub(roots[j]);
if (complex_abs_dec(diff).lt(tiny)) {
diff = diff.cadd(Complex(tiny, tiny));
}
rep = rep.cadd(ONE.cdiv(diff));
}
var denom = ONE.csub(delta.cmul(rep));
if (complex_abs_dec(denom).lt(tiny)) {
denom = denom.cadd(Complex(tiny, tiny));
}
var step = delta.cdiv(denom);
next[i] = roots[i].csub(step);
var st = complex_abs_dec(step);
if (st.gt(maxStep)) maxStep = st;
}
roots = next;
var maxRoot = Decimal(1);
for (var i = 0; i < N; i++) {
var ar = complex_abs_dec(roots[i]);
if (ar.gt(maxRoot)) maxRoot = ar;
}
if (it >= minIter &&
maxStep.lte(tolStep.mul(maxRoot)) &&
maxRes.lte(tolRes)) {
break;
}
}
roots.sort(function(a, b) {
var dr = a.re.sub(b.re).toNumber();
if (Math.abs(dr) > 1e-14) return dr;
return a.im.sub(b.im).toNumber();
});
return {
roots: roots,
maxStep: maxStep,
maxRes: maxRes
};
}
function output_taylor_roots_from_coeffs(result, factor, centre, F) {
var coeffs = extract_taylor_coeffs_for_roots(result, factor);
var N = effective_poly_degree(coeffs);
var str = "";
var rootsZ = [];
str += "\nNullstellen des Taylorpolynoms\n";
str += "--------------------------------\n";
if (N < 1) {
str += "Das Taylorpolynom besitzt nach der aktuellen Trunkation keinen positiven Grad.\n";
return { text: str, rootsZ: rootsZ };
}
var RC = cauchy_radius_for_local_coeffs(coeffs, N);
var rho = next_power_two_decimal(RC);
var K = next_power_two_int_gt(2 * N);
var d = kron_dft_coeffs_from_taylor(coeffs, N, rho, K);
var sol = aberth_roots_coeffs(d, N);
var rhoC = Complex(rho, Decimal(0));
str += "Berechnet werden die Nullstellen des abgeschnittenen Taylorpolynoms\n";
str += "T(h)=a_0+a_1*h+...+a_N*h^N mit h=z-z_0.\n\n";
str += "Entwicklungspunkt z_0 = " + fmt_complex16(centre) + "\n";
str += "effektiver Grad N = " + N + "\n";
str += "Cauchy-Radius R_C = " + RC.toExponential(16, 4) + "\n";
str += "Kronradius rho = " + rho.toExponential(16, 4) + "\n";
str += "Kronpunkte K = " + K + "\n";
str += "max. skaliertes Res. = " + sol.maxRes.toExponential(16, 4) + "\n\n";
str += "Gefundene Nullstellen z_i = z_0 + rho*q_i:\n";
var maxT = Decimal(0);
var maxF = Decimal(0);
for (var i = 0; i < sol.roots.length; i++) {
var q = sol.roots[i];
var h = rhoC.cmul(q);
var z = centre.cadd(h);
rootsZ.push(z);
var tv = poly_eval_coeffs(coeffs, h, N);
var at = complex_abs_dec(tv);
if (at.gt(maxT)) maxT = at;
var fv = safeEvalComplex(F, z);
if (fv !== null) {
var af = complex_abs_dec(fv);
if (af.gt(maxF)) maxF = af;
}
str += "z_" + (i + 1) + " = " + fmt_complex16(z) + "\n";
}
str += "\nKontrolle am Taylorpolynom:\n";
str += "max. |T(z_i-z_0)| = " + maxT.toExponential(16, 4) + "\n";
str += "\nHinweis:\n";
str += "Dies sind Nullstellen des abgeschnittenen Taylorpolynoms. Sie sind nur dann\n";
str += "auch gute Nullstellen der Originalfunktion, wenn das Taylorpolynom im betreffenden\n";
str += "Gebiet eine hinreichend gute Näherung der Funktion ist.\n";
if (maxF.gt(Decimal(0))) {
str += "max. |F(z_i)| = " + maxF.toExponential(16, 4) + "\n";
}
return {
text: str,
rootsZ: rootsZ
};
}
// ============================================================================
// ENDE: Nullstellen des Taylorpolynoms per Kron-DFT-Aberth
// ============================================================================
function compute_coeffs() {
var deg, dif, eps, fac, fct, inp = [], kno = [], pee = [], pol = [], quu = [], o = Decimal(0), x = o, index = document.fft.b.selectedIndex;
var ell, emm = 0, inn, lli, lev, llr, lol, num, pfi, pfr, poi, por, shi, shr, str = "", uli, ulr, upl, vpf, vpo, vsh, f = 1, z;
var fun = [
"log(z+1)",
"atan(z)",
"1/(1-z)",
"2*i*exp(z)",
"pi + e*z^3",
"cos(2*z)",
"abs(z)",
"z/abs(z)",
"abs(sin(z))",
"sin(1/z)"
];
// check function and initialise variables
if (index > 0) document.fft.pbsize.value = fun[index - 1];
fct = document.fft.pbsize.value.toString();
// Avoid alert on first load / empty input
if (!fct || !fct.trim()) {
return true;
}
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));
var F;
try {
F = compileUserFunction(fct, vpf, vsh);
} catch (e) {
return alert(e.message);
}
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 != 3) dif = -Number(document.fft.ite.value); else document.fft.ite.value = 0;
if (Number(document.fft.exp.value) > 0) exp = Math.pow(2, Number(document.fft.exp.value));
switch (opt) {
// compute (reverse) Taylor series
case 1 :
case 2 :
initialise_tft(inp, exp * 2, F, vpo, opt, eps);
if (inp[0] !== null) {
compute_fft(inp);
compute_taylor(inp, exp * 2, dif, eps);
if (opt == 2) {
z = inp[0].cneg();
compute_reverse(inp, exp * 2, vpo, eps);
// Nullstelle der Umkehrung
f = inp[exp * 2 - 1];
for (var l = exp * 2 - 2; l >= 0; l--) {
f = f.cmul(z).cadd(inp[l]);
}
vpo = z.cneg();
}
var taylorRootInfo = null;
var taylorRootReport = "";
if (opt == 1 && document.fft.troot && document.fft.troot.checked) {
taylorRootInfo = output_taylor_roots_from_coeffs(inp, dif, vpo, F);
taylorRootReport = taylorRootInfo.text;
}
document.fft.result.value = output_taylor(inp, fct, dif, inn, upl, lol, f, F, taylorRootReport);
quu[0] = f;
// Taylor-Nullstellen für die Grafik anhängen
if (taylorRootInfo && taylorRootInfo.rootsZ) {
quu.taylorRoots = taylorRootInfo.rootsZ;
} else {
delete quu.taylorRoots;
}
pee = inp;
ell = exp;
}
break;
// compute Fourier series
case 3 :
exp /= 2;
initialise_fft(inp, exp, F);
if (inp[0] !== null) {
compute_fft(inp);
document.fft.result.value = compute_fourier(inp, exp * 2, 0);
quu[0] = f;
pee = inp;
ell = exp;
}
exp /= 2;
break;
// compute zero of function
case 4 :
exp /= 2;
lev = Number(document.fft.lev.value);
// Für reine Iterationen (lev <= 1) genügt eine Auswertung; sonst Taylor-Daten via TFT/FFT
initialise_tft(inp, exp, F, vpo, opt, eps);
if (inp[0] !== null) {
if (lev > 1) compute_fft(inp);
quu[0] = compute_simpson(inp, exp, vpo, lev, F);
pee = quu[0].split("|");
document.fft.result.value = pee[0];
quu[0] = Complex(
Decimal(pee[1]),
Decimal((pee.length > 2) ? pee[2] : "0")
);
pee = inp;
ell = exp;
}
exp *= 2;
break;
// Compute pole of function via root search applied to 1/F
case 5 :
exp /= 2;
lev = Number(document.fft.lev.value);
var G = makeReciprocalFunction(F);
// Same machinery as for roots, but applied to G = 1/F.
initialise_tft(inp, exp, G, vpo, 4, eps);
if (inp[0] !== null) {
if (lev > 1) compute_fft(inp);
quu[0] = compute_simpson(inp, exp, vpo, lev, G);
pee = quu[0].split("|");
var poleCandidate = Complex(
Decimal(pee[1]),
Decimal((pee.length > 2) ? pee[2] : "0")
);
var out = "";
out += "Polsuche über die Nullstellenbestimmung von 1/F.\n";
out += "-------------------------------------------------\n\n";
out += pee[0] + "\n";
out += "Berechneter Polkandidat:\n";
out += "Re(z) = " + poleCandidate.re.toString() + "\n";
out += "Im(z) = " + poleCandidate.im.toString() + "\n\n";
out += verifyPoleCandidate(F, poleCandidate);
document.fft.result.value = out;
quu[0] = poleCandidate;
pee = inp;
ell = exp;
}
exp *= 2;
break;
// Berechne Padé-Approximation (Erweiterter Schutz gegen Singularitäten/NaN)
case 6 :
// Ermittle native komplexe Nullstellen mittels des neuen Toleranzschlauch-Scans
var rootsF = findRootsInStrip(F, false, 2000);
var rootsG = findRootsInStrip(F, true, 2000);
var Nf = rootsF.length, Ng = rootsG.length;
// Komplexe Deflations-Logik
var F_deflated = function(z) {
var tol = Decimal("1e-11"), ch = new Complex(Decimal("1e-7"), Decimal(0)), nearSingularity = false;
for (var i = 0; i < rootsF.length; i++) if (z.csub(rootsF[i]).re.abs().add(z.csub(rootsF[i]).im.abs()).lt(tol)) nearSingularity = true;
for (var i = 0; i < rootsG.length; i++) if (z.csub(rootsG[i]).re.abs().add(z.csub(rootsG[i]).im.abs()).lt(tol)) nearSingularity = true;
var evalZ = nearSingularity ? z.cadd(ch) : z;
var val = F(evalZ), prodF = ONE, prodG = ONE;
// Direkte Anwendung der Complex-Objekte
for (var i = 0; i < rootsF.length; i++) prodF = prodF.cmul(evalZ.csub(rootsF[i]));
for (var i = 0; i < rootsG.length; i++) prodG = prodG.cmul(evalZ.csub(rootsG[i]));
return val.cmul(prodG).cdiv(prodF);
};
initialise_tft(inp, exp * 2, F_deflated, vpo, 1, eps);
if (inp[0] !== null) {
compute_fft(inp);
compute_taylor(inp, exp * 2, 0, eps);
var ui_ell = Number(document.fft.ell.value), ui_emm = Number(document.fft.emm.value);
// Strukturelle Analyse
var rest_ell_opt = 0, noise = Decimal("1e-10");
for(var k = 0; k < 30; k++) if (inp[k].re.abs().gt(noise) || inp[k].im.abs().gt(noise)) rest_ell_opt = k;
var sug_L = (rest_ell_opt < 25) ? rest_ell_opt : 4;
var sug_M = (rest_ell_opt < 25) ? 0 : 4;
var opt_L = sug_L + Nf, opt_M = sug_M + Ng;
var outStr = "VORGEZOGENE ABSPALTUNG VON LINEARFAKTOREN\n-----------------------------------------\n\n";
outStr += "1. ERKANNTE STRUKTUR (heuristische Analyse)\n";
outStr += "Vorgeschlagene Mindestgrade: Zählergrad " + opt_L +
" und Nennergrad " + opt_M + "\n\n";
var rest_L = ui_ell - Nf, rest_M = ui_emm - Ng;
if (rest_L < 0 || rest_M < 0) {
outStr += "2. FEHLER: Die ausgewählten Grade sind für die gefundenen Nullstellen/Pole zu klein.\n";
document.fft.result.value = outStr;
} else if (rest_L + rest_M > 31) {
outStr += "2. FEHLER: Die ausgewählten Padé-Grade erfordern zu viele Taylorkoeffizienten.\n";
outStr += "Für die aktuelle Taylor-Ordnung müssen die abgespaltenen Grade folgende Bedingung erfüllen: ";
outStr += "Restzählergrad + Restnennergrad <= 31.\n";
outStr += "Aktuell: " + rest_L + " + " + rest_M + " = " + (rest_L + rest_M) + ".\n";
document.fft.result.value = outStr;
} else {
var sel_p = [], sel_q = [];
for(var i=0; i<64; i++) { sel_p[i] = ZER; sel_q[i] = ZER; }
var pade_msg = output_pade(inp, 0, rest_L, rest_M, sel_p, sel_q);
var hasNaN = false;
for(var i=0; i<=rest_L; i++) if(sel_p[i].re.isNaN()) hasNaN = true;
for(var i=0; i<=rest_M; i++) if(sel_q[i].re.isNaN()) hasNaN = true;
var isSolvable = (pade_msg.indexOf("nicht lösbar") === -1) && !hasNaN;
if (!isSolvable) {
for(var i=0; i<64; i++) {
sel_p[i] = (i <= rest_L) ? (inp[i] || ZER) : ZER;
sel_q[i] = (i === 0) ? ONE : ZER;
}
}
var reconstruct = function(p_coeffs, q_coeffs, rF, rG, centre) {
var resP = [];
for (var i = 0; i < 64; i++) {
resP[i] = (p_coeffs[i]) ? p_coeffs[i] : ZER;
}
// Multipliziere Zähler mit den lokalen Faktoren (w - r_loc),
// wobei w = z - centre und r_loc = r - centre sind.
for (var i = 0; i < rF.length; i++) {
var rloc = rF[i].csub(centre);
var temp = [];
for (var k = 0; k < 64; k++) {
temp[k] = ((k > 0) ? resP[k - 1] : ZER).csub(resP[k].cmul(rloc));
}
resP = temp;
}
var resQ = [];
for (var i = 0; i < 64; i++) {
resQ[i] = (q_coeffs[i]) ? q_coeffs[i] : ZER;
}
// Multipliziere Nenner mit den lokalen Faktoren (w - p_loc).
for (var i = 0; i < rG.length; i++) {
var ploc = rG[i].csub(centre);
var temp = [];
for (var k = 0; k < 64; k++) {
temp[k] = ((k > 0) ? resQ[k - 1] : ZER).csub(resQ[k].cmul(ploc));
}
resQ = temp;
}
return { p: resP, q: resQ };
};
var finalSel = reconstruct(sel_p, sel_q, rootsF, rootsG, vpo);
outStr += "2. KOEFFIZIENTEN DER GEWÄHLTEN APPROXIMATION\n";
if (!isSolvable) outStr += "Hinweis: System ist singulär (Überanpassung). Verwende Taylor-Fallback.\n";
outStr += "Auf sechzehn Dezimalstellen gerundete Koeffizienten (Grad " + ui_ell + "/" + ui_emm + "):\n";
var formatHighPrec = function(val, label, k) {
if (!val || val.re.isNaN()) val = ZER;
var s, res = (k < 10) ? label + "[0" + k + "] = " : label + "[" + k + "] = ";
s = ((val.re.gte(o)) ? "+" : "") + val.re.toExponential(50, 4).replace("e", ")e");
s = s.padEnd(59);
res += s.substring(0, 17) + "(" + s.substring(17);
s = ((val.im.gte(o)) ? " +" : " ") + val.im.toExponential(50, 4).replace("e", ")e");
res += s.substring(0, 17) + "(" + s.substring(17) + "im\n";
return res;
};
for (var k = 0; k <= ui_ell; k++) outStr += formatHighPrec(finalSel.p[k], "p", k);
for (var k = 0; k <= ui_emm; k++) outStr += formatHighPrec(finalSel.q[k], "q", k);
document.fft.result.value = outStr;
ell = ui_ell; emm = ui_emm;
for(var i = 0; i < 64; i++) { pee[i] = finalSel.p[i] || ZER; quu[i] = finalSel.q[i] || ZER; }
}
}
break;
// compute spline approximation
case 7 :
initialise_tft(inp, exp * 2, F, vpo, 1, eps);
if (inp[0] !== null) {
deg = Number(document.fft.deg.value);
num = Number(document.fft.num.value);
compute_fft(inp);
compute_taylor(inp, exp * 2, 0, eps);
kno[0] = new Complex(LZ.neg(), o);
lev = new Complex(LZ.mul(Decimal(2).div(Decimal(num - 1))), o);
for (var i = 0; i < num; i++) {
pol[i] = [];
initialise_tft(pol[i], exp * 2, F, kno[i], 1, eps);
compute_fft(pol[i]);
compute_taylor(pol[i], exp * 2, 0, eps);
kno[i+1] = kno[i].cadd(lev);
}
document.fft.result.value = output_taylor(inp, fct, 0, inn, upl, lol, 1, F) + output_spline(pol, kno, deg, num);
quu[0] = f;
pee = pol;
ell = num;
emm = deg;
}
break;
default :
}
if (inp[0] === null) document.fft.result.value = str + "Auswertung ist nicht möglich.";
if (fct > "") draw_graph(exp, F, com, opt, ell, emm, pee, quu, vpo, dif, kno, pol);
return true;
}
function draw_graph(asize, F, compl, optio, eldeg, emdeg, numer, denom, value, diffe, knots, polys) {
function drawRootCrosses(rootList, colour) {
if (!rootList || !rootList.length) return;
for (var i = 0; i < rootList.length; i++) {
var mark = rootList[i];
if (!mark) continue;
// Wie bei der einfachen Nullstellenbestimmung:
// Nur nahe reelle Nullstellen auf dieser reellen Grafik markieren.
if (mark.im.abs().gt(Decimal("1e-8"))) continue;
var xpix = CENTRE + Math.round(CENTRE * mark.re.toNumber() / Math.PI);
if (xpix < 0 || xpix > twice) continue;
canvas.text(xpix, CENTRE, "X").attr({fill: colour});
}
}
var bnd = 1e5, nul = Decimal(0), fac = new Complex(Decimal(-diffe), nul), fra = new Complex(LZ.div(Decimal(asize)), nul), graph = [], inp, padeg = [];
var max = -bnd, min = bnd, inc, len, meths = [" ", "Taylor", "Umkehrung", "Fourier", "Nullstelle", "Pol ", "Padé ", "Spline"], k = -twice, m = 2 * asize;
var den, hlp, tst, r = asize - 0.5, quo = 2 * r / twice, z = Complex(LZ, nul).cneg(), end_x, end_y, horizontal_axis, vertical_axis, start_x, start_y;
switch (optio) {
// compute (reverse) Taylor series
case 1:
case 2:
z = z.csub(value);
inc = Math.max(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
case 3:
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(nul, Decimal(-l))).cexp().cmul(numer[m - l]));
padeg[n] = hlp;
z = z.cadd(fra);
}
break;
// compute Padé approximant
case 6:
var z_pade = z.csub(value); // Lokale Kopie von z für den Padé-Scan
for (var n = 0; n < 2 * asize; n++) {
var hlp = numer[eldeg] || ZER;
for (var l = eldeg - 1; l >= 0; l--) hlp = hlp.cmul(z_pade).cadd(numer[l] || ZER);
var den = denom[emdeg] || ONE;
// Variable 'm_loop' statt 'm' verwenden, um Kollision zu vermeiden
for (var m_loop = emdeg - 1; m_loop >= 0; m_loop--) den = den.cmul(z_pade).cadd(denom[m_loop] || ZER);
if (den.cabs().re.lt(Decimal("1e-25"))) den = ONE; // Schutz gegen Division durch 0
padeg[n] = hlp.cdiv(den);
z_pade = z_pade.cadd(fra);
}
break;
// compute spline approximation
case 7:
fac = new Complex(LZ.mul(Decimal(2).div(Decimal(eldeg - 1))), nul);
eldeg = 0;
for (var n = 0; n < 2 * asize; n++) {
hlp = polys[eldeg][emdeg];
z = z.csub(knots[eldeg]);
// FIX: m_loop statt m nutzen
for (var m_loop = emdeg - 1; m_loop >= 0; m_loop--) hlp = hlp.cmul(z).cadd(polys[eldeg][m_loop]);
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_loop = emdeg - 1; m_loop >= 0; m_loop--) hlp = hlp.cmul(z).cadd(polys[eldeg][m_loop]);
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;
default :
}
// determine minimum and maximum
z = Complex(LZ.neg(), nul);
for (var n = 0; n < 2 * asize; n++) {
graph[n] = F(z);
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);
if (bnd < 1 || !isFinite(bnd)) bnd = 1;
len = bnd.toString().length;
// draw axes
canvas.clear();
start_x = start_y = 0;
end_x = 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});
var legendPrefix = (optio == 4 || optio == 5) ? "X " : "--- ";
canvas.text(
meths[optio].length * 3 + 5,
20,
legendPrefix + meths[optio]
).attr({fill: "#FF0000"});
// draw curve of 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 function root in green, then second result in red if available
if ((optio == 4 || optio == 5) && denom && denom[0]) {
var mark = denom[0];
// Only mark it on this real-axis graph if the imaginary part is tiny.
if (mark.im.abs().lt(Decimal("1e-8"))) {
var markColour = (optio == 5) ? "#FFFF00" : "#00FF00";
canvas.text(
CENTRE + Math.round(CENTRE * mark.re.toNumber() / Math.PI),
CENTRE,
"X"
).attr({fill: markColour});
}
}
// optio == 0 means: draw only the original function.
// optio == 4 means: root-finding.
// optio == 5 means: pole-finding via 1/F.
// In both cases no red approximation curve is required.
if (optio == 0 || optio == 4 || optio == 5) return;
// RESET FÜR ROTE KURVE
start_x = 0;
k = -twice;
inc = quo * k / 2 + r;
while (++k < twice) {
if (k % 2 == 0) {
var val_pad = padeg[Math.round(inc)];
if (val_pad && !val_pad.re.isNaN()) {
start_y = (compl == 1) ? CENTRE - Math.round(bnd * val_pad.re.toNumber()) : CENTRE - Math.round(bnd * val_pad.im.toNumber());
inc += quo;
var next_pad = padeg[Math.round(inc)] || val_pad;
end_y = (compl == 1) ? CENTRE - Math.round(bnd * next_pad.re.toNumber()) : CENTRE - Math.round(bnd * next_pad.im.toNumber());
var red_path = canvas.path("M" + (start_x - 1) + " " + start_y + "L" + ++start_x + " " + end_y);
red_path.attr({"stroke": "#FF0000", "stroke-width": 1.5 }); // Dickere rote Linie zur besseren Sichtbarkeit
} else {
inc += quo;
start_x++;
}
}
}
if (optio == 1 && denom && denom.taylorRoots) {
drawRootCrosses(denom.taylorRoots, "#00FF00");
}
}
// main programme
build_form();
document.fft.result.value =
"Kurze Bedienungsanleitung\n" +
"=========================\n\n" +
"1. Funktion auswählen oder eingeben\n" +
"Wählen Sie links oben unter „Beispiele“ eine Beispielfunktion aus oder geben Sie\n" +
"eine eigene Formel in der Variablen z ein, zum Beispiel:\n" +
" log(z+1)\n" +
" atan(z)\n" +
" 1/(1-z)\n" +
" 2*i*exp(z)\n" +
" pi + e*z^3\n" +
" sin(1/z)\n\n" +
"Erlaubte Konstanten:\n" +
" pi, e, i\n\n" +
"Erlaubte Rechenzeichen:\n" +
" +, -, *, /, ^ oder **\n" +
"Auch implizite Multiplikation ist möglich, zum Beispiel 2z oder 2(z+1).\n\n" +
"Erlaubte Funktionen:\n" +
" exp, log, ln, sqrt, abs,\n" +
" sin, cos, tan, cot, sec, csc,\n" +
" asin, acos, atan, acot, asec, acsc,\n" +
" sinh, cosh, tanh, coth, sech, csch,\n" +
" asinh, acosh, atanh, acoth, asech, acsch.\n\n" +
"2. Realteil oder Imaginärteil\n" +
"Mit „Re f“ wird der Realteil der Funktion gezeichnet und ausgewertet.\n" +
"Mit „Im f“ wird der Imaginärteil der Funktion gezeichnet und ausgewertet.\n\n" +
"3. Verfahren\n" +
" Taylorreihe berechnet Taylor-Koeffizienten am Punkt; mit Checkbox auch Nullstellen des Taylorpolynoms.\n" +
" Umkehrung berechnet die Umkehrreihe einer Taylorreihe.\n" +
" Fourierreihe berechnet Fourier-Koeffizienten auf [-pi, pi].\n" +
" Nullstelle sucht eine Nullstelle der Funktion.\n" +
" Pol sucht einen Pol über die Nullstellensuche von 1/F.\n" +
" Padé-Approximation berechnet eine rationale Padé-Approximation.\n" +
" Spline-Approximation berechnet lokale Taylor-/Spline-Approximationen.\n\n" +
"4. Parameter\n" +
" Grad Grad der lokalen Spline-Polynome.\n" +
" Stufe Iterations- bzw. Konvergenzstufe bei Nullstelle und Pol.\n" +
" -lb(Epsilon) bestimmt Epsilon = 2^(-Wert).\n" +
" Ableitung Ordnung der Ableitung.\n" +
" Integral Ordnung der Integration, falls keine Ableitung gewählt ist.\n" +
" Exponent steuert die Auflösung der Grafik.\n" +
" Knoten Anzahl der Knoten für die Spline-Approximation.\n" +
" Anzahl Anzahl der Teilintervalle für Summen/Integrale.\n" +
" Zählergrad Grad des Padé-Zählers.\n" +
" Nennergrad Grad des Padé-Nenners.\n\n" +
"5. Komplexe Eingabewerte\n" +
" Punkt Entwicklungs- bzw. Startpunkt.\n" +
" Obergrenze obere Grenze für Auswertungen, Summen oder Integrale.\n" +
" Untergrenze untere Grenze für Auswertungen, Summen oder Integrale.\n" +
" Vorfaktor Vorfaktor der affinen Einbettung.\n" +
" Verschiebung Verschiebung der affinen Einbettung.\n\n" +
"Alle diese Werte besitzen jeweils einen Realteil und einen Imaginärteil.\n" +
"Intern wird die Funktion als f(Vorfaktor*z - Verschiebung) ausgewertet.\n\n" +
"6. Grafik\n" +
"Die weiße Kurve zeigt die eingegebene Funktion.\n" +
"Die rote Kurve zeigt die berechnete Approximation, also Taylorreihe,\n" +
"Umkehrung, Fourierreihe, Padé-Approximation oder Spline-Approximation.\n" +
"Bei Nullstelle und Pol wird keine rote Approximation gezeichnet; stattdessen wird der\n" +
"gefundene Punkt grün bzw. gelb markiert, sofern er nahe genug zur reellen Achse liegt.\n\n" +
"7. Hinweise\n" +
"Bei Singularitäten, Verzweigungspunkten oder sehr nahe an Polen können\n" +
"Taylorreihe und Padé-Approximation empfindlich reagieren.\n" +
"Für Funktionen wie sin(1/z) ist oft ein Punkt abseits von z = 0 sinnvoll,\n" +
"zum Beispiel Punkt: Realteil = 1, Imaginärteil = 0.\n" +
"Die Spline-Approximation ist bei stark oszillierenden Funktionen oft robuster\n" +
"als eine einzelne globale Padé-Approximation.\n\n" +
"Nach Eingabe oder Änderung eines Parameters bitte „Berechnen“ drücken.";
compute_coeffs();