Run this code instantly in your browser: https://instacode.app/gist/ce1443125cdb152b81b7261fcb85c810

zeta-primes-approximation.js

// https://www.youtube.com/live/H7jdq0elNoY?si=NNHutyJiNx5E8IPz&t=2085
// it is not possible to link this gist on youtube, bcs their links policy,
// so linked this gist with YT video

import { canvas, onResize } from 'canvas';
import { mobius } from '@gkucmierz/utils@3.2.0';

const ctx = canvas.getContext('2d');

// --- Configuration Constants ---
const BG_COLOR = '#0f172a'; // Background color of the canvas (Slate 900)
const LINE_OPACITY = 0.7;  // Opacity of all plotted function lines (0.0 to 1.0)
const MAX_X = 150;
const MAX_Y = 0;           // Maximum limit on Y axis (set to 0 or -1 for auto-scaling)
const STEP_X = 20;         // Grid step size for X axis (set to 0 for auto-scaling)
const STEP_Y = 0;          // Grid step size for Y axis (set to 0 for auto-scaling)

// --- Zeta Zeros Configuration ---
const RIEMANN_ZEROS_COUNT = 50; // Number of non-trivial zeta zeros to use (1 to 100)

// Imaginary parts of the first 100 non-trivial zeros of Riemann Zeta function
const ZETA_ZEROS = [
  14.13472514, 21.02203964, 25.01085758, 30.42487613, 32.93506159,
  37.58617816, 40.91871901, 43.32707328, 48.00515088, 49.77383248,
  52.97032148, 56.44624770, 59.34704400, 60.83177852, 65.11254405,
  67.07981053, 69.54640171, 72.06715767, 75.70469070, 77.14484007,
  79.33737502, 82.91038085, 84.73549298, 87.42527461, 88.80911121,
  92.49189927, 94.65134404, 95.87063423, 98.83119422, 101.31785101,
  103.72553804, 105.44662305, 107.16861118, 111.02953554, 111.87465918,
  114.32022092, 116.22668032, 118.79075, 121.37015, 122.94683,
  124.25682, 126.07530, 127.51664, 129.57889, 131.08779,
  133.48661, 134.78637, 138.11687, 139.73620, 141.12739,
  143.11182, 146.00972, 147.42265, 150.05332, 151.35397,
  153.02476, 156.10863, 157.55332, 158.84993, 161.18930,
  163.03668, 165.53755, 167.18434, 169.10279, 171.36616,
  172.76643, 174.55960, 178.23073, 179.91720, 225.84588,
  228.63117, 229.36372, 231.25033, 231.98737, 233.40263,
  234.62841, 236.09340, 237.34000, 239.59600, 240.62343,
  242.33230, 243.67055, 245.59107, 247.92591, 250.00713,
  251.43105, 252.78203, 254.99120, 256.33132, 257.69750,
  259.33735, 261.25032, 262.53934, 264.91264, 266.52423,
  268.16335, 269.38356, 271.13196, 272.88057, 274.11730
];

// --- Mathematical Helper Functions ---

// Check if a number is prime
function isPrime(n) {
  if (n < 2) return false;
  for (let i = 2; i <= Math.sqrt(n); i++) {
    if (n % i === 0) return false;
  }
  return true;
}

// Prime counting function pi(x)
function pi(x) {
  let count = 0;
  for (let i = 2; i <= x; i++) {
    if (isPrime(i)) count++;
  }
  return count;
}

// Logarithmic integral function Li(x) approximated numerically
function Li(x) {
  if (x <= 1.0001) return 0;
  const Li2 = 1.04516378;
  if (x === 2) return Li2;
  
  let sum = 0;
  const steps = 500;
  const dt = (x - 2) / steps;
  for (let i = 0; i < steps; i++) {
    const t = 2 + i * dt + dt / 2;
    sum += 1 / Math.log(t);
  }
  return Li2 + sum * dt;
}

// Riemann's prime-counting approximation R(x)
function R(x) {
  if (x <= 1.0001) return 0;
  let sum = 0;
  for (let n = 1; n <= 10; n++) {
    const mu = mobius(n);
    if (mu === 0) continue;
    const power = Math.pow(x, 1 / n);
    sum += (mu / n) * Li(power);
  }
  return sum;
}

// Calculates 2 * Re( Li(x^z) ) for complex z = a + i*b
// Using the analytical approximation: Re( Li(x^z) ) ≈ Re( x^z / (z * ln x) )
function ReLi_complex(x, a, b) {
  if (x <= 1.01) return 0; // Avoid near-singularity at x = 1
  const lnX = Math.log(x);
  const denom = a * a + b * b;
  return (Math.pow(x, a) / (lnX * denom)) * (a * Math.cos(b * lnX) + b * Math.sin(b * lnX));
}

