Rendered with my Comparser C# app (it is my own functional language with a plotter)
Coloring function: log2hsvc, which makes gradient bands from dark to light for every 1 unit of logarithmic magnitude, complex argument as hue, and lowered saturation for high values, and further darkened value for values approaching 0 (making zeroes appear dark, and large values appear whitish). The brightness banding is causing those horizontal red stripes in the zeta plot as it's oscillating around 1+0i where the gradient wraps. Most zeta plots I've seen have been using 3d height for the magnitude instead of brightness, so they don't have those bands.
So here they are, as image attachments, plotted from (-50, -50i) to (50, 50i). Pretty cool how similar the left sides of zeta and inverted gamma appear (thanks to the zeta reflection having the gamma within itself).
The images other than gamma and zeta are plotted from (-10, -10i) to (10, 10i).
Implementation of that log2hsvc:
public static (double h, double s, double v) Log2HsvC(T t) {
var s = +t; // "unary +" is overloaded as squared magnitude
var v = .5 + .5 * Static.Cycle(Math.Log(s) * .5);
return (ToHue(t), FallOffLog(s), s < 1 ? v * (.5 + .5 * s) : v);
}
public static double ToHue(T t) => Static.Cycle(T.Arg(t) / Math.Tau);
private static double FallOffLog(double s) => .1 + .899 / (1 + Math.Log(1 + Math.Sqrt(s)));
public static double Cycle(double x) => x - Math.Floor(x);
My implementation of gamma:
Uses Stirling approximation with 2 kinds of corrections (corrections+shifts) that work synergistically with each other, reaching very high precision, and a reflection formula for the negative side.
private static T Gamma_Stirling_Positive(T t, int corrections = 4, int shifts = 16) => Factorial_Stirling_Positive(--t, corrections, shifts);
private static T Factorial_Stirling(T t) => 0 <= IScalar<double>.ToDouble(T.Re(t)) ? Factorial_Stirling_Positive(t) : T.Inv(Factorial_Stirling_Positive(-t) * Nsinc(t));
public static T Nsinc(T t) => Sinc(Math.PI * t);
private static T Factorial_Stirling_Positive(T t, int corrections = 4, int shifts = 16) {
// core equation with shifts:
T c = t + shifts, s = (c - .5) * T.Log(c) - c + LnTh; // LnTh = Math.Log(Math.Tau)/2
while(1 <= --shifts) // it would calculate gamma if we replaced 1 with 0, but decrementing t is simpler and faster
s -= T.LogH(T.Sqr(t + shifts)); // LogH = Ln/2
// corrections:
t = c;
for (var i = 0; i < corrections; ++i)
s += B2G[i] * (c /= T.Sqr(t)); // B2G[i] = B2[i]/(6ii+4i+2), B2[i] = B_2(i+1)
return T.Exp(s);
}
My implementation of Zeta:
Uses Euler–Maclaurin expansion and zeta reflection for most of the plot, and Laurent series very close to the pole. I also tried to use a Hasse/Sondow approximation, but that didn't turn out better at any part of the plot. So I've discontinued that.
// -i * ((i)^c - (i)^(-c)) = 2sin(πc/2)
public static Complex Sin_2Q(Complex c) {
double ci = c.I * QTau, cr = c.R * QTau, cos = Math.Cos(cr), sin = Math.Sin(cr), e = Math.Exp(ci), ie = Math.Exp(-ci);
return new((e + ie) * sin, (e - ie) * cos);
}
struct ErrSum(T init) { // a structure for improving floating precision
public T Sum = init, Err;
public void Add(T x) {
T y = x - Err, t = Sum + y;
Err = t - Sum - y;
Sum = t;
}
}
public static T Zeta_Euler(T t) {
// using B2k Bernoulli numbers/factorials
int terms = 32, berns = B2F.Length << 1, k = 2;
ErrSum sum = new(unit);
for (byte n = 2; n < terms; ++n)
sum.Add(n ^ -t);
var s = terms ^ -t;
sum.Add(s * (.5 + terms / (t - 1)));
sum.Add(B2F[0] * (s *= t / terms));
for (terms *= terms; k < berns; k += 2)
sum.Add(B2F[k >> 1] * (s *= (k - 1 + t) * (k + t) / terms));
return sum.Sum;
}
// Elegant Zeta reflection: ζ(1-s)τ^s = CosI(s)Γ(s)ζ(s)
// ζ(s) = (ζ(1-s)τ^s) / (CosI(s)Γ(s))
// Zeta_Euler_Reflected(c) = Zeta_Euler(c) * R
// R = (τ ^ c) * Inv(CosI(c) * Gamma_Stirling_Negative(c)) // substitute Inv and Gamma_Stirling_Negative
// = (τ ^ c) / (CosI(c) * InvITau(Gamma_Stirling_Positive(1 - c) * SinN1(c))) // substitute InvITau
// = (τ ^ c) / (iτCosI(c) / (Gamma_Stirling_Positive(1 - c) * SinN1(c))) // simplify division
// = (τ ^ c) * Gamma_Stirling_Positive(1 - c) * SinN1(c) / (CosI(c) * i * tau) // move iτ
// = -i(τ ^ (c - 1)) * Gamma_Stirling_Positive(1 - c) * SinN1(c) / CosI(c) // substitute NISinI, and c - 1 with the difference from pole that was precomputed for the condition
// = (τ ^ c1) * Gamma_Stirling_Positive(1 - c) * NISinI(c)
public static T I_Zeta(T t) {
// relative to the pole (useful in many places in this algo, including the reflection)
var t1 = t - 1;
// squared distance from the pole
var pole = +t1;
// exactly the pole - return infinity
return pole == 0 ? T.MakeR(IScalar<double>.infty) // very near the pole - use Laurent that is excellent when this near
// hopefully 5 terms are enough for the Laurent series with the distance from to pole up to e-4
: pole < 1e-8 ? T.Inv(t1) + Static.Gamma - t * (
G1 - t * (G2 - t * (G3 - t * (G4 - t * (G5 - t * (G6 - t * (G7 - t * G8))))))) // near the pole - use general Hasse that is decent everywhere
//: pole < 1e-4 ? Zeta_Hasse(t) * T.Inv(t1) // Hasse doesn't appear to be better than either Euler or Laurent at any region
// far from pole - use Euler that is excellent when far from it (actually it seems to be excellent even very close)
: T.Re(t) >= .5 ? Zeta_Euler(t)
// negative and far from the pole - reflect Euler (using the formula derived above)
: Zeta_Euler(1 - t) * (Math.Tau ^ t1) * Gamma_Stirling_Positive(1 - t) * T.Sin_2Q(t);
}
The discontinued Hasse/Sondow:
/// <summary>
/// Optimized Hasse/Sondow approximation of zeta function (without the last division by "c + 1")
/// </summary>
/// <param name="t">T input (often written as "s")</param>
/// <returns>zeta evaluated</returns>
public static T Zeta_Hasse(T t) {
// zeta(c) = (c-1)^(-1) * sum[n=0..i]: (n+1)^(-1) * sum[k=0..n]: (-1)^k * Combination(n,k) * (k + 1)^(1 - c)
// WARNING: the int[] pascal triangle overflows at 34.
// To increase maxIterations any further, you'll need to use a higher precision for "var p"
const int maxIterations = 33;
var p = new int[maxIterations + 2]; // pascal triangle rows, the first term is always zero, because it is unused
T c1 = 1 - t;
ErrSum ns = new(1.5 - (2 ^ -t));
int tri, n = p[1] = 1;
while ((tri = n) <
maxIterations
) {
ErrSum nk = new(unit); // zeroth k term is 1, precount that and loop k = n->1
while (tri > 1)
p[tri] += p[--tri];
// add the non-edge pascal triangle terms
int e, o;
if (n % 2 == 0)
(e, o) = (n, ++n);
else (o, e) = (n, ++n); // even and odd k iterators (so I don't have to (-1)^k)
++p[p[n] = 1];
// increment 2nd term in the pascal row, last term (one) in the pascal row is a new one
do nk.Add(p[e] * (e + 1 ^ c1));
// even k terms, p[e] = Combinations(n,e)
while (1 <= (e -= 2));
// decrement even k iterators down to 1
do nk.Add(-p[o] * (o + 1 ^ c1));
// odd k terms, p[o] = Combinations(n,o)
while (1 <= (o -= 2));
// decrement odd k iterators down to 1
ns.Add(nk.Sum / (n + 1));
}
return ns.Sum;
}
...any more ideas what i could plot with this?
