/* 2022-09-07 */

/* inspired by https://twitter.com/Marina_Costant/status/1560219380482539522 */

/* gcc -Wall -o norm_infinity norm_infinity.c -lm -O3 -ffast-math -fomit-frame-pointer -march=native -fopenmp
 * ./norm_infinity > XX
 * ffmpeg -f rawvideo -s 640x360 -r 25 -pix_fmt bgra -i XX -vf palettegen palette.png
 * ffmpeg -f rawvideo -s 640x360 -r 25 -pix_fmt bgra -i XX -i palette.png -filter_complex paletteuse movie.gif
 *
 * or:
 * ffmpeg -f rawvideo -s 640x360 -r 25 -pix_fmt bgra -i XX -plays 0 movie.apng
 */

/* not realtime at all, set delta = 0.05 or higher for faster but uglier */
/* or improve dig_t(): we look for the smallest t>1 for which pnorm() = 1,
 * if it exists of course
 * also classic lookup tables, caching, and so on, can be done
 */

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

/* linear interpolation, x is in [x0..x1] we want y in [y0..y1] such that
 * when x is x0 then y is y0 and when x is x1 then y is y1
 */
/* x0   y0
 * x    y
 * x1   y1
 * (x-x0) / (x1-x0) = (y-y0) / (y1-y0)
 * y = y0 + (x-x0) / (x1-x0) * (y1-y0)
 */
#define P3(x, x0, x1, y0, y1) ((double)(y0) + ((double)(x)-(double)(x0)) / ((double)(x1)-(double)(x0)) * ((double)(y1)-(double)(y0)))

#define rotate_x(x, y, z, angle) do {          \
  double y2 = y * cos(angle) - z * sin(angle); \
  double z2 = z * cos(angle) + y * sin(angle); \
  y = y2;                                      \
  z = z2;                                      \
} while (0)

#define rotate_y(x, y, z, angle) do {          \
  double z2 = z * cos(angle) - x * sin(angle); \
  double x2 = x * cos(angle) + z * sin(angle); \
  x = x2;                                      \
  z = z2;                                      \
} while (0)

#define rotate_z(x, y, z, angle) do {          \
  double x2 = x * cos(angle) - y * sin(angle); \
  double y2 = y * cos(angle) + x * sin(angle); \
  x = x2;                                      \
  y = y2;                                      \
} while (0)

#define clip(v, a, b) ((v) < (a) ? (a) : (v) > (b) ? (b) : (v))
#define max(a, b) ((a) > (b) ? (a) : (b))

/****************************************************************************/
/* bezier curve begin                                                       */
/****************************************************************************/

/* bezier curves are super nice for trajectories */

typedef struct {
  double x;
  double y;
  double z;
} point_t;

typedef struct {
  point_t p0;
  point_t p1;
  point_t p2;
  point_t p3;
} bezier_t;

/* cubic bezier (4 control points) */
point_t bezier(bezier_t *b, double t)
{
  double _1_t = 1 - t;
  double _1_t_2 = _1_t * _1_t;
  double _1_t_3 = _1_t_2 * _1_t;
  double t_2 = t * t;
  double t_3 = t * t_2;
  return (point_t){
    x: _1_t_3 * b->p0.x + 3 * _1_t_2 * t * b->p1.x + 3 * _1_t * t_2 * b->p2.x + t_3 * b->p3.x,
    y: _1_t_3 * b->p0.y + 3 * _1_t_2 * t * b->p1.y + 3 * _1_t * t_2 * b->p2.y + t_3 * b->p3.y,
    z: _1_t_3 * b->p0.z + 3 * _1_t_2 * t * b->p1.z + 3 * _1_t * t_2 * b->p2.z + t_3 * b->p3.z
  };
}

/****************************************************************************/
/* bezier curve end                                                         */
/****************************************************************************/

/****************************************************************************/
/* camera management begin                                                  */
/****************************************************************************/

typedef struct {
  point_t eye;           /* position of the eye */
  point_t center;        /* center of screen of camera */
  point_t vx, vy;        /* unit vectors on screen */
} camera_t;

