// [url=https://www.thomann.de/de/moog_synthesizer.html]Moog[/url]-Ladder als Schaltungsmodell: implizite Trapezregel, je Sample per Newton gelöst.
// Ein Template für zwei Lanetypen:
// V = double -> eine Stimme
// V = Vec -> vier Stimmen je AVX2-Register
// Beide Instanzen führen je Lane dieselben IEEE-Operationen in derselben Reihenfolge aus
// (kein FMA) und liefern damit bitgleiche Ergebnisse.
#include <algorithm>
#include <array>
#include <cmath>
// tanh als rationale Minimax-Näherung x*P(x^2)/Q(x^2), Grad 13/6 (Koeffizienten aus Eigens
// generic_fast_tanh_float, hier in double ausgewertet), |Fehler| <= 2.2e-7, ungerade, monoton,
// |t| < 1. Estrin statt Horner: kürzere Abhängigkeitskette, denn die Newton-Schleife ist
// latenz-, nicht durchsatzgebunden.
inline constexpr double fastTanhLimit = 7.99881172180175781;
template <class V>
inline V fastTanhClamped(V x) noexcept
{
constexpr double a1 = 4.89352455891786e-03, a3 = 6.37261928875436e-04,
a5 = 1.48572235717979e-05, a7 = 5.12229709037114e-08,
a9 = -8.60467152213735e-11, a11 = 2.00018790482477e-13,
a13 = -2.76076847742355e-16;
constexpr double b0 = 4.89352518554385e-03, b2 = 2.26843463243900e-03,
b4 = 1.18534705686654e-04, b6 = 1.19825839466702e-06;
const V x2 = x * x, x4 = x2 * x2, x8 = x4 * x4;
const V p = (a1 + a3 * x2) + x4 * (a5 + a7 * x2) + x8 * ((a9 + a11 * x2) + x4 * a13);
const V q = (b0 + b2 * x2) + x4 * (b4 + b6 * x2);
return x * p / q;
}
// Skalare Lane-Primitive; der AVX2-Typ stellt dieselben Funktionen bereit.
inline double laneAbs(double x) noexcept { return std::abs(x); }
inline double laneClamp(double x, double lo, double hi) noexcept { return std::clamp(x, lo, hi); }
inline double laneSelect(bool m, double a, double b) noexcept { return m ? a : b; }
inline bool laneAny(bool m) noexcept { return m; }
inline double laneMax(double change, double x) noexcept { return std::max(change, x); }
inline bool laneLess(double a, double b) noexcept { return a < b; }
inline bool laneAndNot(bool a, bool b) noexcept { return a && !b; }
template <class V>
V laneTanh(V x) noexcept
{
return fastTanhClamped(laneClamp(x, -fastTanhLimit, fastTanhLimit));
}
// Startwert aus der linearen Extrapolation der beiden vorigen Schritte:
// Fehler zweiter statt erster Ordnung in der Schrittweite.
template <class V>
void start(std::array<V, 4>& y, std::array<V, 4>& previous, int count) noexcept
{
for (int i = 0; i < count; ++i) {
const V guess = 2 * y[i] - previous[i];
previous[i] = y[i];
y[i] = guess;
}
}
// Abbruch bei |dy|^2 < 0.01 * tol: Bei quadratischer Konvergenz liegt der Restfehler danach
// bei C*|dy|^2, weit unter der Schwelle. Der reine Bestätigungsschritt entfällt
// (im Mittel 2 statt 3 Iterationen).
template <class V, class M>
M converged(M active, V change, double tolerance) noexcept
{
return laneAndNot(active, laneLess(change * change, .01 * tolerance));
}
// Zustandsupdate der Trapezregel; sehr kleine Werte gehen auf null, bevor sie denormal werden.
template <class V>
void integrate(std::array<V, 4>& s, const std::array<V, 4>& y, int count) noexcept
{
for (int i = 0; i < count; ++i) {
s[i] = 2 * y[i] - s[i];
s[i] = laneSelect(laneLess(laneAbs(s[i]), 1e-30), V(0.0), s[i]);
}
}
// Ladder: dy_i/dt = w * (u_i - tanh(y_i)), u_0 = tanh(x - k*y_3), u_i = tanh(y_(i-1)).
// Trapezregel => y = s + g * (u(y) - tanh(y)).
//
// Die Jacobi-Matrix ist bidiagonal plus ein Rückkopplungseintrag (0,3). Statt allgemeiner LU:
// Vorwärtssubstitution mit dy_3 als Unbekannter (dy_i = p_i + q_i * dy_3), danach
// dy_3 = p_3 / (1 - q_3). Exakt, O(n), ohne Pivotsuche. Die Kehrwerte der Diagonale hängen
// nur von y ab und liegen nicht auf der seriellen Kette p_(i-1) -> p_i.
//
// Jede Lane hat eine eigene Konvergenzmaske: konvergierte Lanes frieren ein. Das Ergebnis einer
// Stimme hängt damit nicht davon ab, mit welchen Stimmen sie ein Register teilt.
template <class V, class M>
V ladder(std::array<V, 4>& s, std::array<V, 4>& y, std::array<V, 4>& previous,
V x, V g, V k, double tolerance, M active) noexcept
{
start(y, previous, 4);
for (int iteration = 0; iteration < 8 && laneAny(active); ++iteration) {
std::array<V, 4> t{}, slope{}, p{}, q{};
for (int i = 0; i < 4; ++i) {
t[i] = laneTanh(y[i]);
slope[i] = 1 - t[i] * t[i]; // tanh' aus dem bereits berechneten tanh
}
const V feedback = laneTanh(x - k * y[3]);
for (int i = 0; i < 4; ++i) {
const V drive = i == 0 ? feedback : t[i - 1];
const V inverse = 1 / (1 + g * slope[i]);
const V residual = (y[i] - s[i] - g * (drive - t[i])) * inverse;
if (i == 0) {
p[0] = -residual;
q[0] = -g * k * (1 - feedback * feedback) * inverse;
} else {
const V lower = g * slope[i - 1] * inverse;
p[i] = lower * p[i - 1] - residual;
q[i] = lower * q[i - 1];
}
}
const V delta3 = p[3] / (1 - q[3]);
V change(0.0);
for (int i = 0; i < 4; ++i) {
const V delta = p[i] + q[i] * delta3;
y[i] = laneSelect(active, y[i] + delta, y[i]);
change = laneMax(change, laneAbs(delta));
}
active = converged(active, change, tolerance);
}
integrate(s, y, 4);
return y[3];
}
Selbst wenn man von einer künstlerischen Intention ausgeht und es sich nicht einfach nur um Gebrauchs-/Verbrauchsmusik handelt, darf das Werkzeug imho auch mit K.I. Hilfe erstellt worden sein.
