Mark rational points on a sphere

MX914152b.cc

// http://math.stackexchange.com/a/914871/35416

#include <algorithm>
#include <cassert>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <vector>

constexpr int size = 600;
constexpr int size1 = size - 1;
constexpr int size2 = size*size;
constexpr int steps = 6;
constexpr int limit = 1500;

char img[size2];
unsigned short count[size2];

struct Quad {
  unsigned short d, a, b, c;
  Quad(unsigned short d, unsigned short a, unsigned short b, unsigned short c)
    : d(d), a(a), b(b), c(c) { }
  bool operator<(const Quad& q) const { return d < q.d; }
};

std::vector<Quad> quads;

inline int gcd(int a, int b) {
  while (b) {
    int c = a%b;
    a = b;
    b = c;
  }
  return a;
}

void gen() {
  // Following http://dx.doi.org/10.2307/2312125
  int du, dv, dw, d;
  for (int u = 1; (du = u*u) <= limit; ++u) { // (c1)
    for (int v = 0; (dv = du + v*v) <= limit; ++v) { // (c1)
      for (int w = 1; (dw = dv + w*w) <= limit; ++w) { // (c2)
        for (int t = (v ? 0 : 1); (d = dw + t*t) <= limit; ++t) { // (c2,c3)
          if (((u ^ v ^ w ^ t) & 1) != 1) continue; // (d)
          if ((t == 0 && u > v) || (v == 0 && w > t)) continue; // (f,g)
          int c = 2*dv - d;
          if (c <= 0) break; // (b)
          int a = 2*(u*w - v*t);
          if (a <= 0) break; // (a)
          int b = 2*(u*t + v*w);
          if (gcd(u*t + v*w, gcd(dv, d)) != 1) continue; // (e)
          quads.emplace_back(d, a, b, c);
          if (d % 2 == 0) {
            std::cout << "u=" << u << ", v=" << v << ", w=" << w << ", t=" << t 
                      << ", a=" << a << ", b=" << b << ", c=" << c
                      << ", d=" << d << std::endl;
            std::exit(1);
          }
        }
      }
    }
  }

  // Since the above generates only positive points, let's generate
  // points on the circle as well.
  for (int u = 0; (du = u*u) <= limit; ++u) {
    for (int v = u + 1; (d = du + v*v) <= limit; v += 2) {
      if (gcd(u, v) != 1) continue;
      quads.emplace_back(d, d - 2*du, 2*u*v, 0);
    }
  }
  std::sort(quads.begin(), quads.end());
}

void write(unsigned short d) {
  std::ostringstream oss;
  oss << std::setfill('0') << std::setw(8) << d << ".pgm";
  std::ofstream f(oss.str());
  f << "P5\n" << size << " " << size << "\n255\n";
  f.write(img, size2);
}

void write() {
  std::fill_n(img, size2, '\xff');
  std::fill_n(count, size2, 0);
  unsigned short d = 0;
  for (const Quad& q : quads) {
    if (q.d != d) {
      if (d) write(d);
      d = q.d;
      if (d % 2 == 0) {
        std::cout << ", a=" << q.a << ", b=" << q.b << ", c=" << q.c
                  << ", d=" << q.d << std::endl;
        std::exit(1);
      }
    }
    unsigned short abc[] = { q.a, q.b, q.c };
    std::sort(abc, abc + 3);
    if (false)
      std::cout << std::setw(3) << abc[0] << "² + "
                << std::setw(3) << abc[1] << "² + "
                << std::setw(3) << abc[2] << "² = "
                << std::setw(3) << q.d << "²\n";
    do {
      auto x = (size*abc[0])/d;
      auto y = (size*abc[1])/d;
      if (x == size) x = size - 1;
      if (y == size) y = size - 1;
      auto xy = size*(size - y - 1) + x;
      auto c = ++count[xy];
      if (c > steps) c = steps;
      img[xy] = 255 - (255*c + steps/2)/steps;
    } while (std::next_permutation(abc, abc + 3));
  }
  write(d);
}

int main() {
  gen();
  std::cout << "Generated pythagorean quadruples." << std::endl;
  write();
  return 0;
}

MX914152b.sh

#!/bin/bash

imgs=( )
for ((i=1; i<1500; i+=2)); do
    input=$(printf "%08d.pgm" $i)
    cnt=$(printf "%04d" $i)
    # test -e $cnt.png || \
    convert $input -family "DejaVu Sans" -pointsize 15 \
        -annotate +300+20 "d ≤ $i" $cnt.png || exit $?
    imgs+=( $cnt.png )
done
echo "Frames converted, will now convert whole image."
convert 00001499.pgm -colors 7 map.gif || exit $?
convert -delay 30 ${imgs[@]} +dither -map map.gif gif:- \
| convert gif:- -layers OptimizeFrame gif:- \
| convert gif:- -layers OptimizeTransparency gif:- \
| convert gif:- +map MX914152b.gif