diff --git a/README.md b/README.md
index 044ed00..7c09c0a 100644
--- a/README.md
+++ b/README.md
@@ -84,6 +84,27 @@ Encoder for [augmented Hadamard codes](https://en.wikipedia.org/wiki/Hadamard_co
 
 [SIMD](https://en.wikipedia.org/wiki/SIMD) intra-frame accelerated [Low-density parity-check](https://en.wikipedia.org/wiki/Low-density_parity-check_code) layered decoder.
 
+### [polar_freezer.hh](polar_freezer.hh)
+
+Bit freezers for the construction of [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
+
+* PolarFreezer: Constructs code for a given erasure probability without the need for storing nor sorting of erasure probabilities for the virtual channels.
+* PolarCodeConst0: Constructs code by choosing the K best virtual channels computed from a given erasure probability.
+
+### [polar_encoder.hh](polar_encoder.hh)
+
+Encoders for [non-systematic and systematic](https://en.wikipedia.org/wiki/Systematic_code) [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
+
+### [polar_decoder.hh](polar_decoder.hh)
+
+Successive cancellation decoding of [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
+
+### [polar_list_decoder.hh](polar_list_decoder.hh)
+
+Successive cancellation [list decoding](https://en.wikipedia.org/wiki/List_decoding) of [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
+
+List size depends on used SIMD type. Decoding performance of the fixed-point implementation is better with shorter codes, as the list size is larger but a deterioration can be observed with larger codes.
+
 ### [osd.hh](osd.hh)
 
 Ordered statistics decoding allows the practical [soft-decision decoding](https://en.wikipedia.org/wiki/Soft-decision_decoder) of short to medium sized [linear codes](https://en.wikipedia.org/wiki/Linear_code).
diff --git a/polar_decoder.hh b/polar_decoder.hh
new file mode 100644
index 0000000..206263d
--- /dev/null
+++ b/polar_decoder.hh
@@ -0,0 +1,99 @@
+/*
+Successive cancellation decoding of polar codes
+
+Copyright 2020 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include <algorithm>
+#include "polar_helper.hh"
+
+namespace CODE {
+
+template <typename TYPE, int M>
+struct PolarTree
+{
+	typedef PolarHelper<TYPE> PH;
+	static const int N = 1 << M;
+	static void decode(TYPE **message, TYPE *hard, TYPE *soft, const uint8_t *frozen)
+	{
+		for (int i = 0; i < N/2; ++i)
+			soft[i+N/2] = PH::prod(soft[i+N], soft[i+N/2+N]);
+		PolarTree<TYPE, M-1>::decode(message, hard, soft, frozen);
+		for (int i = 0; i < N/2; ++i)
+			soft[i+N/2] = PH::madd(hard[i], soft[i+N], soft[i+N/2+N]);
+		PolarTree<TYPE, M-1>::decode(message, hard+N/2, soft, frozen+N/2);
+		for (int i = 0; i < N/2; ++i)
+			hard[i] = PH::qmul(hard[i], hard[i+N/2]);
+	}
+};
+
+template <typename TYPE>
+struct PolarTree<TYPE, 0>
+{
+	typedef PolarHelper<TYPE> PH;
+	static void decode(TYPE **message, TYPE *hard, TYPE *soft, const uint8_t *frozen)
+	{
+		if (*frozen) {
+			*hard = PH::one();
+		} else {
+			*hard = PH::signum(soft[1]);
+			*(*message)++ = *hard;
+		}
+	}
+};
+
+template <typename TYPE, int MAX_M>
+class PolarDecoder
+{
+	typedef PolarHelper<TYPE> PH;
+	static const int MAX_N = 1 << MAX_M;
+	TYPE soft[2*MAX_N];
+	TYPE hard[MAX_N];
+public:
+	void operator()(TYPE *message, const TYPE *codeword, const uint8_t *frozen, int level)
+	{
+		assert(level <= MAX_M);
+		int length = 1 << level;
+		for (int i = 0; i < length; ++i)
+			soft[length+i] = codeword[i];
+
+		switch (level) {
+		case 0: PolarTree<TYPE, 0>::decode(&message, hard, soft, frozen); break;
+		case 1: PolarTree<TYPE, 1>::decode(&message, hard, soft, frozen); break;
+		case 2: PolarTree<TYPE, 2>::decode(&message, hard, soft, frozen); break;
+		case 3: PolarTree<TYPE, 3>::decode(&message, hard, soft, frozen); break;
+		case 4: PolarTree<TYPE, 4>::decode(&message, hard, soft, frozen); break;
+		case 5: PolarTree<TYPE, 5>::decode(&message, hard, soft, frozen); break;
+		case 6: PolarTree<TYPE, 6>::decode(&message, hard, soft, frozen); break;
+		case 7: PolarTree<TYPE, 7>::decode(&message, hard, soft, frozen); break;
+		case 8: PolarTree<TYPE, 8>::decode(&message, hard, soft, frozen); break;
+		case 9: PolarTree<TYPE, 9>::decode(&message, hard, soft, frozen); break;
+		case 10: PolarTree<TYPE, 10>::decode(&message, hard, soft, frozen); break;
+		case 11: PolarTree<TYPE, 11>::decode(&message, hard, soft, frozen); break;
+		case 12: PolarTree<TYPE, 12>::decode(&message, hard, soft, frozen); break;
+		case 13: PolarTree<TYPE, 13>::decode(&message, hard, soft, frozen); break;
+		case 14: PolarTree<TYPE, 14>::decode(&message, hard, soft, frozen); break;
+		case 15: PolarTree<TYPE, 15>::decode(&message, hard, soft, frozen); break;
+		case 16: PolarTree<TYPE, 16>::decode(&message, hard, soft, frozen); break;
+		case 17: PolarTree<TYPE, 17>::decode(&message, hard, soft, frozen); break;
+		case 18: PolarTree<TYPE, 18>::decode(&message, hard, soft, frozen); break;
+		case 19: PolarTree<TYPE, 19>::decode(&message, hard, soft, frozen); break;
+		case 20: PolarTree<TYPE, 20>::decode(&message, hard, soft, frozen); break;
+		case 21: PolarTree<TYPE, 21>::decode(&message, hard, soft, frozen); break;
+		case 22: PolarTree<TYPE, 22>::decode(&message, hard, soft, frozen); break;
+		case 23: PolarTree<TYPE, 23>::decode(&message, hard, soft, frozen); break;
+		case 24: PolarTree<TYPE, 24>::decode(&message, hard, soft, frozen); break;
+		case 25: PolarTree<TYPE, 25>::decode(&message, hard, soft, frozen); break;
+		case 26: PolarTree<TYPE, 26>::decode(&message, hard, soft, frozen); break;
+		case 27: PolarTree<TYPE, 27>::decode(&message, hard, soft, frozen); break;
+		case 28: PolarTree<TYPE, 28>::decode(&message, hard, soft, frozen); break;
+		case 29: PolarTree<TYPE, 29>::decode(&message, hard, soft, frozen); break;
+		default: assert(false);
+		}
+	}
+};
+
+}
+
diff --git a/polar_encoder.hh b/polar_encoder.hh
new file mode 100644
index 0000000..8eb27b0
--- /dev/null
+++ b/polar_encoder.hh
@@ -0,0 +1,66 @@
+/*
+Polar encoder for non-systematic and systematic codes
+
+Copyright 2020 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include "polar_helper.hh"
+
+namespace CODE {
+
+template <typename TYPE>
+class PolarEncoder
+{
+	typedef PolarHelper<TYPE> PH;
+public:
+	void operator()(TYPE *codeword, const TYPE *message, const uint8_t *frozen, int level)
+	{
+		int length = 1 << level;
+		for (int i = 0; i < length; i += 2) {
+			TYPE msg0 = frozen[i] ? PH::one() : *message++;
+			TYPE msg1 = frozen[i+1] ? PH::one() : *message++;
+			codeword[i] = PH::qmul(msg0, msg1);
+			codeword[i+1] = msg1;
+		}
+		for (int h = 2; h < length; h *= 2)
+			for (int i = 0; i < length; i += 2 * h)
+				for (int j = i; j < i + h; ++j)
+					codeword[j] = PH::qmul(codeword[j], codeword[j+h]);
+	}
+};
+
+template <typename TYPE>
+class PolarSysEnc
+{
+	typedef PolarHelper<TYPE> PH;
+public:
+	void operator()(TYPE *codeword, const TYPE *message, const uint8_t *frozen, int level)
+	{
+		int length = 1 << level;
+		for (int i = 0; i < length; i += 2) {
+			TYPE msg0 = frozen[i] ? PH::one() : *message++;
+			TYPE msg1 = frozen[i+1] ? PH::one() : *message++;
+			codeword[i] = PH::qmul(msg0, msg1);
+			codeword[i+1] = msg1;
+		}
+		for (int h = 2; h < length; h *= 2)
+			for (int i = 0; i < length; i += 2 * h)
+				for (int j = i; j < i + h; ++j)
+					codeword[j] = PH::qmul(codeword[j], codeword[j+h]);
+		for (int i = 0; i < length; i += 2) {
+			TYPE msg0 = frozen[i] ? PH::one() : codeword[i];
+			TYPE msg1 = frozen[i+1] ? PH::one() : codeword[i+1];
+			codeword[i] = PH::qmul(msg0, msg1);
+			codeword[i+1] = msg1;
+		}
+		for (int h = 2; h < length; h *= 2)
+			for (int i = 0; i < length; i += 2 * h)
+				for (int j = i; j < i + h; ++j)
+					codeword[j] = PH::qmul(codeword[j], codeword[j+h]);
+	}
+};
+
+}
+
diff --git a/polar_freezer.hh b/polar_freezer.hh
new file mode 100644
index 0000000..35a1a09
--- /dev/null
+++ b/polar_freezer.hh
@@ -0,0 +1,67 @@
+/*
+Bit freezers for polar codes
+
+Copyright 2020 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include <algorithm>
+
+namespace CODE {
+
+class PolarFreezer
+{
+	static void freeze(uint8_t *bits, long double pe, long double th, int i, int h)
+	{
+		if (h) {
+			freeze(bits, pe * (2-pe), th, i, h/2);
+			freeze(bits, pe * pe, th, i+h, h/2);
+		} else {
+			bits[i] = pe > th;
+		}
+	}
+public:
+	int operator()(uint8_t *frozen_bits, int level, long double erasure_probability = 0.5L, long double freezing_threshold = 0.5L)
+	{
+		int length = 1 << level;
+		freeze(frozen_bits, erasure_probability, freezing_threshold, 0, length / 2);
+		int K = 0;
+		for (int i = 0; i < length; ++i)
+			K += !frozen_bits[i];
+		return K;
+	}
+};
+
+template <int MAX_M>
+class PolarCodeConst0
+{
+	void compute(long double pe, int i, int h)
+	{
+		if (h) {
+			compute(pe * (2-pe), i, h/2);
+			compute(pe * pe, i+h, h/2);
+		} else {
+			prob[i] = pe;
+		}
+	}
+	long double prob[1<<MAX_M];
+	int index[1<<MAX_M];
+public:
+	void operator()(uint8_t *frozen_bits, int level, int K, long double erasure_probability = std::exp(-1.L))
+	{
+		assert(level <= MAX_M);
+		int length = 1 << level;
+		compute(erasure_probability, 0, length / 2);
+		for (int i = 0; i < length; ++i)
+			index[i] = i;
+		std::nth_element(index, index+K, index+length, [this](int a, int b){ return prob[a] < prob[b]; });
+		for (int i = 0; i < K; ++i)
+			frozen_bits[index[i]] = 0;
+		for (int i = K; i < length; ++i)
+			frozen_bits[index[i]] = 1;
+	}
+};
+
+}
+
diff --git a/polar_helper.hh b/polar_helper.hh
new file mode 100644
index 0000000..b31b5d4
--- /dev/null
+++ b/polar_helper.hh
@@ -0,0 +1,206 @@
+/*
+SIMD wrapper used by polar encoder and decoder
+
+Copyright 2020 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include "simd.hh"
+
+namespace CODE {
+
+template <typename TYPE>
+struct PolarHelper
+{
+	typedef TYPE PATH;
+	static TYPE one()
+	{
+		return 1;
+	}
+	static TYPE zero()
+	{
+		return 0;
+	}
+	static TYPE signum(TYPE v)
+	{
+		return (v > 0) - (v < 0);
+	}
+	template <typename IN>
+	static TYPE quant(IN in)
+	{
+		return in;
+	}
+	static TYPE qabs(TYPE a)
+	{
+		return std::abs(a);
+	}
+	static TYPE qmin(TYPE a, TYPE b)
+	{
+		return std::min(a, b);
+	}
+	static TYPE qadd(TYPE a, TYPE b)
+	{
+		return a + b;
+	}
+	static TYPE qmul(TYPE a, TYPE b)
+	{
+		return a * b;
+	}
+	static TYPE prod(TYPE a, TYPE b)
+	{
+		return signum(a) * signum(b) * qmin(qabs(a), qabs(b));
+	}
+	static TYPE madd(TYPE a, TYPE b, TYPE c)
+	{
+		return a * b + c;
+	}
+};
+
+template <typename VALUE, int WIDTH>
+struct PolarHelper<SIMD<VALUE, WIDTH>>
+{
+	typedef SIMD<VALUE, WIDTH> TYPE;
+	typedef VALUE PATH;
+	typedef SIMD<typename TYPE::uint_type, WIDTH> MAP;
+	static TYPE one()
+	{
+		return vdup<TYPE>(1);
+	}
+	static TYPE zero()
+	{
+		return vzero<TYPE>();
+	}
+	static TYPE signum(TYPE a)
+	{
+		return vsignum(a);
+	}
+	static TYPE qabs(TYPE a)
+	{
+		return vabs(a);
+	}
+	static TYPE qmin(TYPE a, TYPE b)
+	{
+		return vmin(a, b);
+	}
+	static TYPE qadd(TYPE a, TYPE b)
+	{
+		return vadd(a, b);
+	}
+	static TYPE qmul(TYPE a, TYPE b)
+	{
+		return vmul(a, b);
+	}
+	static TYPE prod(TYPE a, TYPE b)
+	{
+		return vmul(vmul(vsignum(a), vsignum(b)), vmin(vabs(a), vabs(b)));
+	}
+	static TYPE madd(TYPE a, TYPE b, TYPE c)
+	{
+		return vadd(vmul(a, b), c);
+	}
+};
+
+template <int WIDTH>
+struct PolarHelper<SIMD<int8_t, WIDTH>>
+{
+	typedef SIMD<int8_t, WIDTH> TYPE;
+	typedef int PATH;
+	typedef SIMD<uint8_t, WIDTH> MAP;
+	static TYPE one()
+	{
+		return vdup<TYPE>(1);
+	}
+	static TYPE zero()
+	{
+		return vzero<TYPE>();
+	}
+	static TYPE signum(TYPE a)
+	{
+		return vsignum(a);
+	}
+	static TYPE qabs(TYPE a)
+	{
+		return vqabs(a);
+	}
+	static TYPE qadd(TYPE a, TYPE b)
+	{
+		return vqadd(a, b);
+	}
+	static TYPE qmul(TYPE a, TYPE b)
+	{
+#ifdef __ARM_NEON__
+		return vmul(a, b);
+#else
+		return vsign(a, b);
+#endif
+	}
+	static TYPE prod(TYPE a, TYPE b)
+	{
+#ifdef __ARM_NEON__
+		return vmul(vmul(vsignum(a), vsignum(b)), vmin(vqabs(a), vqabs(b)));
+#else
+		return vsign(vmin(vqabs(a), vqabs(b)), vsign(vsignum(a), b));
+#endif
+	}
+	static TYPE madd(TYPE a, TYPE b, TYPE c)
+	{
+#ifdef __ARM_NEON__
+		return vqadd(vmul(a, vmax(b, vdup<TYPE>(-127))), c);
+#else
+		return vqadd(vsign(vmax(b, vdup<TYPE>(-127)), a), c);
+#endif
+	}
+};
+
+template <>
+struct PolarHelper<int8_t>
+{
+	typedef int PATH;
+	static int8_t one()
+	{
+		return 1;
+	}
+	static int8_t zero()
+	{
+		return 0;
+	}
+	static int8_t signum(int8_t v)
+	{
+		return (v > 0) - (v < 0);
+	}
+	template <typename IN>
+	static int8_t quant(IN in)
+	{
+		return std::min<IN>(std::max<IN>(std::nearbyint(in), -128), 127);
+	}
+	static int8_t qabs(int8_t a)
+	{
+		return std::abs(std::max<int8_t>(a, -127));
+	}
+	static int8_t qmin(int8_t a, int8_t b)
+	{
+		return std::min(a, b);
+	}
+	static int8_t qadd(int8_t a, int8_t b)
+	{
+		return std::min<int16_t>(std::max<int16_t>(int16_t(a) + int16_t(b), -128), 127);
+	}
+	static int8_t qmul(int8_t a, int8_t b)
+	{
+		// return std::min<int16_t>(std::max<int16_t>(int16_t(a) * int16_t(b), -128), 127);
+		// only used for hard decision values anyway
+		return a * b;
+	}
+	static int8_t prod(int8_t a, int8_t b)
+	{
+		return signum(a) * signum(b) * qmin(qabs(a), qabs(b));
+	}
+	static int8_t madd(int8_t a, int8_t b, int8_t c)
+	{
+		return std::min<int16_t>(std::max<int16_t>(int16_t(a) * int16_t(b) + int16_t(c), -128), 127);
+	}
+};
+
+}
+
diff --git a/polar_list_decoder.hh b/polar_list_decoder.hh
new file mode 100644
index 0000000..7810e28
--- /dev/null
+++ b/polar_list_decoder.hh
@@ -0,0 +1,146 @@
+/*
+Successive cancellation list decoding of polar codes
+
+Copyright 2020 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include <algorithm>
+#include "polar_helper.hh"
+
+namespace CODE {
+
+template <typename TYPE, int M>
+struct PolarListTree
+{
+	typedef PolarHelper<TYPE> PH;
+	typedef typename PH::PATH PATH;
+	typedef typename PH::MAP MAP;
+	static const int N = 1 << M;
+	static MAP decode(PATH *metric, TYPE *message, MAP *maps, int *count, TYPE *hard, TYPE *soft, const uint8_t *frozen)
+	{
+		for (int i = 0; i < N/2; ++i)
+			soft[i+N/2] = PH::prod(soft[i+N], soft[i+N/2+N]);
+		MAP lmap = PolarListTree<TYPE, M-1>::decode(metric, message, maps, count, hard, soft, frozen);
+		for (int i = 0; i < N/2; ++i)
+			soft[i+N/2] = PH::madd(hard[i], vshuf(soft[i+N], lmap), vshuf(soft[i+N/2+N], lmap));
+		MAP rmap = PolarListTree<TYPE, M-1>::decode(metric, message, maps, count, hard+N/2, soft, frozen+N/2);
+		for (int i = 0; i < N/2; ++i)
+			hard[i] = PH::qmul(vshuf(hard[i], rmap), hard[i+N/2]);
+		return vshuf(lmap, rmap);
+	}
+};
+
+template <typename TYPE>
+struct PolarListTree<TYPE, 0>
+{
+	typedef PolarHelper<TYPE> PH;
+	typedef typename PH::PATH PATH;
+	typedef typename PH::MAP MAP;
+	static MAP decode(PATH *metric, TYPE *message, MAP *maps, int *count, TYPE *hard, TYPE *soft, const uint8_t *frozen)
+	{
+		MAP map;
+		TYPE hrd, sft = soft[1];
+		if (*frozen) {
+			for (int k = 0; k < TYPE::SIZE; ++k)
+				if (sft.v[k] < 0)
+					metric[k] -= sft.v[k];
+			hrd = PH::one();
+			for (int k = 0; k < TYPE::SIZE; ++k)
+				map.v[k] = k;
+		} else {
+			PATH fork[2*TYPE::SIZE];
+			for (int k = 0; k < TYPE::SIZE; ++k)
+				fork[k] = fork[k+TYPE::SIZE] = metric[k];
+			for (int k = 0; k < TYPE::SIZE; ++k)
+				if (sft.v[k] < 0)
+					fork[k] -= sft.v[k];
+				else
+					fork[k+TYPE::SIZE] += sft.v[k];
+			int perm[2*TYPE::SIZE];
+			for (int k = 0; k < 2*TYPE::SIZE; ++k)
+				perm[k] = k;
+			std::nth_element(perm, perm+TYPE::SIZE, perm+2*TYPE::SIZE, [fork](int a, int b){ return fork[a] < fork[b]; });
+			for (int k = 0; k < TYPE::SIZE; ++k)
+				metric[k] = fork[perm[k]];
+			for (int k = 0; k < TYPE::SIZE; ++k)
+				map.v[k] = perm[k] % TYPE::SIZE;
+			for (int k = 0; k < TYPE::SIZE; ++k)
+				hrd.v[k] = perm[k] < TYPE::SIZE ? 1 : -1;
+			message[*count] = hrd;
+			maps[*count] = map;
+			++*count;
+		}
+		*hard = hrd;
+		return map;
+	}
+};
+
+template <typename TYPE, int MAX_M>
+class PolarListDecoder
+{
+	typedef PolarHelper<TYPE> PH;
+	typedef typename TYPE::value_type VALUE;
+	typedef typename PH::PATH PATH;
+	typedef typename PH::MAP MAP;
+	static const int MAX_N = 1 << MAX_M;
+	TYPE soft[2*MAX_N];
+	TYPE hard[MAX_N];
+	MAP maps[MAX_N];
+public:
+	void operator()(PATH *metric, TYPE *message, const VALUE *codeword, const uint8_t *frozen, int level)
+	{
+		assert(level <= MAX_M);
+		int count = 0;
+		metric[0] = 0;
+		for (int k = 1; k < TYPE::SIZE; ++k)
+			metric[k] = 1000;
+		int length = 1 << level;
+		for (int i = 0; i < length; ++i)
+			soft[length+i] = vdup<TYPE>(codeword[i]);
+
+		switch (level) {
+		case 0: PolarListTree<TYPE, 0>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 1: PolarListTree<TYPE, 1>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 2: PolarListTree<TYPE, 2>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 3: PolarListTree<TYPE, 3>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 4: PolarListTree<TYPE, 4>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 5: PolarListTree<TYPE, 5>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 6: PolarListTree<TYPE, 6>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 7: PolarListTree<TYPE, 7>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 8: PolarListTree<TYPE, 8>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 9: PolarListTree<TYPE, 9>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 10: PolarListTree<TYPE, 10>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 11: PolarListTree<TYPE, 11>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 12: PolarListTree<TYPE, 12>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 13: PolarListTree<TYPE, 13>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 14: PolarListTree<TYPE, 14>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 15: PolarListTree<TYPE, 15>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 16: PolarListTree<TYPE, 16>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 17: PolarListTree<TYPE, 17>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 18: PolarListTree<TYPE, 18>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 19: PolarListTree<TYPE, 19>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 20: PolarListTree<TYPE, 20>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 21: PolarListTree<TYPE, 21>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 22: PolarListTree<TYPE, 22>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 23: PolarListTree<TYPE, 23>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 24: PolarListTree<TYPE, 24>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 25: PolarListTree<TYPE, 25>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 26: PolarListTree<TYPE, 26>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 27: PolarListTree<TYPE, 27>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 28: PolarListTree<TYPE, 28>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		case 29: PolarListTree<TYPE, 29>::decode(metric, message, maps, &count, hard, soft, frozen); break;
+		default: assert(false);
+		}
+
+		MAP acc = maps[count-1];
+		for (int i = count-2; i >= 0; --i) {
+			message[i] = vshuf(message[i], acc);
+			acc = vshuf(maps[i], acc);
+		}
+	}
+};
+
+}
+
diff --git a/tests/polar_list_regression_test.cc b/tests/polar_list_regression_test.cc
new file mode 100644
index 0000000..dd40282
--- /dev/null
+++ b/tests/polar_list_regression_test.cc
@@ -0,0 +1,202 @@
+/*
+Regression Test for the Polar Encoder and List Decoders
+
+Copyright 2020 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <limits>
+#include <random>
+#include <chrono>
+#include <cassert>
+#include <iomanip>
+#include <iostream>
+#include <algorithm>
+#include <functional>
+#include "polar_helper.hh"
+#include "polar_list_decoder.hh"
+#include "polar_encoder.hh"
+#include "polar_freezer.hh"
+
+int main()
+{
+	const int M = 16;
+	const int N = 1 << M;
+	const bool systematic = true;
+#if 0
+	typedef int8_t code_type;
+	double SCALE = 2;
+#else
+	typedef float code_type;
+	double SCALE = 1;
+#endif
+
+#ifdef __AVX2__
+	const int SIZEOF_SIMD = 32;
+#else
+	const int SIZEOF_SIMD = 16;
+#endif
+	const int SIMD_WIDTH = SIZEOF_SIMD / sizeof(code_type);
+	typedef SIMD<code_type, SIMD_WIDTH> simd_type;
+
+	std::random_device rd;
+	typedef std::default_random_engine generator;
+	typedef std::uniform_int_distribution<int> distribution;
+	auto data = std::bind(distribution(0, 1), generator(rd()));
+	auto frozen = new uint8_t[N];
+	auto codeword = new code_type[N];
+	auto temp = new simd_type[N];
+
+	long double erasure_probability = 1. / 3.;
+	int K = (1 - erasure_probability) * N;
+	double design_SNR = 10 * std::log10(-std::log(erasure_probability));
+	std::cerr << "design SNR: " << design_SNR << std::endl;
+	if (0) {
+		CODE::PolarFreezer freeze;
+		long double freezing_threshold = 0 ? 0.5 : std::numeric_limits<float>::epsilon();
+		K = freeze(frozen, M, erasure_probability, freezing_threshold);
+	} else {
+		auto freeze = new CODE::PolarCodeConst0<M>;
+		std::cerr << "sizeof(PolarCodeConst0<M>) = " << sizeof(CODE::PolarCodeConst0<M>) << std::endl;
+		double better_SNR = design_SNR + 1.59175;
+		std::cerr << "better SNR: " << better_SNR << std::endl;
+		long double probability = std::exp(-pow(10.0, better_SNR / 10));
+		(*freeze)(frozen, M, K, probability);
+		delete freeze;
+	}
+	std::cerr << "Polar(" << N << ", " << K << ")" << std::endl;
+	auto message = new code_type[K];
+	auto decoded = new simd_type[K];
+	CODE::PolarHelper<simd_type>::PATH metric[SIMD_WIDTH];
+	std::cerr << "sizeof(PolarListDecoder<simd_type, M>) = " << sizeof(CODE::PolarListDecoder<simd_type, M>) << std::endl;
+	auto decode = new CODE::PolarListDecoder<simd_type, M>;
+
+	auto orig = new code_type[N];
+	auto noisy = new code_type[N];
+	auto symb = new double[N];
+	double low_SNR = std::floor(design_SNR-3);
+	double high_SNR = std::ceil(design_SNR+5);
+	double min_SNR = high_SNR, max_mbs = 0;
+	int count = 0;
+	std::cerr << "SNR BER Mbit/s Eb/N0" << std::endl;
+	for (double SNR = low_SNR; count <= 3 && SNR <= high_SNR; SNR += 0.1, ++count) {
+		//double mean_signal = 0;
+		double sigma_signal = 1;
+		double mean_noise = 0;
+		double sigma_noise = std::sqrt(sigma_signal * sigma_signal / (2 * std::pow(10, SNR / 10)));
+
+		typedef std::normal_distribution<double> normal;
+		auto awgn = std::bind(normal(mean_noise, sigma_noise), generator(rd()));
+
+		int64_t awgn_errors = 0;
+		int64_t quantization_erasures = 0;
+		int64_t uncorrected_errors = 0;
+		int64_t ambiguity_erasures = 0;
+		double avg_mbs = 0;
+		int64_t loops = 0;
+		while (uncorrected_errors < 1000 && ++loops < 100) {
+			for (int i = 0; i < K; ++i)
+				message[i] = 1 - 2 * data();
+
+			if (systematic) {
+				CODE::PolarSysEnc<code_type> sysenc;
+				sysenc(codeword, message, frozen, M);
+				for (int i = 0, j = 0; i < N; ++i)
+					if (!frozen[i])
+						assert(codeword[i] == message[j++]);
+			} else {
+				CODE::PolarEncoder<code_type> encode;
+				encode(codeword, message, frozen, M);
+			}
+
+			for (int i = 0; i < N; ++i)
+				orig[i] = codeword[i];
+
+			for (int i = 0; i < N; ++i)
+				symb[i] = codeword[i];
+
+			for (int i = 0; i < N; ++i)
+				symb[i] += awgn();
+
+			// $LLR=log(\frac{p(x=+1|y)}{p(x=-1|y)})$
+			// $p(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$
+			double DIST = 2; // BPSK
+			double fact = SCALE * DIST / (sigma_noise * sigma_noise);
+			for (int i = 0; i < N; ++i)
+				codeword[i] = CODE::PolarHelper<code_type>::quant(fact * symb[i]);
+
+			for (int i = 0; i < N; ++i)
+				noisy[i] = codeword[i];
+
+			auto start = std::chrono::system_clock::now();
+			(*decode)(metric, decoded, codeword, frozen, M);
+			auto end = std::chrono::system_clock::now();
+			auto usec = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
+			double mbs = (double)K / usec.count();
+			avg_mbs += mbs;
+
+			if (systematic) {
+				CODE::PolarEncoder<simd_type> encode;
+				encode(temp, decoded, frozen, M);
+				for (int i = 0, j = 0; i < N; ++i)
+					if (!frozen[i])
+						decoded[j++] = temp[i];
+			}
+
+			int best = 0;
+			if (1) {
+				for (int k = 0; k < SIMD_WIDTH; ++k)
+					if (metric[k] < metric[best])
+						best = k;
+			} else {
+				int errs[SIMD_WIDTH] = { 0 };
+				for (int i = 0; i < K; ++i)
+					for (int k = 0; k < SIMD_WIDTH; ++k)
+						errs[k] += message[i] != decoded[i].v[k];
+				for (int k = 0; k < SIMD_WIDTH; ++k)
+					if (errs[k] < errs[best])
+						best = k;
+			}
+
+			for (int i = 0; i < N; ++i)
+				awgn_errors += noisy[i] * (orig[i] < 0);
+			for (int i = 0; i < N; ++i)
+				quantization_erasures += !noisy[i];
+			for (int i = 0; i < K; ++i)
+				uncorrected_errors += decoded[i].v[best] * message[i] <= 0;
+			for (int i = 0; i < K; ++i)
+				ambiguity_erasures += !decoded[i].v[best];
+		}
+
+		avg_mbs /= loops;
+
+		max_mbs = std::max(max_mbs, avg_mbs);
+		double bit_error_rate = (double)uncorrected_errors / (double)(K * loops);
+		if (!uncorrected_errors)
+			min_SNR = std::min(min_SNR, SNR);
+		else
+			count = 0;
+
+		int MOD_BITS = 1; // BPSK
+		double code_rate = (double)K / (double)N;
+		double spectral_efficiency = code_rate * MOD_BITS;
+		double EbN0 = 10 * std::log10(sigma_signal * sigma_signal / (spectral_efficiency * 2 * sigma_noise * sigma_noise));
+
+		if (0) {
+			std::cerr << SNR << " Es/N0 => AWGN with standard deviation of " << sigma_noise << " and mean " << mean_noise << std::endl;
+			std::cerr << EbN0 << " Eb/N0, using spectral efficiency of " << spectral_efficiency << " from " << code_rate << " code rate and " << MOD_BITS << " bits per symbol." << std::endl;
+			std::cerr << awgn_errors << " errors caused by AWGN." << std::endl;
+			std::cerr << quantization_erasures << " erasures caused by quantization." << std::endl;
+			std::cerr << uncorrected_errors << " errors uncorrected." << std::endl;
+			std::cerr << ambiguity_erasures << " ambiguity erasures." << std::endl;
+			std::cerr << bit_error_rate << " bit error rate." << std::endl;
+			std::cerr << avg_mbs << " megabit per second." << std::endl;
+		} else {
+			std::cout << SNR << " " << bit_error_rate << " " << avg_mbs << " " << EbN0 << std::endl;
+		}
+	}
+	std::cerr << "QEF at: " << min_SNR << " SNR, speed: " << max_mbs << " Mb/s." << std::endl;
+	double QEF_SNR = design_SNR + 0.2;
+	assert(min_SNR < QEF_SNR);
+	std::cerr << "Polar list regression test passed!" << std::endl;
+	return 0;
+}
diff --git a/tests/polar_regression_test.cc b/tests/polar_regression_test.cc
new file mode 100644
index 0000000..07974a4
--- /dev/null
+++ b/tests/polar_regression_test.cc
@@ -0,0 +1,178 @@
+/*
+Regression Test for the Polar Encoder and Decoders
+
+Copyright 2020 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <limits>
+#include <random>
+#include <chrono>
+#include <cassert>
+#include <iomanip>
+#include <iostream>
+#include <algorithm>
+#include <functional>
+#include "polar_helper.hh"
+#include "polar_decoder.hh"
+#include "polar_encoder.hh"
+#include "polar_freezer.hh"
+
+int main()
+{
+	const int M = 20;
+	const int N = 1 << M;
+	const bool systematic = true;
+#if 1
+	typedef int8_t code_type;
+	double SCALE = 2;
+#else
+	typedef float code_type;
+	double SCALE = 1;
+#endif
+
+	std::random_device rd;
+	typedef std::default_random_engine generator;
+	typedef std::uniform_int_distribution<int> distribution;
+	auto data = std::bind(distribution(0, 1), generator(rd()));
+	auto frozen = new uint8_t[N];
+	auto codeword = new code_type[N];
+	auto temp = new code_type[N];
+
+	long double erasure_probability = 1. / 3.;
+	int K = (1 - erasure_probability) * N;
+	double design_SNR = 10 * std::log10(-std::log(erasure_probability));
+	std::cerr << "design SNR: " << design_SNR << std::endl;
+	if (0) {
+		CODE::PolarFreezer freeze;
+		long double freezing_threshold = 0 ? 0.5 : std::numeric_limits<float>::epsilon();
+		K = freeze(frozen, M, erasure_probability, freezing_threshold);
+	} else {
+		auto freeze = new CODE::PolarCodeConst0<M>;
+		std::cerr << "sizeof(PolarCodeConst0<M>) = " << sizeof(CODE::PolarCodeConst0<M>) << std::endl;
+		double better_SNR = design_SNR + 0.5;//1.59175;
+		std::cerr << "better SNR: " << better_SNR << std::endl;
+		long double probability = std::exp(-pow(10.0, better_SNR / 10));
+		(*freeze)(frozen, M, K, probability);
+		delete freeze;
+	}
+	std::cerr << "Polar(" << N << ", " << K << ")" << std::endl;
+	auto message = new code_type[K];
+	auto decoded = new code_type[K];
+	std::cerr << "sizeof(PolarDecoder<code_type, M>) = " << sizeof(CODE::PolarDecoder<code_type, M>) << std::endl;
+	auto decode = new CODE::PolarDecoder<code_type, M>;
+
+	auto orig = new code_type[N];
+	auto noisy = new code_type[N];
+	auto symb = new double[N];
+	double low_SNR = std::floor(design_SNR-3);
+	double high_SNR = std::ceil(design_SNR+5);
+	double min_SNR = high_SNR, max_mbs = 0;
+	int count = 0;
+	std::cerr << "SNR BER Mbit/s Eb/N0" << std::endl;
+	for (double SNR = low_SNR; count <= 3 && SNR <= high_SNR; SNR += 0.1, ++count) {
+		//double mean_signal = 0;
+		double sigma_signal = 1;
+		double mean_noise = 0;
+		double sigma_noise = std::sqrt(sigma_signal * sigma_signal / (2 * std::pow(10, SNR / 10)));
+
+		typedef std::normal_distribution<double> normal;
+		auto awgn = std::bind(normal(mean_noise, sigma_noise), generator(rd()));
+
+		int64_t awgn_errors = 0;
+		int64_t quantization_erasures = 0;
+		int64_t uncorrected_errors = 0;
+		int64_t ambiguity_erasures = 0;
+		double avg_mbs = 0;
+		int64_t loops = 0;
+		while (uncorrected_errors < 1000 && ++loops < 100) {
+			for (int i = 0; i < K; ++i)
+				message[i] = 1 - 2 * data();
+
+			if (systematic) {
+				CODE::PolarSysEnc<code_type> sysenc;
+				sysenc(codeword, message, frozen, M);
+				for (int i = 0, j = 0; i < N; ++i)
+					if (!frozen[i])
+						assert(codeword[i] == message[j++]);
+			} else {
+				CODE::PolarEncoder<code_type> encode;
+				encode(codeword, message, frozen, M);
+			}
+
+			for (int i = 0; i < N; ++i)
+				orig[i] = codeword[i];
+
+			for (int i = 0; i < N; ++i)
+				symb[i] = codeword[i];
+
+			for (int i = 0; i < N; ++i)
+				symb[i] += awgn();
+
+			// $LLR=log(\frac{p(x=+1|y)}{p(x=-1|y)})$
+			// $p(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$
+			double DIST = 2; // BPSK
+			double fact = SCALE * DIST / (sigma_noise * sigma_noise);
+			for (int i = 0; i < N; ++i)
+				codeword[i] = CODE::PolarHelper<code_type>::quant(fact * symb[i]);
+
+			for (int i = 0; i < N; ++i)
+				noisy[i] = codeword[i];
+
+			auto start = std::chrono::system_clock::now();
+			(*decode)(decoded, codeword, frozen, M);
+			auto end = std::chrono::system_clock::now();
+			auto usec = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
+			double mbs = (double)K / usec.count();
+			avg_mbs += mbs;
+
+			if (systematic) {
+				CODE::PolarEncoder<code_type> encode;
+				encode(temp, decoded, frozen, M);
+				for (int i = 0, j = 0; i < N; ++i)
+					if (!frozen[i])
+						decoded[j++] = temp[i];
+			}
+
+			for (int i = 0; i < N; ++i)
+				awgn_errors += noisy[i] * (orig[i] < 0);
+			for (int i = 0; i < N; ++i)
+				quantization_erasures += !noisy[i];
+			for (int i = 0; i < K; ++i)
+				uncorrected_errors += decoded[i] * message[i] <= 0;
+			for (int i = 0; i < K; ++i)
+				ambiguity_erasures += !decoded[i];
+		}
+
+		avg_mbs /= loops;
+
+		max_mbs = std::max(max_mbs, avg_mbs);
+		double bit_error_rate = (double)uncorrected_errors / (double)(K * loops);
+		if (!uncorrected_errors)
+			min_SNR = std::min(min_SNR, SNR);
+		else
+			count = 0;
+
+		int MOD_BITS = 1; // BPSK
+		double code_rate = (double)K / (double)N;
+		double spectral_efficiency = code_rate * MOD_BITS;
+		double EbN0 = 10 * std::log10(sigma_signal * sigma_signal / (spectral_efficiency * 2 * sigma_noise * sigma_noise));
+
+		if (0) {
+			std::cerr << SNR << " Es/N0 => AWGN with standard deviation of " << sigma_noise << " and mean " << mean_noise << std::endl;
+			std::cerr << EbN0 << " Eb/N0, using spectral efficiency of " << spectral_efficiency << " from " << code_rate << " code rate and " << MOD_BITS << " bits per symbol." << std::endl;
+			std::cerr << awgn_errors << " errors caused by AWGN." << std::endl;
+			std::cerr << quantization_erasures << " erasures caused by quantization." << std::endl;
+			std::cerr << uncorrected_errors << " errors uncorrected." << std::endl;
+			std::cerr << ambiguity_erasures << " ambiguity erasures." << std::endl;
+			std::cerr << bit_error_rate << " bit error rate." << std::endl;
+			std::cerr << avg_mbs << " megabit per second." << std::endl;
+		} else {
+			std::cout << SNR << " " << bit_error_rate << " " << avg_mbs << " " << EbN0 << std::endl;
+		}
+	}
+	std::cerr << "QEF at: " << min_SNR << " SNR, speed: " << max_mbs << " Mb/s." << std::endl;
+	double QEF_SNR = design_SNR + 0.2;
+	assert(min_SNR < QEF_SNR);
+	std::cerr << "Polar regression test passed!" << std::endl;
+	return 0;
+}