point_t vectorial_product(point_t a, point_t b)
{
  /* a ^ b = ( a.y b.z - a.z b.y )
   *         ( a.z b.x - a.x b.z )
   *         ( a.x b.y - a.y b.x )
   */
  return (point_t){ x: a.y * b.z - a.z * b.y,
                    y: a.z * b.x - a.x * b.z,
                    z: a.x * b.y - a.y * b.x };
}

camera_t compute_camera(int frame)
{
  camera_t ret;
  double d;
  double t;
  double angle;

  frame %= 400;

  /* basic idea is to start the camera in front of cube (0,0,0) then move
   * forward with some rotations for nice visual effect and end in front
   * of cube (0,0,-24) so that we can loop the animation (the cube (0,0,-24)
   * has the same neighbors as the cube (0,0,0), so we can loop)
   */

  /* position the eye */
  /* eye movement is snakish, we start in front of cube (0,0,0) and end
   * in front of cube (0,0,-24) turning around center of each cube on the path
   */
  if (frame < 100) {
    /* turn on the left */
    angle = P3(frame, 0, 100, -M_PI/2, -3*M_PI/2);
  } else if (frame < 200) {
    /* right */
    angle = P3(frame, 100, 200, -M_PI/2, M_PI/2);
  } else if (frame < 300) {
    /* left */
    angle = P3(frame, 200, 300, -M_PI/2, -3*M_PI/2);
  } else {
    /* right */
    angle = P3(frame, 300, 400, -M_PI/2, M_PI/2);
  }
  /* we start with a point pointing to the right: (3,0,0) */
  point_t p0 = { x:3, y:0, z:0 };
  rotate_y(p0.x, p0.y, p0.z, angle);
  /* we force z to move forward, the rotation-based stuff is ugly for z */
  p0.z = P3(frame, 0, 100, 3, -3);
  /* move a bit closer to the shape */
  p0.z += - .5;
  /* tilt trajectory */
  rotate_z(p0.x, p0.y, p0.z, M_PI/6);

  ret.eye = p0;

  fprintf(stderr, "%d p0 %g %g %g\n", frame, p0.x, p0.y, p0.z);

  /* position p1 (direction of camera from the eye) */
  /* more or less going forward but some backward for not too much boring */
  bezier_t p1b = {
    p0: { x: 0,   y: 0,  z: 0 },     /* start point */
    p1: { x: 5,   y: 3,  z: 2 },
    p2: { x: -5,  y: -3, z: -46 },
    p3: { x: 0,   y: 0,  z: -24 },   /* end point */
  };
  point_t p1 = bezier(&p1b, P3(frame, 0, 400, 0, 1));

  fprintf(stderr, "   p1 %g %g %g\n", p1.x, p1.y, p1.z);

  /* center is on vector (p1-p0) at right distance of the eye */
  d = sqrt((p1.x-p0.x)*(p1.x-p0.x) + (p1.y-p0.y)*(p1.y-p0.y) + (p1.z-p0.z)*(p1.z-p0.z));
  /* right distance is 0.5 in front of eye (quite close, wide camera angle) */
  t = 0.5 / d;
  ret.center.x = p0.x + t * (p1.x - p0.x);
  ret.center.y = p0.y + t * (p1.y - p0.y);
  ret.center.z = p0.z + t * (p1.z - p0.z);

  p1 = ret.center;

  fprintf(stderr, "   center %g %g %g\n", ret.center.x, ret.center.y, ret.center.z);

  /* position p2 (anchor for angle of camera) (not sure if that means
   * something, but it's clear in my head: 'eye' is where is the eye, 'center',
   * is where the eye looks at, but then the camera can rotate around this
   * center, so we need a third point to define this rotation, this is p2)
   */
  /* rotate a bit, not too much, more or less randomly chosen */
  bezier_t p2b = {
    p0: { x: 10, y: 0,  z: 0 },      /* start: no rotation, horizontal cam */
    p1: { x: -10, y: -10, z: 0 },
    p2: { x: 10,  y: 10,  z: -20 },
    p3: { x: 10,  y: 0, z: -24 },    /* end: no rotation, horizontal cam */
  };
  point_t p2 = bezier(&p2b, P3(frame, 0, 400, 0, 1));

  fprintf(stderr, "   p2 %g %g %g\n", p2.x, p2.y, p2.z);

  /* project p2 on plane perpendicular to vector (p1-p0) containing p1 */
  /* vector (p2p-p2) and (p1-p0) are parallel: (p2p-p2) = t * (p1-p0)
   * and scalar product (p2p-p1) x (p1-p0) = 0
   * (p2 + t * (p1-p0) - p1) x (p1-p0) = 0
   * (p2.x + t * (p1.x-p0.x) - p1.x) * (p1.x - p0.x)
   * + [same y]
   * + [same z] = 0
   * ((p2.x - p1.x) + t * (p1.x-p0.x)) * (p1.x - p0.x)
   * + [same y]
   * + [same z] = 0
   * (p2.x - p1.x) * (p1.x - p0.x) + t * (p1.x-p0.x) * (p1.x - p0.x)
   * + [same y]
   * + [same z] = 0
   * t * (p1.x-p0.x) * (p1.x - p0.x)
   * + [same y]
   * + [same z] = -(p2.x - p1.x) * (p1.x - p0.x)
   * t = (-(p2.x - p1.x) * (p1.x - p0.x) -(p2.y - p1.y) * (p1.y - p0.y) -(p2.z - p1.z) * (p1.z - p0.z)) /
   *     ((p1.x-p0.x) * (p1.x - p0.x) + (p1.y-p0.y) * (p1.y - p0.y) + (p1.z-p0.z) * (p1.z - p0.z))
   */
  point_t p2p;
  t = (-(p2.x - p1.x) * (p1.x - p0.x) -(p2.y - p1.y) * (p1.y - p0.y) -(p2.z - p1.z) * (p1.z - p0.z)) /
      ((p1.x-p0.x) * (p1.x - p0.x) + (p1.y-p0.y) * (p1.y - p0.y) + (p1.z-p0.z) * (p1.z - p0.z));
  p2p.x = p2.x + t * (p1.x - p0.x);
  p2p.y = p2.y + t * (p1.y - p0.y);
  p2p.z = p2.z + t * (p1.z - p0.z);

  fprintf(stderr, "   p2p %g %g %g\n", p2p.x, p2p.y, p2p.z);

  /* vx is vector (p2p-p1) at right distance */
  d = sqrt((p2p.x-p1.x)*(p2p.x-p1.x) + (p2p.y-p1.y)*(p2p.y-p1.y) + (p2p.z-p1.z)*(p2p.z-p1.z));
  /* length (right distance) is 1/150 */
  t = 1./150 / d;
  ret.vx.x = t * (p2p.x - p1.x);
  ret.vx.y = t * (p2p.y - p1.y);
  ret.vx.z = t * (p2p.z - p1.z);

  /* vy is vectorial product of vx and vector (center-p0) */
  ret.vy = vectorial_product(ret.vx, (point_t){x: ret.center.x-p0.x,
                                               y: ret.center.y-p0.y,
                                               z: ret.center.z-p0.z});
  /* at the right length (also 1/150, as for vx) */
  d = sqrt(ret.vy.x*ret.vy.x + ret.vy.y*ret.vy.y + ret.vy.z*ret.vy.z);
  t = 1./150 / d;
  ret.vy.x *= t;
  ret.vy.y *= t;
  ret.vy.z *= t;

  fprintf(stderr, "    vx %g %g %g [%g]  vy %g %g %g [%g]\n",
          ret.vx.x, ret.vx.y, ret.vx.z,
          sqrt(ret.vx.x*ret.vx.x + ret.vx.y*ret.vx.y + ret.vx.z*ret.vx.z),
          ret.vy.x, ret.vy.y, ret.vy.z,
          sqrt(ret.vy.x*ret.vy.x + ret.vy.y*ret.vy.y + ret.vy.z*ret.vy.z));

  return ret;
}

