av1-like decoder

av1.c

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

#define e(a, d) for(a=0;a<d;a++)
int YUV[3][1 << 20], tmpyuv[2048], scratch[166], width;
int  p = 0, r = 1;//reader accu
int  lumaRatio = 0, satuRatio = 0;//decode default

void read8b() {
	int b = getchar();
	r = r << 8;
	p = (p << 8) + (b < 0 ? 0 : b);
}

int readat(int offset) {
	if (r < 256)read8b();
	int*spad = scratch + offset * 2;
	int F = spad[0] + spad[1] + 2;
	int u = r * (spad[0] + 1) / F;
	int k = p >= u;
	if (k) {
		p -= u;
		r -= u;
	} else r = u;
	spad[k]++;
	if (F > 63) {
		spad[0] /= 2;
		spad[1] /= 2;
	}
	return k;
}
// r = 256 , p = 0+23
int n(int step) {
	int a = 0;
	while (readat(step + a) == 0)a++;

	int b = 1;
	while (a--)b = b << 1 | readat(4);
	return b - 1;
}

//twice faster with precomputed cos array
void fft(int *S, int Y, const int *T, int Z, int U, int pass, int g) {
	int H; e(H, pass) {
		int a;e(a, pass) {
			int G = 0;
			int p;e(p, pass) {
				G += T[p * Z + H * U] * lrint(cos(3.1415 / pass * (a + 0.5) * p) * (p ? 1.414 : 1) * 1024);
			}
			S[a * Y + H * U] = G + (1 << g - 1) >> g;
		}
	}
}

void decode(int col, int line, int depth) {
	int a, k;
	if (depth > 5 || depth > 2 && readat(depth - 3)) {
		int c = 1 << --depth;
		e(a, 4)
			decode(col + a % 2 * c, line + a / 2 * c, depth);
		return;
	}
	int size = 1 << depth;
	int surface = size * size;
	int q = n(73);
	int colorspace = 0;e(colorspace, 3) {
		int*yuv = YUV[colorspace] + line * width + col;
		int isColor = colorspace > 0;
		e(a, surface)tmpyuv[a] = 0;
		e(a, surface) {
			if (readat(61 + depth * 2 + isColor))break;
			a += n(5 + isColor * 10);
			k = 1 - 2 * readat(3);
			tmpyuv[a] = k * (n(25 + (isColor + (a < surface / 8) * 2) * 10) + 1) * (colorspace ? satuRatio : lumaRatio);
		}
		if (!q) {
			int v = 0;
			e(a, size) {
				v += line ? yuv[a - width    ] : 0;
				v += col  ? yuv[a * width - 1] : 0;
			}
			*tmpyuv += col && line ? v / 2 : v;
		}
		fft(tmpyuv + surface,     1, tmpyuv          , 1   , size, size, 10);
		fft(yuv             , width, tmpyuv + surface, size,    1, size, 10 + depth);
		if (!q) continue;
		int C = q < 17;
		int w = C ? 9 - q : q - 25;
		e(a, size){
			int j;e(j, size) {
				int I,J;
				e(I, 2) {
					int h = a * w + w;
					J = h & 7;
					h = (h >> 3) + j + I;
					k = h < 0;
					if (k)h = (h * 8 + w / 2) / w - 2;
					h = h < size ? h : size - 1;
					tmpyuv[I] = k ^ C ? yuv[h * width - 1] : yuv[-width + h];
				}
				yuv[C ? j * width + a : a * width + j] += *tmpyuv * (8 - J) + tmpyuv[1] * J + 4 >> 3;
			}
		}
	}
}

void output(int b) {
	putchar(b < 0 ? 0 : b > 255 ? 255 : b);
}

int main(int argc) {
	read8b();// r = 256 , p = 0+23
	int depth = n(5);
	width = 1 << depth;
	int height = width - n(5);
	lumaRatio = n(5);
	satuRatio = n(5);
	decode(0, 0, depth);
	printf("P6 %d %d 255 ", width, height);
	for(int a=0,n=height * width;a < n;a++) {
		int l = YUV[0][a];//-33..283
		int L = YUV[1][a];//-53..43
		int M = YUV[2][a];//-79..122
		output(l - L + M);
		output(l + L);
		output(l - L - M);
	}
	return 0;
}