It's cursed: https://www.desmos.com/calculator/asjxggbglo. What if you had only one discrete variable that you could plug into a sum sign, and *everything else had to be perfectly continuous* - no floor, no mod, no tricks. Could you even DO a generalized 2d riemann sum that converges as k -> infinity. Find out now (spoiler: you kinda can) - modi…

cursed-convergence.c

#include <stdio.h>
#include <math.h>
#include <time.h>

#define PI 3.14159265358979323846

double s_inv(double x) {
    return asin(2.0 * (x - 0.5)) / PI + 0.5;
}

double r(double x, double s) {
    return tanh(s * sin(2.0 * PI * (x - 0.25)));
}

double s_awtooth(double x, double s) {
    double inner = sin(PI * (x - 0.5) * r((x - 0.5) / 2.0, s)) / 2.0 + 0.5;
    return s_inv(inner);
}

double f_loor(double x, double s) {
    return x + s_awtooth(x, s) - 1.0;
}

double m_od(double x, double n, double s) {
    return x - n * f_loor(x/n, s);
}

double x_p(double n, double k) {
    return (m_od(n, k + 1, 2*k) - 1) / (k - 1) * 2.0 + (-1.0);
}

double y_p(double n, double k) {
    return f_loor(n / (k + 1), k) / k * 2.0 + (-1.0);
}

double C(double x, double y) {
    return x*x + y*y - 1.0;
}

double A(double x, double k) {
    return 1.0 / (1.0 + exp(-k * (-x)));
}

double V(double k) {
    double total = 0.0;
    const double A_bbox = 4.0;
    long n_max = (long)((k + 1) * (k + 1));
    
    for(long n = 0; n < n_max; n++) {
        double x = x_p(n, k);
        double y = y_p(n, k);
        total += A(C(x, y), k);
    }
    
    return (total / ((k + 1)*(k + 1))) * A_bbox;
}

int main() {
    // Print CSV header
    printf("k,area,error,time_ms\n");
    
    // Calculate for logarithmically spaced k values
    const int num_k = 100;
    double k_min = 10;
    double k_max = 30000.0;
    double log_k_min = log(k_min);
    double log_k_max = log(k_max);
    
    for(int i = 0; i < num_k; i++) {
        // Calculate k value
        double log_k = log_k_min + (log_k_max - log_k_min) * i / (num_k - 1);
        long k = (long)round(exp(log_k));
        
        // Time the calculation
        clock_t start = clock();
        double area = V(k);
        clock_t end = clock();
        
        double time_ms = 1000.0 * (end - start) / CLOCKS_PER_SEC;
        double error = fabs(area - PI);
        
        // Output CSV row
        printf("%ld,%.10f,%.10f,%.2f\n", k, area, error, time_ms);
        fflush(stdout);  // Ensure immediate output
    }
    
    return 0;
}