/****************************************************************************/
/* camera management end                                                    */
/****************************************************************************/

typedef struct {
  unsigned char r, g, b;
} color_t;

unsigned char data[160*90*4];

double pnorm(double p, double t, double x0, double y0, double z0, double x1, double y1, double z1,
             int frame,
             double p2, double p3, double p4)
{
  /* the space is covered by cubes of size 6x6x6, repeating in all directions,
   * each one containing a shape.
   * We want 4 different shapes on screen for less boring, so each cube is
   * numbered in [0, 1, 2, 3] and the numbering repeats like this in all
   * directions, each numbering receives a shape (different shapes have
   * different p for a given frame)
   */

  /* (x,y,z) is the point to test */
  double x = x0 + t * (x1 - x0);
  double y = y0 + t * (y1 - y0);
  double z = z0 + t * (z1 - z0);

  /* find in which cube we are */
  int px = round(x / 6); px = (px % 4 + 4) % 4;
  int py = round(y / 6); py = (py % 4 + 4) % 4;
  int pz = round(z / 6); pz = (pz % 4 + 4) % 4;

  /* put (x,y,z) in the cube (0,0,0) for processing */
  x = fmod(fmod(x-3,6) + 6, 6)-3;
  y = fmod(fmod(y-3,6) + 6, 6)-3;
  z = fmod(fmod(z-3,6) + 6, 6)-3;

  /* rotate the shape */
  double angle = 2 * M_PI * frame / 100;
  rotate_y(x, y, z, angle);

  /* select p depending on the cube we're in */
  if ((px+py+pz) % 4 == 1) p = p2;
  else if ((px+py+pz) % 4 == 2) p = p3;
  else if ((px+py+pz) % 4 == 3) p = p4;

  /* no need for final pow(1/p), we care about distance = 1 */
  return   pow(fabs(x), p)
         + pow(fabs(y), p)
         + pow(fabs(z), p);
}

