# Needlework

# 1 Video

# 2 Source Code

A variant on a program that calculates the number of components of the Mandelbrot set that are on the real axis (see OEIS sequences A000048 and A051841).

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

#include <mandelbrot-graphics.h>
#include <mandelbrot-numerics.h>
#include <mandelbrot-symbolics.h>

#include <sndfile.h>

int main(int argc, char **argv)
{
  (void) argc;
  (void) argv;
  int width = 640;
  int height = 360;
  int maxiters = 100100;
  double er = exp(2 * 3.141592653);
  m_pixel_t black = m_pixel_rgba(0,0,0,1);
  m_pixel_t grey = m_pixel_rgba(0.9,0.9,0.9,1);
  m_pixel_t white = m_pixel_rgba(1,1,1,1);
  m_d_colour_t *colour = m_d_colour_minimal_grid(white, black, white, grey);
  m_image *image = m_image_new(width, height);
  m_binangle b;
  mpq_t q, qq, one;
  m_binangle_init(&b);
  mpq_init(q);
  mpq_init(qq);
  mpq_init(one);
  mpq_set_ui(one, 1, 1);
  long long period = 17;
  long long mask = (1 << (period >> 1)) - 1;
  long long total = 0;
  long long real = 0;
  long long bulb = 0;
  long long den = (1ll << period) - 1;
  int frame = 0;
  float buffer[1600][2];
  float sintab[1600][period];
  float costab[1600][period];
  for (int i = 0; i < 1600; ++i)
  {
    float window = (1 - cos(2 * M_PI * (i + 0.5) / 1600));
    for (int j = 0; j < period; ++j)
    {
      sintab[i][j] = window * sin(2 * M_PI * (1 + j) * (i + 0.5) / 1600) / sqrt(1 + j);
      costab[i][j] = window * cos(2 * M_PI * (1 + j) * (i + 0.5) / 1600) / sqrt(1 + j);
    }
  }
  SF_INFO info = {0};
  info.samplerate = 48000;
  info.channels = 2;
  info.format = SF_FORMAT_WAV | SF_FORMAT_FLOAT;
  SNDFILE *out = sf_open("needlework.wav", SFM_WRITE, &info);
  for (long long num = 1; 2 * num < den; ++num)
  {
    b.per.length = period;
    mpz_set_si(b.per.bits, num);
    m_binangle_to_rational(q, &b);
    m_binangle_from_rational(&b, q);
    if (b.per.length == period)
    {
      total++;
      m_q_conjugate(qq, q);
      mpq_sub(qq, one, qq);
      if (mpq_equal(q, qq))
      {
        real++;
        if ((period & 1) == 0 && ((~(num >> (period >> 1))) & mask) == (num & mask))
        {
          bulb++;
        }
        else
        {
          double _Complex c = m_d_exray_in_do(q, 8, 2 * 8 * period, 16);
          m_d_nucleus(&c, c, period, 64);
          double r = 16 * cabs(m_d_size(c, period));
          m_d_transform *transform = m_d_transform_rectangular(width, height, c, r);
          m_d_render_scanline(image, transform, er, maxiters, colour);
          m_d_transform_delete(transform);
          char filename[100];
          snprintf(filename, sizeof(filename), "%04d.png", frame++);
          m_image_save_png(image, filename);
          memset(buffer, 0, sizeof(buffer));
          for (int i = 0; i < 1600; ++i)
          {
            for (int j = 0; j < period; ++j)
            {
              if ((num >> (period-1 - j)) & 1)
              {
                buffer[i][0] += sintab[i][j];
                buffer[i][1] += costab[i][j];
              }
            }
            buffer[i][0] = sin(buffer[i][0]);
            buffer[i][1] = sin(buffer[i][1]);
          }
          sf_writef_float(out, &buffer[0][0], 1600);
        }
      }
    }
  }
  sf_close(out);
  return 0;
}

# 3 Video Encoding

for pass in 1 2 3
do
  ffmpeg -r 30 -i %04d.png \
    -pix_fmt yuv420p -profile:v high -level:v 4.1 \
    -b:v 2M -bf:v 2 -g:v 30 -movflags +faststart -pass $pass \
    -y v.mp4
done
fdkaac -b 192k needlework.wav 
ffmpeg -i v.mp4 -i needlework.m4a \
  -codec:a copy -codec:v copy -movflags +faststart needlework.mp4
ffmpeg -ss 10 -i needlework.mp4 -frames:v 1 needlework.jpg