NES palette functional, works

#include <cmath>
#include <algorithm>

// Scaled signal levels (see Level() for comments).
// At 12*8 are the original 8 signal levels.
// The first 12*8 are the same, except multiplied by 1.7
// for increased saturation, and multiplied by sin(phase).
static constexpr unsigned char Levels[] =
    "YYYYYYYY"             "TYgykx\206\206"       "PYq\221x\216\247\247"
    "NYu\231}\227\263\263" "PYq\221x\216\247\247" "TYgykx\206\206"
    "YYYYYYYY"             "^YK9G:,,"             "bYA!:#\013\013"
    "cY=\0305\033\000\000"  "bYA!:#\013\013"       "^YK9G:,,"
    "SYi n}\216\216";
static const double mul = 53.04827*12, add = 88.900051*12;

// Determine the momentary signal level on NES PPU
static constexpr unsigned char Level(unsigned pixel, unsigned phase, unsigned mode)
{
    // Voltage levels, relative to synch voltage:
    // Signal lows:   0.350, 0.518, 0.962, 1.550
    // Signal highs:  1.094, 1.506, 1.962, 1.962

    // Normalized:    -0.116, 0.00, 0.307, 0.715
    //                0.399, 0.684, 1.000, 1.000

    // The "add", "mul" values are chosen in such manner that
    // Levels[] produces the Normalized values divided by 12.

    return Levels[
        mode*8
      + ((pixel&0xF)<14 ? (pixel>>4)&3 : 1) // For colors 14..15, level 1 is forced.
      + 4*((pixel&0xF) <= 12 * (((pixel&0xF) + phase)%12 < 6))];
}

// Calculate the sum of parameters
constexpr unsigned Sum() { return 0; }
template<typename... T2> constexpr unsigned Sum(unsigned v, T2... values)
    { return v + Sum(values...); }

// This gets the NTSC signal total (Y, I or Q) for given pixel color,
// by summing the momentary signal values at given phases together.
template<typename... T2>
static constexpr unsigned GetC(unsigned mode, unsigned pixel, T2... phase)
{
    return Sum( Level(pixel, phase, mode<12 ? (mode+phase)%12 : mode)... );
}

// YIQ components for each NES color:
static constexpr unsigned Y(unsigned p) { return GetC(12, p, 0,1,2,3,4,5,6,7,8,9,10,11); }
static constexpr unsigned I(unsigned p) { return GetC(0,  p, 0,1,2,3,4,5,6,7,8,9,10,11); }
static constexpr unsigned Q(unsigned p) { return GetC(3,  p, 0,1,2,3,4,5,6,7,8,9,10,11); }

// YIQ to RGB matrix coefficients:
static constexpr double RY = 1, RI =  0.946882, RQ =  0.623557;
static constexpr double GY = 1, GI = -0.274788, GQ = -0.635691;
static constexpr double BY = 1, BI = -1.108545, BQ =  1.709007;

// RGB translations of the YIQ value:
static constexpr double R(unsigned p) { return (RY * Y(p) + RI * I(p) + RQ * Q(p) - add*(RY+RI+RQ)) / mul; }
static constexpr double G(unsigned p) { return (GY * Y(p) + GI * I(p) + GQ * Q(p) - add*(GY+GI+GQ)) / mul; }
static constexpr double B(unsigned p) { return (BY * Y(p) + BI * I(p) + BQ * Q(p) - add*(BY+BI+BQ)) / mul; }

// Utility functions for clamping
static constexpr double Clamp01(double a) { return a<0 ? 0 : (a>1 ? 1 : a); }

// Get the linear RGB color (not gamma-corrected) for given NES pixel color
constexpr unsigned MakeNEScolor(unsigned pixel)
{
    return 0x010000 * (int)(255 * Clamp01( R(pixel) ))
         + 0x000100 * (int)(255 * Clamp01( G(pixel) ))
         + 0x000001 * (int)(255 * Clamp01( B(pixel) ));
}

#include <cstdio>

template<int... values>
static void Colors()
{
    static const unsigned colors[sizeof...(values)] = { MakeNEScolor(values) ... };
    for(auto c: colors)
        std::printf("%06X\n", c);
}

int main()
{
    Colors<0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,
           32,33,34,35,36,37,38,38,39,40,41,42,43,44,45,46,47,48,49,
           50,51,52,53,54,55,56,57,58,59,60,61,62,63>();
}