double get_p(int frame)
{
  double p;
  frame %= 200;
  /* we get p from [1/2 .. 4], linear choice is ugly, so we use several
   * linear ranges, looks more or less nice, there may be better choices
   */
  /* 1/2   1   2    4 */
  if (frame < 50)       p = P3(frame, 0, 50, 1./2, 1);
  else if (frame < 75)  p = P3(frame, 50, 75, 1, 2);
  else if (frame < 100) p = P3(frame, 75, 100, 2, 4);
  else if (frame < 125) p = P3(frame, 100, 125, 4, 2);
  else if (frame < 150) p = P3(frame, 125, 150, 2, 1);
  else                  p = P3(frame, 150, 200, 1, 1./2);
  return p;
}

double dig_t(double p, double x0, double y0, double z0, double x1, double y1, double z1, int frame, double p2, double p3, double p4)
{
  double t0, t1;
  double n0, n1;
  /* delta = .01 for nice result but very very slow, bigger value for faster
   * computation but uglier on screen
   */
  double delta = .01;

  /* norm is abs(x)^p + abs(y)^p + abs(z)^p (we omit final ^1/p)
   * we want norm = 1 (we want to be on the "unit sphere" of the norm)
   * current ray is from (x0, y0, z0) to (x1,y1,z1)
   * a point (x,y,z) is on the ray if there is t such that:
   * (x,y,z) = t * ray
   * for t=0 we decide to be at (x0,y0,z0)
   * for t=1 we decide to be at (x1,y1,z1)
   * so for generic point (x,y,z), we have:
   * 0 x0  y0  z0
   * 1 x1  y1  z1
   * t x   y   z
   * so:
   * t = (x - x0) / (x1 - x0)
   * same for y and z
   * and so:
   * x = x0 + t * (x1 - x0)
   * y = y0 + t * (y1 - y0)
   * z = z0 + t * (z1 - z0)
   */

  /* start from t = 1 to t = 70, step is delta (this is very slow) */
  t0 = 1;
  n0 = pnorm(p, t0, x0, y0, z0, x1, y1, z1, frame, p2, p3, p4);
  if (n0 <= 1) { fprintf(stderr, "bad start\n"); exit(1); }
  t1 = t0 + delta;
  n1 = pnorm(p, t1, x0, y0, z0, x1, y1, z1, frame, p2, p3, p4);
  /* we stop if we find a point with n = 1 or if we go too far */
  while ((   (n0 > 1 && n1 > 1)
          || (n0 < 1 && n1 < 1))
         && t1 < 70) {
    n0 = n1;
    t0 = t1;
    t1 += delta;
    n1 = pnorm(p, t1, x0, y0, z0, x1, y1, z1, frame, p2, p3, p4);
  }
  if (t1 >= 70) return 1000;
  /* this is not correct, t for which distance = 1 is not necessarily t1
   * but it looks ok on screen, so no big deal at this point
   */
  return t1;
}

