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

 /*
  * Do a finer sampling the higher up we get in the picture,
  * since there's much more detail there.
  */
 #define SAMPLES_FROM_Y(y) ((16384 * 512) / (y*y + 512))
 #define MAX_SAMPLES SAMPLES_FROM_Y(0)

 float weights[MAX_SAMPLES][MAX_SAMPLES];

 static inline int color(double x, double y)
 {
         double t, z;
         int i, j, k;

         x = x * (1.0 / 128.0) - 0.5;
         y = y * (1.0 / 2048.0);

         t = 1.0 / (y + 0.001);
         z = t * x;
         i = floor(t);
         j = floor(z);
         k = i + j;
         return ((k % 2) != 0);
 }

 static double sinc(double x)
 {
         static const double cutoff = 1.220703668e-4;  /* sqrt(sqrt(eps)) */

         if (abs(x) < cutoff) {
                 /* For small |x|, use Taylor series instead */
                 const double x2 = x * x;
                 const double x4 = x2 * x2;

                 return 1.0 - x2 / 6.0 + x4 / 120.0;
         } else {
                 return sin(x) / x;
         }
 }

 static double lanczos_tap(double x)
 {
         if (x < -3.0 || x > 3.0)
                 return 0.0;
         if (x < 0.0)
                 return sinc(-x*M_PI) * sinc(-x*M_PI / 3.0);
         else
                 return sinc(x*M_PI) * sinc(x*M_PI / 3.0);
 }

 static void precalculate_weights(unsigned num_samples)
 {
         double total = 0.0;
         for (int yi = 0; yi < num_samples; ++yi) {
                 double sy = 6.0 * ((double)yi / num_samples - 0.5);
                 double wy = lanczos_tap(sy);
                 for (int xi = 0; xi < num_samples; ++xi) {
                         double sx = 6.0 * ((double)xi / num_samples - 0.5);
                         weights[yi][xi] = wy * lanczos_tap(sx);
                         total += weights[yi][xi];
                 }
         }

         double inv_total = 1.0 / total;
         for (int yi = 0; yi < num_samples; ++yi) {
                 for (int xi = 0; xi < num_samples; ++xi) {
                         weights[yi][xi] *= inv_total;
                 }
         }
 }

 static double antialiased_color(double x, double y, unsigned num_samples)
 {
         double acc = 0.0;
         double delta = 6.0 / num_samples;
         int xi, yi;
         double sx, sy;

         for (yi = 0, sy = y - 3.0; yi < num_samples; ++yi, sy += delta) {
                 for (xi = 0, sx = x - 3.0; xi < num_samples; ++xi, sx += delta) {
                         if (color(sx, sy)) {
                                 acc += weights[yi][xi];
                         }
                 }
         }
         return acc;
 }