// Calculates the contribution of a single zero pair to the Riemann prime-counting formula
function R_zero_pair(x, gamma, maxN = 5) {
  let sum = 0;
  for (let n = 1; n <= maxN; n++) {
    const mu = mobius(n);
    if (mu === 0) continue;
    // For a complex zero rho_k / n = (1/2n) + i * (gamma/n)
    const term = (mu / n) * 2 * ReLi_complex(x, 0.5 / n, gamma / n);
    sum += term;
  }
  return sum;
}

// Riemann's explicit prime-counting formula approximation
function piExplicit(x, numZeros) {
  if (x <= 1.5) return 0;
  let sum = R(x);
  // Subtract the oscillatory terms from nontrivial zeros
  const limit = Math.min(numZeros, ZETA_ZEROS.length);
  for (let k = 0; k < limit; k++) {
    sum -= R_zero_pair(x, ZETA_ZEROS[k], 5);
  }
  return sum;
}

// --- Canvas Rendering logic ---

onResize((width, height) => {
  canvas.width = width;
  canvas.height = height;

  // Clear background
  ctx.fillStyle = BG_COLOR;
  ctx.fillRect(0, 0, width, height);

  // Layout dimensions
  const paddingLeft = 50;
  const paddingRight = 130; // Extra space for the legend on the right
  const paddingTop = 65;
  const paddingBottom = 45;

  const graphWidth = width - paddingLeft - paddingRight;
  const graphHeight = height - paddingTop - paddingBottom;

  if (graphWidth <= 0 || graphHeight <= 0) return;

  // --- Dynamic Axis Calibration ---

  // 1. Calculate MAX_Y automatically if needed
  let maxYValue = MAX_Y;
  if (maxYValue <= 0) {
    const computedMaxY = Math.max(MAX_X / Math.log(MAX_X), Li(MAX_X), pi(MAX_X));
    maxYValue = Math.ceil(computedMaxY * 1.1);
  }

  // 2. Calculate STEP_X automatically if needed
  let stepXValue = STEP_X;
  if (stepXValue <= 0) {
    const rawStep = MAX_X / 5;
    if (rawStep <= 1) stepXValue = 1;
    else if (rawStep <= 2) stepXValue = 2;
    else if (rawStep <= 5) stepXValue = 5;
    else if (rawStep <= 10) stepXValue = 10;
    else if (rawStep <= 20) stepXValue = 20;
    else if (rawStep <= 50) stepXValue = 50;
    else if (rawStep <= 100) stepXValue = 100;
    else stepXValue = Math.ceil(rawStep / 50) * 50;
  }

  // 3. Calculate STEP_Y automatically if needed
  let stepYValue = STEP_Y;
  if (stepYValue <= 0) {
    const rawStep = maxYValue / 5;
    if (rawStep <= 1) stepYValue = 1;
    else if (rawStep <= 2) stepYValue = 2;
    else if (rawStep <= 5) stepYValue = 5;
    else if (rawStep <= 10) stepYValue = 10;
    else if (rawStep <= 20) stepYValue = 20;
    else if (rawStep <= 50) stepYValue = 50;
    else if (rawStep <= 100) stepYValue = 100;
    else stepYValue = Math.ceil(rawStep / 50) * 50;
  }

  // Coordinate conversion helpers
  const toScreenX = (valX) => paddingLeft + (valX / MAX_X) * graphWidth;
  const toScreenY = (valY) => height - paddingBottom - (valY / maxYValue) * graphHeight;

  // 1. Draw Title
  ctx.fillStyle = '#f8fafc';
  ctx.font = 'bold 15px sans-serif';
  ctx.textAlign = 'left';
  ctx.fillText('Primes & Riemann Zeta Zeros Relation', paddingLeft, 30);

  // Draw Nontrivial Zeros Subtitle
  ctx.fillStyle = '#94a3b8';
  ctx.font = '12px monospace';
  ctx.fillText(`${RIEMANN_ZEROS_COUNT} nontrivial zeros`, paddingLeft, 50);

  // 2. Draw Axes & Grid
  ctx.strokeStyle = '#334155'; // Slate 700
  ctx.lineWidth = 1;
  
  // X grid and labels
  ctx.fillStyle = '#94a3b8';
  ctx.font = '10px monospace';
  ctx.textAlign = 'center';
  ctx.textBaseline = 'top';
  for (let x = 0; x <= MAX_X; x += stepXValue) {
    const sx = toScreenX(x);
    // Grid line
    ctx.strokeStyle = x === 0 ? '#64748b' : '#1e293b';
    ctx.beginPath();
    ctx.moveTo(sx, toScreenY(0));
    ctx.lineTo(sx, toScreenY(maxYValue));
    ctx.stroke();
    // Label
    ctx.fillStyle = '#94a3b8';
    ctx.fillText(String(x), sx, toScreenY(0) + 6);
  }

  // Y grid and labels
  ctx.textAlign = 'right';
  ctx.textBaseline = 'middle';
  for (let y = 0; y <= maxYValue; y += stepYValue) {
    const sy = toScreenY(y);
    // Grid line
    ctx.strokeStyle = y === 0 ? '#64748b' : '#1e293b';
    ctx.beginPath();
    ctx.moveTo(toScreenX(0), sy);
    ctx.lineTo(toScreenX(MAX_X), sy);
    ctx.stroke();
    // Label
    ctx.fillStyle = '#94a3b8';
    ctx.fillText(String(y), toScreenX(0) - 8, sy);
  }

  // 3. Draw Functions (with global opacity)
  ctx.globalAlpha = LINE_OPACITY;

  // A. Gauss approximation: x / ln(x) (Green Line)
  ctx.strokeStyle = '#22c55e'; // Green 500
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  let first = true;
  for (let x = 2.0; x <= MAX_X; x += 0.5) {
    const y = x / Math.log(x);
    const sx = toScreenX(x);
    const sy = toScreenY(y);
    if (first) {
      ctx.moveTo(sx, sy);
      first = false;
    } else {
      ctx.lineTo(sx, sy);
    }
  }
  ctx.stroke();

  // B. Logarithmic integral: Li(x) (Purple Line)
  ctx.strokeStyle = '#d946ef'; // Fuchsia 500
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  first = true;
  for (let x = 2.0; x <= MAX_X; x += 0.5) {
    const y = Li(x);
    const sx = toScreenX(x);
    const sy = toScreenY(y);
    if (first) {
      ctx.moveTo(sx, sy);
      first = false;
    } else {
      ctx.lineTo(sx, sy);
    }
  }
  ctx.stroke();

  // C. Prime counting function: pi(x) (Red Step-like Line)
  ctx.strokeStyle = '#ef4444'; // Red 500
  ctx.lineWidth = 2.5;
  ctx.beginPath();
  ctx.moveTo(toScreenX(0), toScreenY(0));
  
  let currentY = 0;
  for (let x = 1; x <= MAX_X; x++) {
    const nextY = pi(x);
    const sxCurrent = toScreenX(x);
    const syCurrent = toScreenY(currentY);
    const syNext = toScreenY(nextY);
    
    // Draw horizontal step to current X, then vertical jump to next Y value
    ctx.lineTo(sxCurrent, syCurrent);
    ctx.lineTo(sxCurrent, syNext);
    currentY = nextY;
  }
  ctx.stroke();

  // D. Riemann Explicit Formula with Zeros: piExplicit(x) (Blue Line - drawn last to be on top)
  ctx.strokeStyle = '#3b82f6'; // Blue 500
  ctx.lineWidth = 2.2;
  ctx.beginPath();
  first = true;
  for (let x = 2.0; x <= MAX_X; x += 0.25) { // Higher resolution for oscillations
    const y = piExplicit(x, RIEMANN_ZEROS_COUNT);
    const sx = toScreenX(x);
    const sy = toScreenY(y);
    if (first) {
      ctx.moveTo(sx, sy);
      first = false;
    } else {
      ctx.lineTo(sx, sy);
    }
  }
  ctx.stroke();

  // Reset opacity for text, UI borders, and legend
  ctx.globalAlpha = 1.0;

  // 5. Draw Legend Box on the Right
  const legendX = toScreenX(MAX_X) + 15;
  const legendY = paddingTop + 10;
  const itemHeight = 22;

  ctx.fillStyle = '#1e293b'; // Slate 800
  ctx.strokeStyle = '#334155';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.roundRect(legendX - 5, legendY - 8, 105, 100, 4);
  ctx.fill();
  ctx.stroke();

  const legendItems = [
    { label: 'π(x)', color: '#ef4444' },
    { label: 'x/ln(x)', color: '#22c55e' },
    { label: 'Li(x)', color: '#d946ef' },
    { label: 'Riemann', color: '#3b82f6' }
  ];

  ctx.textAlign = 'left';
  ctx.textBaseline = 'middle';
  ctx.font = '11px "JetBrains Mono", monospace';

  legendItems.forEach((item, idx) => {
    const cy = legendY + idx * itemHeight;
    // Draw color line indicator
    ctx.strokeStyle = item.color;
    ctx.lineWidth = 3;
    ctx.beginPath();
    ctx.moveTo(legendX, cy);
    ctx.lineTo(legendX + 15, cy);
    ctx.stroke();
    // Draw label
    ctx.fillStyle = '#cbd5e1';
    ctx.fillText(item.label, legendX + 22, cy);
  });
});