color_t raytrace_pnorm(double x, double y, int frame, camera_t c)
{
  double x0, y0, z0, x1, y1, z1;
  double p;
  double t;

  /* ray is from eye to (x,y) on screen (center of screen is (0,0)) */
  /* start of ray: eye */
  x0 = c.eye.x;
  y0 = c.eye.y;
  z0 = c.eye.z;

  /* end of ray */
  /* (x,y) of screen is vector (x*vx + y*vx) starting at point center in
   * our 3d space
   */
  x1 = c.center.x + x * c.vx.x + y * c.vy.x;
  y1 = c.center.y + x * c.vx.y + y * c.vy.y;
  z1 = c.center.z + x * c.vx.z + y * c.vy.z;

  /* get the p in [1/2 .. 4], four of them for more diversity on screen */
  p = get_p(frame);
  double p2 = get_p(frame + 100);
  double p3 = get_p(frame + 50);
  double p4 = get_p(frame + 150);

  /* blue is actually black, looks better ('blue' name kept, lazy me) */
  color_t blue = { r : 0, g : 0, b : 0 };
  t = dig_t(p, x0, y0, z0, x1, y1, z1, frame, p2, p3, p4);
  if (t == 1000) return blue;

  /* a palette based on euclidian distance from (0,0,0) to the point */
  /* close is yellow, farther is red, far is white */
  double d = pnorm(2, t, x0, y0, z0, x1, y1, z1, frame, 2, 2, 2);
  /* g+rb define yellow */
  double g = P3(d, 0, 0.6, 125, 10);
  double rb = P3(d, 0, 0.6, 125, 10);
  /* ra defines red */
  double ra = P3(d, 0, 2, 0, 255);
  g = clip(g, 0, 255);
  rb = clip(rb, 0, 255);
  /* w defines white */
  double w = clip(P3(d, 1, 1.4, 0, 120), 0, 255);
  /* alpha is used to fade to 'blue' when very far away, looks nice */
  double alpha = clip(P3(t, 10, 70, 1, 0), 0, 1);
  return (color_t){
    r: alpha * max(clip(ra+rb, 0, 255), w) + (1-alpha) * blue.r,
    g : alpha * max(g, w) + (1-alpha) * blue.g,
    b : alpha * max(20, w) + (1-alpha) * blue.b
  };
}

color_t raytrace(int x, int y, int frame, camera_t c)
{
  /* center (x,y), reverse y for negative values to be at bottom not top */
  x -= 80;
  y = 45 - y;
  return raytrace_pnorm(x, y, frame, c);
}

void point(int x, int y, color_t color)
{
  data[(y * 160 + x) * 4] = color.b;
  data[(y * 160 + x) * 4 + 1] = color.g;
  data[(y * 160 + x) * 4 + 2] = color.r;
  data[(y * 160 + x) * 4 + 3] = 255;
}

void draw(int frame)
{
  camera_t c;
  int x, y;
  c = compute_camera(frame);
  for (y = 0; y < 90; y++) {
/* some openmp speed up (from > 2h computation I go down to 5 minutes on
 * some powerful computer with > 30 CPUs (shame on me, that thing should
 * be realtime on a single CPU...))
 */
#pragma omp parallel for
    for (x = 0; x < 160; x++) {
      point(x, y, raytrace(x, y, frame, c));
    }
  }
}

unsigned char bigscreen[160*90*4*4*4];

void dump(void)
{
  unsigned char *b = bigscreen;
  unsigned char *screen = data;
  int x, y;
  int i;
  /* initial image is 160x90, let's output something 640x360, so each pixel
   * becomes 4x4 instead of 1x1
   */
  for (y = 0; y < 90; y++) for (i = 0; i < 4; i++) for (x = 0; x < 160; x++) {
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
  }
  fwrite(bigscreen, 160*90*4*4*4, 1, stdout);
}

int main(void)
{
  int i;

//  while (1)
  for (i = 0; i < 400; i++) {
    fprintf(stderr, "%d\n", i);
    draw(i);
    dump();
  }
  fflush(stdout);
  return 0;
}