text:
#ifndef		HAS_QEF_SIMD_H_BEEN_INCLUDED
#define		HAS_QEF_SIMD_H_BEEN_INCLUDED

//
// Quadric Error Function / Singluar Value Decomposition SSE2 implementation
// Public domain
//
// Input is a set of positions / vertices of a surface and the surface normals
// at these positions (i.e. Hermite data). A "best fit" position is calculated
// from these positions along with an error value.
//
// Exmaple use cases:
//   - placing a vertex inside a voxel/node for algorithms like Dual Contouring
//     (the positions and normals corresponding to edge/isosurface intersections)
//   - placing the vertex resulting from collapsing an edge when simplifying a mesh
//     (the positions and normals of the collapsed edge should be supplied)
//
// Usage:
//
// SSE registers:
//
//	__m128 position[2] = { ... };
//	__m128 normal[2] = { ... };
//	__m128 solvedPos;
//	float error = qef_solve_from_points(position, normal, 2, &solvedPos);
//
// 4D vectors:
//
//	__declspec(align(16)) float positions[2 * 4] = { ... };
//	__declspec(align(16)) float normals[2 * 4] = { ... };
//	__declspec(align(16)) float solved[4];
//	float error = qef_solve_from_points_4d(positions, normals, 2, solvedPos);
//
// 3D vectors (or struct with 3 float members, e.g. glm::vec3):
//
//	glm::vec3 positions[2] = { ... };
//	glm::vec3 normals[2] = { ... };
//	glm::vec3 solvedPos;
//	float error = qef_solve_from_points_3d(&positions[0].x, &normals[0].x, 2, &solvedPos.x);
//

#include	<xmmintrin.h>
#include	<immintrin.h>

const int QEF_MAX_INPUT_COUNT = 12;

// Ideally the data would already be in SSE registers & returned in a SEE register
float qef_solve_from_points(
	const __m128* positions,
	const __m128* normals,
	const int count,
	__m128* solved_position);

// Expects 4d vectors contiguous in memory for the positions/normals
// Addresses pointed to by positions, normals and solved_position MUST be 16 byte aligned
// Writes result to 4d vector.
float qef_solve_from_points_4d(
	const float* positions,
	const float* normals,
	const int count,
	float* solved_position);

float qef_solve_from_points_4d_interleaved(
	const float* data,
	const size_t stride,
	const int count,
	float* solved_position);

// Expects 3d vectors contiguous in memory for the positions/normals
// No alignment requirements.
// Writes result to 3d vector.
float qef_solve_from_points_3d(
	const float* positions,
	const float* normals,
	const int count,
	float* solved_position);

#ifdef QEF_INCLUDE_IMPL

union Mat4x4
{
	float	m[4][4];
	__m128	row[4];
};

#define SVD_NUM_SWEEPS 5
const float PSUEDO_INVERSE_THRESHOLD = 0.001f;

// ----------------------------------------------------------------------------

static inline __m128 vec4_abs(const __m128& x)
{
	static const __m128 mask = _mm_set1_ps(-0.f);
	return _mm_andnot_ps(mask, x);
}

// ----------------------------------------------------------------------------

static inline float vec4_dot(const __m128& a, const __m128& b)
{
	// apparently _mm_dp_ps isn't well implemented in hardware so this "basic" version is faster
	__m128 mul = _mm_mul_ps(a, b);
	__m128 s0  = _mm_shuffle_ps(mul, mul, _MM_SHUFFLE(2, 3, 0, 1));
	__m128 add = _mm_add_ps(mul, s0);
	__m128 s1  = _mm_shuffle_ps(add, add, _MM_SHUFFLE(0, 1, 2, 3));
	__m128 res = _mm_add_ps(add, s1);
	return _mm_cvtss_f32(res);
}

// ----------------------------------------------------------------------------

static inline __m128 vec4_mul_m4x4(const __m128& a, const Mat4x4& B)
{
	auto bob = _mm_shuffle_ps(a, a, 0x00);

	__m128 result = _mm_mul_ps(_mm_shuffle_ps(a, a, 0x00), B.row[0]);
	result = _mm_add_ps(result, _mm_mul_ps(_mm_shuffle_ps(a, a, 0x55), B.row[1]));
	result = _mm_add_ps(result, _mm_mul_ps(_mm_shuffle_ps(a, a, 0xaa), B.row[2]));
	result = _mm_add_ps(result, _mm_mul_ps(_mm_shuffle_ps(a, a, 0xff), B.row[3]));
	return result;
}

// ----------------------------------------------------------------------------

static inline __m256 avx_vec4_mul_m4x4(const __m256& a, const Mat4x4& B)
{
	__m256 result;
	result = _mm256_mul_ps(_mm256_shuffle_ps(a, a, 0x00), _mm256_broadcast_ps(&B.row[0]));
	result = _mm256_add_ps(result, _mm256_mul_ps(_mm256_shuffle_ps(a, a, 0x55), _mm256_broadcast_ps(&B.row[1])));
	result = _mm256_add_ps(result, _mm256_mul_ps(_mm256_shuffle_ps(a, a, 0xaa), _mm256_broadcast_ps(&B.row[2])));
	result = _mm256_add_ps(result, _mm256_mul_ps(_mm256_shuffle_ps(a, a, 0xff), _mm256_broadcast_ps(&B.row[3])));

	return result;
}

// ----------------------------------------------------------------------------

static void m4x4_mul_m4x4(Mat4x4& out, const Mat4x4& A, const Mat4x4& B)
{
	out.row[0] = vec4_mul_m4x4(A.row[0], B);
	out.row[1] = vec4_mul_m4x4(A.row[1], B);
	out.row[2] = vec4_mul_m4x4(A.row[2], B);
	out.row[3] = vec4_mul_m4x4(A.row[3], B);
}

// ----------------------------------------------------------------------------

static void givens_coeffs_sym(__m128& c_result, __m128& s_result, const Mat4x4& vtav, const int a, const int b)
{
	__m128 simd_pp = _mm_set_ps(
		0.f,
		vtav.row[a].m128_f32[a],
		vtav.row[a].m128_f32[a],
		vtav.row[a].m128_f32[a]);

	__m128 simd_pq = _mm_set_ps(
		0.f,
		vtav.row[a].m128_f32[b],
		vtav.row[a].m128_f32[b],
		vtav.row[a].m128_f32[b]);

	__m128 simd_qq = _mm_set_ps(
		0.f,
		vtav.row[b].m128_f32[b],
		vtav.row[b].m128_f32[b],
		vtav.row[b].m128_f32[b]);

	static const __m128 zeros = _mm_set1_ps(0.f);
	static const __m128 ones  = _mm_set1_ps(1.f);
	static const __m128 twos  = _mm_set1_ps(2.f);

	// tau = (a_qq - a_pp) / (2.f * a_pq);
	__m128 pq2 = _mm_mul_ps(simd_pq, twos);
	__m128 qq_sub_pp = _mm_sub_ps(simd_qq, simd_pp);
	__m128 tau = _mm_div_ps(qq_sub_pp, pq2);

	// stt = sqrt(1.f + tau * tau);
	__m128 tau_sq = _mm_mul_ps(tau, tau);
	__m128 tau_sq_1 = _mm_add_ps(tau_sq, ones);
	__m128 stt = _mm_sqrt_ps(tau_sq_1);

	// tan = 1.f / ((tau >= 0.f) ? (tau + stt) : (tau - stt));
	__m128 tan_gt = _mm_add_ps(tau, stt);
	__m128 tan_lt = _mm_sub_ps(tau, stt);
	__m128 tan_cmp = _mm_cmpge_ps(tau, zeros);
	__m128 tan_cmp_gt = _mm_and_ps(tan_cmp, tan_gt);
	__m128 tan_cmp_lt = _mm_andnot_ps(tan_cmp, tan_lt);
	__m128 tan_inv = _mm_or_ps(tan_cmp_gt, tan_cmp_lt);
	__m128 tan = _mm_div_ps(ones, tan_inv);

	// c = rsqrt(1.f + tan * tan);
	__m128 tan_sq = _mm_mul_ps(tan, tan);
	__m128 tan_sq_1 = _mm_add_ps(ones, tan_sq);
	__m128 c = _mm_rsqrt_ps(tan_sq_1);

	// s = tan * c;
	__m128 s = _mm_mul_ps(tan, c);

	// if pq == 0.0: c = 1.f, s = 0.f
	__m128 pq_cmp = _mm_cmpeq_ps(simd_pq, zeros);

	__m128 c_true = _mm_and_ps(pq_cmp, ones);
	__m128 c_false = _mm_andnot_ps(pq_cmp, c);
	c_result = _mm_or_ps(c_true, c_false);

	__m128 s_true = _mm_and_ps(pq_cmp, zeros);
	__m128 s_false = _mm_andnot_ps(pq_cmp, s);
	s_result = _mm_or_ps(s_true, s_false);
}

// ----------------------------------------------------------------------------

static void rotateq_xy(Mat4x4& vtav, const __m128& c, const __m128& s, const int a, const int b)
{
	__m128 u = _mm_set_ps(
		0.f,
		vtav.row[a].m128_f32[a],
		vtav.row[a].m128_f32[a],
		vtav.row[a].m128_f32[a]);

	__m128 v = _mm_set_ps(
		0.f,
		vtav.row[b].m128_f32[b],
		vtav.row[b].m128_f32[b],
		vtav.row[b].m128_f32[b]);

	__m128 A = _mm_set_ps(
		0.f,
		vtav.row[a].m128_f32[b],
		vtav.row[a].m128_f32[b],
		vtav.row[a].m128_f32[b]);

	static const __m128 twos = _mm_set1_ps(2.f);

	__m128 cc = _mm_mul_ps(c, c);
	__m128 ss = _mm_mul_ps(s, s);

	// mx = 2.0 * c * s * A;
	__m128 c2 = _mm_mul_ps(twos, c);
	__m128 c2s = _mm_mul_ps(c2, s);
	__m128 mx = _mm_mul_ps(c2s, A);

	// x = cc * u - mx + ss * v;
	__m128 x0 = _mm_mul_ps(cc, u);
	__m128 x1 = _mm_sub_ps(x0, mx);
	__m128 x2 = _mm_mul_ps(ss, v);
	__m128 x  = _mm_add_ps(x1, x2);

	// y = ss * u + mx + cc * v;
	__m128 y0 = _mm_mul_ps(ss, u);
	__m128 y1 = _mm_add_ps(y0, mx);
	__m128 y2 = _mm_mul_ps(cc, v);
	__m128 y  = _mm_add_ps(y1, y2);

	vtav.row[a].m128_f32[a] = x.m128_f32[0];
	vtav.row[b].m128_f32[b] = y.m128_f32[0];
}

// ----------------------------------------------------------------------------

static void rotate_xy(Mat4x4& vtav, Mat4x4& v, float c, float s, const int& a, const int& b)
{
	 __m128 simd_u = _mm_set_ps(
		vtav.row[0].m128_f32[3-b],
		v.row[2].m128_f32[a],
		v.row[1].m128_f32[a],
		v.row[0].m128_f32[a]);

	__m128 simd_v = _mm_set_ps(
		vtav.row[1-a].m128_f32[2],
		v.row[2].m128_f32[b],
		v.row[1].m128_f32[b],
		v.row[0].m128_f32[b]);

	__m128 simd_c = _mm_load1_ps(&c);
	__m128 simd_s = _mm_load1_ps(&s);

	__m128 x0 = _mm_mul_ps(simd_c, simd_u);
	__m128 x1 = _mm_mul_ps(simd_s, simd_v);
	__m128 x = _mm_sub_ps(x0, x1);

	__m128 y0 = _mm_mul_ps(simd_s, simd_u);
	__m128 y1 = _mm_mul_ps(simd_c, simd_v);
	__m128 y = _mm_add_ps(y0, y1);

	v.row[0].m128_f32[a] = x.m128_f32[0];
	v.row[1].m128_f32[a] = x.m128_f32[1];
	v.row[2].m128_f32[a] = x.m128_f32[2];
	vtav.row[0].m128_f32[3-b] = x.m128_f32[3];

	v.row[0].m128_f32[b] = y.m128_f32[0];
	v.row[1].m128_f32[b] = y.m128_f32[1];
	v.row[2].m128_f32[b] = y.m128_f32[2];
	vtav.row[1-a].m128_f32[2] = y.m128_f32[3];

	vtav.row[a].m128_f32[b] = 0.f;
}

// ----------------------------------------------------------------------------

static __m128 svd_solve_sym(Mat4x4& v, const Mat4x4& a)
{
	Mat4x4 vtav = a;

	for (int i = 0; i < SVD_NUM_SWEEPS; ++i)
	{
		__m128 c, s;

		if (vtav.row[0].m128_f32[1] != 0.f)
		{
			givens_coeffs_sym(c, s, vtav, 0, 1);
			rotateq_xy(vtav, c, s, 0, 1);
			rotate_xy(vtav, v, c.m128_f32[1], s.m128_f32[1], 0, 1);
			vtav.row[0].m128_f32[1] = 0.f;
		}

		if (vtav.row[0].m128_f32[2] != 0.f)
		{
			givens_coeffs_sym(c, s, vtav, 0, 2);
			rotateq_xy(vtav, c, s, 0, 2);
			rotate_xy(vtav, v, c.m128_f32[1], s.m128_f32[1], 0, 2);
			vtav.row[0].m128_f32[2] = 0.f;
		}

		if (vtav.row[1].m128_f32[2] != 0.f)
		{
			givens_coeffs_sym(c, s, vtav, 1, 2);
			rotateq_xy(vtav, c, s, 1, 2);
			rotate_xy(vtav, v, c.m128_f32[2], s.m128_f32[2], 1, 2);
			vtav.row[1].m128_f32[2] = 0.f;
		}
	}

	return _mm_set_ps(
		0.f,
		vtav.row[2].m128_f32[2],
		vtav.row[1].m128_f32[1],
		vtav.row[0].m128_f32[0]);
}

// ----------------------------------------------------------------------------

static inline __m128 svd_invdet(const __m128& x)
{
	static const __m128 ones = _mm_set1_ps(1.f);
	static const __m128 tol = _mm_set1_ps(PSUEDO_INVERSE_THRESHOLD);

	__m128 abs_x = vec4_abs(x);
	__m128 one_over_x = _mm_div_ps(ones, x);
	__m128 abs_one_over_x = vec4_abs(one_over_x);
	__m128 min_abs = _mm_min_ps(abs_x, abs_one_over_x);
	__m128 cmp = _mm_cmpge_ps(min_abs, tol);
	return _mm_and_ps(cmp, one_over_x);
}

// ----------------------------------------------------------------------------

static void svd_pseudoinverse(Mat4x4& o, const __m128& sigma, const Mat4x4& v)
{
	const __m128 invdet = svd_invdet(sigma);

	Mat4x4 m;
	m.row[0] = _mm_mul_ps(v.row[0], invdet);
	m.row[1] = _mm_mul_ps(v.row[1], invdet);
	m.row[2] = _mm_mul_ps(v.row[2], invdet);
	m.row[3] = _mm_set1_ps(0.f);

	o.row[0].m128_f32[0] = vec4_dot(m.row[0], v.row[0]);
	o.row[0].m128_f32[1] = vec4_dot(m.row[1], v.row[0]);
	o.row[0].m128_f32[2] = vec4_dot(m.row[2], v.row[0]);
	o.row[0].m128_f32[3] = 0.f;

	o.row[1].m128_f32[0] = vec4_dot(m.row[0], v.row[1]);
	o.row[1].m128_f32[1] = vec4_dot(m.row[1], v.row[1]);
	o.row[1].m128_f32[2] = vec4_dot(m.row[2], v.row[1]);
	o.row[1].m128_f32[3] = 0.f;

	o.row[2].m128_f32[0] = vec4_dot(m.row[0], v.row[2]);
	o.row[2].m128_f32[1] = vec4_dot(m.row[1], v.row[2]);
	o.row[2].m128_f32[2] = vec4_dot(m.row[2], v.row[2]);
	o.row[2].m128_f32[3] = 0.f;

	o.row[3] = m.row[3];
}

// ----------------------------------------------------------------------------

static void svd_solve_ATA_ATb(const Mat4x4& ATA, const __m128& ATb, __m128& x)
{
	Mat4x4 V;
	V.row[0] = _mm_set_ps(0.f, 0.f, 0.f, 1.f);
	V.row[1] = _mm_set_ps(0.f, 0.f, 1.f, 0.f);
	V.row[2] = _mm_set_ps(0.f, 1.f, 0.f, 0.f);
	V.row[3] = _mm_set_ps(0.f, 0.f, 0.f, 0.f);

	const __m128 sigma = svd_solve_sym(V, ATA);

	// A = UEV^T; U = A / (E*V^T)
	Mat4x4 Vinv;
	svd_pseudoinverse(Vinv, sigma, V);
	x = vec4_mul_m4x4(ATb, Vinv);
}

// ----------------------------------------------------------------------------

void qef_simd_add(
	const __m128& p, const __m128& n,
	Mat4x4& ATA,
	__m128& ATb,
	__m128& pointaccum)
{
	__m128 nX = _mm_mul_ps(_mm_shuffle_ps(n, n, _MM_SHUFFLE(0, 0, 0, 0)), n);
	__m128 nY = _mm_mul_ps(_mm_shuffle_ps(n, n, _MM_SHUFFLE(1, 1, 1, 1)), n);
	__m128 nZ = _mm_mul_ps(_mm_shuffle_ps(n, n, _MM_SHUFFLE(2, 2, 2, 2)), n);

	ATA.row[0] = _mm_add_ps(ATA.row[0], nX);
	ATA.row[1] = _mm_add_ps(ATA.row[1], nY);
	ATA.row[2] = _mm_add_ps(ATA.row[2], nZ);

	const float d = vec4_dot(p, n);
	__m128 x = _mm_set_ps(0.f, d, d, d);
	x = _mm_mul_ps(x, n);
	ATb = _mm_add_ps(ATb, x);
	pointaccum = _mm_add_ps(pointaccum, p);
}

// ----------------------------------------------------------------------------

float qef_simd_calc_error(const Mat4x4& A, const __m128& x, const __m128& b)
{
	__m128 tmp =  vec4_mul_m4x4(x, A);
	tmp = _mm_sub_ps(b, tmp);

	return vec4_dot(tmp, tmp);
}

// ----------------------------------------------------------------------------

float qef_simd_solve(
	const Mat4x4& ATA,
	const __m128& ATb,
	const __m128& pointaccum,
	__m128& x)
{
	const __m128 masspoint = _mm_div_ps(pointaccum, _mm_set1_ps(pointaccum.m128_f32[3]));

	__m128 p = vec4_mul_m4x4(masspoint, ATA);
	p = _mm_sub_ps(ATb, p);

	svd_solve_ATA_ATb(ATA, p, x);

	const float error = qef_simd_calc_error(ATA, x, ATb);
	x = _mm_add_ps(x, masspoint);

	return error;
}

// ----------------------------------------------------------------------------

float qef_solve_from_points(
	const __m128* positions,
	const __m128* normals,
	const int count,
	__m128* solved_position)
{
	__m128 pointaccum = _mm_set1_ps(0.f);
	__m128 ATb = _mm_set1_ps(0.f);

	Mat4x4 ATA;
	ATA.row[0] = _mm_set1_ps(0.f);
	ATA.row[1] = _mm_set1_ps(0.f);
	ATA.row[2] = _mm_set1_ps(0.f);
	ATA.row[3] = _mm_set1_ps(0.f);

	for (int i = 0; i < count && i < QEF_MAX_INPUT_COUNT; i++)
	{
		qef_simd_add(positions[i], normals[i], ATA, ATb, pointaccum);
	}

	__declspec(align(16)) float x[4];
	_mm_store_ps(x, ATb);
	_mm_set_ps(0.f, x[2], x[1], x[0]);

	return qef_simd_solve(ATA, ATb, pointaccum, *solved_position);
}

// ----------------------------------------------------------------------------

float qef_solve_from_points_4d(
	const float* positions,
	const float* normals,
	const int count,
	float* solved_position)
{
	if (count < 2 || count > QEF_MAX_INPUT_COUNT)
	{
		solved_position[0] = solved_position[1] = solved_position[2] = solved_position[3] = 0.f;
		return 0.f;
	}

	__m128 p[QEF_MAX_INPUT_COUNT];
	__m128 n[QEF_MAX_INPUT_COUNT];
	for (int i = 0; i < count; i++)
	{
		p[i] = _mm_load_ps(&positions[i * 4]);
		n[i] = _mm_load_ps(&normals[i * 4]);
	}

	__m128 solved;
	const float error = qef_solve_from_points(p, n, count, &solved);
	_mm_store_ps(solved_position, solved);
	return error;
}

// ----------------------------------------------------------------------------

float qef_solve_from_points_4d_interleaved(
	const float* data,
	const size_t stride,
	const int count,
	float* solved_position)
{
	if (count < 2 || count > QEF_MAX_INPUT_COUNT)
	{
		solved_position[0] = solved_position[1] = solved_position[2] = solved_position[3] = 0.f;
		return 0.f;
	}

	__m128 p[QEF_MAX_INPUT_COUNT];
	__m128 n[QEF_MAX_INPUT_COUNT];
	for (int i = 0; i < count; i++)
	{
		p[i] = _mm_load_ps(&data[(i * stride) + 0]);
		n[i] = _mm_load_ps(&data[(i * stride) + 4]);
	}

	__m128 solved;
	const float error = qef_solve_from_points(p, n, count, &solved);
	_mm_store_ps(solved_position, solved);
	return error;
}

// ----------------------------------------------------------------------------

float qef_solve_from_points_3d(
	const float* positions,
	const float* normals,
	const int count,
	float* solved_position)
{
	if (count < 2 || count > QEF_MAX_INPUT_COUNT)
	{
		solved_position[0] = solved_position[1] = solved_position[2] = solved_position[3] = 0.f;
		return 0.f;
	}

	__m128 p[QEF_MAX_INPUT_COUNT];
	__m128 n[QEF_MAX_INPUT_COUNT];
	for (int i = 0; i < count; i++)
	{
		const float* pos = &positions[i * 3];
		const float* nrm = &normals[i * 3];
		p[i] = _mm_set_ps(1.f, pos[2], pos[1], pos[0]);
		n[i] = _mm_set_ps(0.f, nrm[2], nrm[1], nrm[0]);
	}

	__m128 solved;
	const float error = qef_solve_from_points(p, n, count, &solved);

	solved_position[0] = solved.m128_f32[0];
	solved_position[1] = solved.m128_f32[1];
	solved_position[2] = solved.m128_f32[2];
	return error;
}

#endif // QEF_INCLUDE_IMPL

#endif	//	HAS_QEF_SIMD_H_BEEN_INCLUDED


text:
// minimal SVD implementation for calculating feature points from hermite data
// works in C++ and GLSL

// public domain

#ifndef USE_GLSL
#define USE_GLSL 0
#endif

#ifndef USE_WEBGL
#define USE_WEBGL 0
#endif

#ifndef DEBUG_SVD
#define DEBUG_SVD (!USE_GLSL)
#endif

#ifndef SVD_COMPARE_REFERENCE
#define SVD_COMPARE_REFERENCE 1
#endif

#define SVD_NUM_SWEEPS 5

#if USE_GLSL

// GLSL prerequisites

#define IN(t,x) in t x
#define OUT(t, x) out t x
#define INOUT(t, x) inout t x
#define rsqrt inversesqrt

#define SWIZZLE_XYZ(v) v.xyz

#else

// C++ prerequisites

#include <math.h>
#include <glm/glm.hpp>
#include <glm/gtc/swizzle.hpp>
using namespace glm;

#define abs fabs
#define sqrt sqrtf
#define max(x,y) (((x)>(y))?(x):(y))
#define IN(t,x) t x
#define OUT(t, x) t &x
#define INOUT(t, x) t &x

float rsqrt(float x) {
    return 1.0 / sqrt(x);
}

#define SWIZZLE_XYZ(v) vec3(swizzle<X,Y,Z>(v))

#endif

#if DEBUG_SVD

// Debugging
////////////////////////////////////////////////////////////////////////////////

void dump_vec3(vec3 v) {
    printf("(%.5f %.5f %.5f)\n", v[0], v[1], v[2]);
}

void dump_vec4(vec4 v) {
    printf("(%.5f %.5f %.5f %.5f)\n", v[0], v[1], v[2], v[3]);
}

void dump_mat3(mat3 m) {
    printf("(%.5f %.5f %.5f\n %.5f %.5f %.5f\n %.5f %.5f %.5f)\n",
        m[0][0], m[0][1], m[0][2],
        m[1][0], m[1][1], m[1][2],
        m[2][0], m[2][1], m[2][2]);
}

#endif

// SVD
////////////////////////////////////////////////////////////////////////////////

const float Tiny_Number = 1.e-20;

void givens_coeffs_sym(float a_pp, float a_pq, float a_qq, OUT(float,c), OUT(float,s)) {
    if (a_pq == 0.0) {
        c = 1.0;
        s = 0.0;
        return;
    }
    float tau = (a_qq - a_pp) / (2.0 * a_pq);
    float stt = sqrt(1.0 + tau * tau);
    float tan = 1.0 / ((tau >= 0.0) ? (tau + stt) : (tau - stt));
    c = rsqrt(1.0 + tan * tan);
    s = tan * c;
}

void svd_rotate_xy(INOUT(float,x), INOUT(float,y), IN(float,c), IN(float,s)) {
    float u = x; float v = y;
    x = c * u - s * v;
    y = s * u + c * v;
}

void svd_rotateq_xy(INOUT(float,x), INOUT(float,y), INOUT(float,a), IN(float,c), IN(float,s)) {
    float cc = c * c; float ss = s * s;
    float mx = 2.0 * c * s * a;
    float u = x; float v = y;
    x = cc * u - mx + ss * v;
    y = ss * u + mx + cc * v;
}

#if USE_WEBGL

void svd_rotate01(INOUT(mat3,vtav), INOUT(mat3,v)) {
    if (vtav[0][1] == 0.0) return;

    float c, s;
    givens_coeffs_sym(vtav[0][0], vtav[0][1], vtav[1][1], c, s);
    svd_rotateq_xy(vtav[0][0],vtav[1][1],vtav[0][1],c,s);
    svd_rotate_xy(vtav[0][2], vtav[1][2], c, s);
    vtav[0][1] = 0.0;

    svd_rotate_xy(v[0][0], v[0][1], c, s);
    svd_rotate_xy(v[1][0], v[1][1], c, s);
    svd_rotate_xy(v[2][0], v[2][1], c, s);
}

void svd_rotate02(INOUT(mat3,vtav), INOUT(mat3,v)) {
    if (vtav[0][2] == 0.0) return;

    float c, s;
    givens_coeffs_sym(vtav[0][0], vtav[0][2], vtav[2][2], c, s);
    svd_rotateq_xy(vtav[0][0],vtav[2][2],vtav[0][2],c,s);
    svd_rotate_xy(vtav[0][1], vtav[1][2], c, s);
    vtav[0][2] = 0.0;

    svd_rotate_xy(v[0][0], v[0][2], c, s);
    svd_rotate_xy(v[1][0], v[1][2], c, s);
    svd_rotate_xy(v[2][0], v[2][2], c, s);
}

void svd_rotate12(INOUT(mat3,vtav), INOUT(mat3,v)) {
    if (vtav[1][2] == 0.0) return;

    float c, s;
    givens_coeffs_sym(vtav[1][1], vtav[1][2], vtav[2][2], c, s);
    svd_rotateq_xy(vtav[1][1],vtav[2][2],vtav[1][2],c,s);
    svd_rotate_xy(vtav[0][1], vtav[0][2], c, s);
    vtav[1][2] = 0.0;

    svd_rotate_xy(v[0][1], v[0][2], c, s);
    svd_rotate_xy(v[1][1], v[1][2], c, s);
    svd_rotate_xy(v[2][1], v[2][2], c, s);
}

#else

void svd_rotate(INOUT(mat3,vtav), INOUT(mat3,v), IN(int,a), IN(int,b)) {
    if (vtav[a][b] == 0.0) return;

    float c, s;
    givens_coeffs_sym(vtav[a][a], vtav[a][b], vtav[b][b], c, s);
    svd_rotateq_xy(vtav[a][a],vtav[b][b],vtav[a][b],c,s);
    svd_rotate_xy(vtav[0][3-b], vtav[1-a][2], c, s);
    vtav[a][b] = 0.0;

    svd_rotate_xy(v[0][a], v[0][b], c, s);
    svd_rotate_xy(v[1][a], v[1][b], c, s);
    svd_rotate_xy(v[2][a], v[2][b], c, s);
}

#endif

void svd_solve_sym(IN(mat3,a), OUT(vec3,sigma), INOUT(mat3,v)) {
    // assuming that A is symmetric: can optimize all operations for
    // the upper right triagonal
    mat3 vtav = a;
    // assuming V is identity: you can also pass a matrix the rotations
    // should be applied to
    // U is not computed
    for (int i = 0; i < SVD_NUM_SWEEPS; ++i) {
        svd_rotate(vtav, v, 0, 1);
        svd_rotate(vtav, v, 0, 2);
        svd_rotate(vtav, v, 1, 2);
    }
    sigma = vec3(vtav[0][0],vtav[1][1],vtav[2][2]);
}

float svd_invdet(float x, float tol) {
    return (abs(x) < tol || abs(1.0 / x) < tol) ? 0.0 : (1.0 / x);
}

void svd_pseudoinverse(OUT(mat3,o), IN(vec3,sigma), IN(mat3,v)) {
    float d0 = svd_invdet(sigma[0], Tiny_Number);
    float d1 = svd_invdet(sigma[1], Tiny_Number);
    float d2 = svd_invdet(sigma[2], Tiny_Number);
    o = mat3(v[0][0] * d0 * v[0][0] + v[0][1] * d1 * v[0][1] + v[0][2] * d2 * v[0][2],
             v[0][0] * d0 * v[1][0] + v[0][1] * d1 * v[1][1] + v[0][2] * d2 * v[1][2],
             v[0][0] * d0 * v[2][0] + v[0][1] * d1 * v[2][1] + v[0][2] * d2 * v[2][2],
             v[1][0] * d0 * v[0][0] + v[1][1] * d1 * v[0][1] + v[1][2] * d2 * v[0][2],
             v[1][0] * d0 * v[1][0] + v[1][1] * d1 * v[1][1] + v[1][2] * d2 * v[1][2],
             v[1][0] * d0 * v[2][0] + v[1][1] * d1 * v[2][1] + v[1][2] * d2 * v[2][2],
             v[2][0] * d0 * v[0][0] + v[2][1] * d1 * v[0][1] + v[2][2] * d2 * v[0][2],
             v[2][0] * d0 * v[1][0] + v[2][1] * d1 * v[1][1] + v[2][2] * d2 * v[1][2],
             v[2][0] * d0 * v[2][0] + v[2][1] * d1 * v[2][1] + v[2][2] * d2 * v[2][2]);
}

void svd_solve_ATA_ATb(
    IN(mat3,ATA), IN(vec3,ATb), OUT(vec3,x)
) {
    mat3 V = mat3(1.0);
    vec3 sigma;

    svd_solve_sym(ATA, sigma, V);

    // A = UEV^T; U = A / (E*V^T)
#if DEBUG_SVD

    printf("ATA="); dump_mat3(ATA);
    printf("ATb="); dump_vec3(ATb);
    printf("V="); dump_mat3(V);
    printf("sigma="); dump_vec3(sigma);
#endif
    mat3 Vinv;
    svd_pseudoinverse(Vinv, sigma, V);
    x = Vinv * ATb;

#if DEBUG_SVD
    printf("Vinv="); dump_mat3(Vinv);
#endif
}

vec3 svd_vmul_sym(IN(mat3,a), IN(vec3,v)) {
    return vec3(
        dot(a[0],v),
        (a[0][1] * v.x) + (a[1][1] * v.y) + (a[1][2] * v.z),
        (a[0][2] * v.x) + (a[1][2] * v.y) + (a[2][2] * v.z)
    );
}

void svd_mul_ata_sym(OUT(mat3,o), IN(mat3,a))
{
    o[0][0] = a[0][0] * a[0][0] + a[1][0] * a[1][0] + a[2][0] * a[2][0];
    o[0][1] = a[0][0] * a[0][1] + a[1][0] * a[1][1] + a[2][0] * a[2][1];
    o[0][2] = a[0][0] * a[0][2] + a[1][0] * a[1][2] + a[2][0] * a[2][2];
    o[1][1] = a[0][1] * a[0][1] + a[1][1] * a[1][1] + a[2][1] * a[2][1];
    o[1][2] = a[0][1] * a[0][2] + a[1][1] * a[1][2] + a[2][1] * a[2][2];
    o[2][2] = a[0][2] * a[0][2] + a[1][2] * a[1][2] + a[2][2] * a[2][2];
}

void svd_solve_Ax_b(IN(mat3,a), IN(vec3,b), OUT(mat3,ATA), OUT(vec3,ATb), OUT(vec3,x)) {
    svd_mul_ata_sym(ATA, a);
    ATb = b * a; // transpose(a) * b;
    svd_solve_ATA_ATb(ATA, ATb, x);
}

// QEF
////////////////////////////////////////////////////////////////////////////////

void qef_add(
    IN(vec3,n), IN(vec3,p),
    INOUT(mat3,ATA),
    INOUT(vec3,ATb),
    INOUT(vec4,pointaccum))
{
#if DEBUG_SVD
    printf("+plane=");dump_vec4(vec4(n, dot(-p,n)));
#endif
    ATA[0][0] += n.x * n.x;
    ATA[0][1] += n.x * n.y;
    ATA[0][2] += n.x * n.z;
    ATA[1][1] += n.y * n.y;
    ATA[1][2] += n.y * n.z;
    ATA[2][2] += n.z * n.z;

    float b = dot(p, n);
    ATb += n * b;
    pointaccum += vec4(p,1.0);
}

float qef_calc_error(IN(mat3,A), IN(vec3, x), IN(vec3, b)) {
    vec3 vtmp = b - svd_vmul_sym(A, x);
    return dot(vtmp,vtmp);
}

float qef_solve(
    IN(mat3,ATA),
    IN(vec3,ATb),
    IN(vec4,pointaccum),
    OUT(vec3,x)
) {
    vec3 masspoint = SWIZZLE_XYZ(pointaccum) / pointaccum.w;
    ATb -= svd_vmul_sym(ATA, masspoint);
    svd_solve_ATA_ATb(ATA, ATb, x);
    float result = qef_calc_error(ATA, x, ATb);

    x += masspoint;

    return result;
}

// Test
////////////////////////////////////////////////////////////////////////////////

#if USE_GLSL

void main(void) {
}

#else

#if DEBUG_SVD
#undef sqrt
#undef abs
#undef rsqrt
#undef max

#if SVD_COMPARE_REFERENCE
#include "qr_solve.h"
#include "r8lib.h"
#include "qr_solve.c"
#include "r8lib.c"
#endif

int main(void) {
    vec4 pointaccum = vec4(0.0);
    mat3 ATA = mat3(0.0);
    vec3 ATb = vec3(0.0);

    #define COUNT 5
    vec3 normals[COUNT] = {
        normalize(vec3( 1.0,1.0,0.0)),
        normalize(vec3( 1.0,1.0,0.0)),
        normalize(vec3(-1.0,1.0,0.0)),
        normalize(vec3(-1.0,2.0,1.0)),
        //normalize(vec3(-1.0,1.0,0.0)),
        normalize(vec3(-1.0,1.0,0.0)),
    };
    vec3 points[COUNT] = {
        vec3(  1.0,0.0,0.3),
        vec3(  0.9,0.1,-0.5),
        vec3( -0.8,0.2,0.6),
        vec3( -1.0,0.0,0.01),
        vec3( -1.1,-0.1,-0.5),
    };

    for (int i= 0; i < COUNT; ++i) {
        qef_add(normals[i],points[i],ATA,ATb,pointaccum);
    }
    vec3 com = SWIZZLE_XYZ(pointaccum) / pointaccum.w;

    vec3 x;
    float error = qef_solve(ATA,ATb,pointaccum,x);

    printf("masspoint = (%.5f %.5f %.5f)\n", com.x, com.y, com.z);
    printf("point = (%.5f %.5f %.5f)\n", x.x, x.y, x.z);
    printf("error = %.5f\n", error);

#if SVD_COMPARE_REFERENCE
    double a[COUNT*3];
    double b[COUNT];

    for (int i = 0; i < COUNT; ++i) {
        b[i] = (points[i].x - com.x)*normals[i].x
             + (points[i].y - com.y)*normals[i].y
             + (points[i].z - com.z)*normals[i].z;
        a[i] = normals[i].x;
        a[i+COUNT] = normals[i].y;
        a[i+2*COUNT] = normals[i].z;
    }

    double *c = svd_solve(5,3,a,b,0.1);

    vec3 result = com + vec3(c[0], c[1], c[2]);
    r8_free(c);
    printf("reference="); dump_vec3(result);
#endif

    return 0;
}
#endif

#endif

repo_name: nickgildea/qef

repo_language: C

mime_type: text/plain


text:
# qef
Quadric Error Function (QEF) / Single Value Decomposition solvers in various languages


text:
package sudoku

import (
	"strings"
	"testing"
)

func FuzzParseboard(f *testing.F) {
	f.Add("4.....8.5.3..........7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......")
	f.Add(DisplayAsInput(EmptyBoard()))
	f.Add(`
. . . |. . 6 |. . .
. 5 9 |. . . |. . 8
2 . . |. . 8 |. . .
------+------+------
. 4 5 |. . . |. . .
. . 3 |. . . |. . .
. . 6 |. . 3 |. 5 4
------+------+------
. . . |3 2 5 |. . 6
. . . |. . . |. . .
. . . |. . . |. . .`)
	f.Add("1234")
	f.Add("")
	f.Add(strings.Repeat("1", 81))

	f.Fuzz(func(t *testing.T, bstr string) {
		// Mostly checks that there are no panics, and ParseBoard returns either
		// a board or and error, but not both.
		b, err := ParseBoard(bstr, true)
		if b != nil && err != nil {
			t.Fatalf("expect b or err to be nil, got b=%v, err=%v", b, err)
		}
		if b == nil && err == nil {
			t.Fatalf("expect b or err to be non-nil, got b=%v, err=%v", b, err)
		}
	})
}

repo_name: eliben/go-sudoku

repo_language: Go

mime_type: text/x-c


text:
package sudoku

import (
	"log"
	"math/rand"
	"testing"
	"time"

	"golang.org/x/exp/slices"
)

func TestEvaluateDifficulty(t *testing.T) {
	rand.Seed(time.Now().UnixNano())
	getDifficulty := func(board string) float64 {
		v, err := ParseBoard(board, false)
		if err != nil {
			t.Fatal(err)
		}
		vcopy := slices.Clone(v)

		d, err := EvaluateDifficulty(v)
		if err != nil {
			t.Fatal(err)
		}

		if !slices.Equal(v, vcopy) {
			t.Errorf("EvaluateDifficulty modified values")
		}

		return d
	}

	easyD := getDifficulty(easyboard1)
	hardD := getDifficulty(hardboard2)
	hardlongD := getDifficulty(hardlong)

	if easyD > hardD || hardD > hardlongD {
		t.Errorf("got easyD: %v, hardD: %v, hardlongD: %v", easyD, hardD, hardlongD)
	}
}

var filled string = `
3 4 5 |7 9 2 |6 1 8
8 9 6 |3 5 1 |2 7 4
1 7 2 |8 6 4 |9 5 3
------+------+------
2 3 8 |6 4 5 |1 9 7
7 5 4 |9 1 8 |3 2 6
6 1 9 |2 3 7 |4 8 5
------+------+------
5 6 7 |4 2 9 |8 3 1
4 2 1 |5 8 3 |7 6 9
9 8 3 |1 7 6 |5 4 2
`

func TestFilled(t *testing.T) {
	v, err := ParseBoard(filled, false)
	if err != nil {
		log.Fatal(err)
	}

	d, err := EvaluateDifficulty(v)
	if err != nil {
		log.Fatal(err)
	}

	if d != 1.0 {
		t.Errorf("got d=%v; expect difficulty of filled board to be 1.0", d)
	}
}


text:
package sudoku

import (
	"fmt"

	"golang.org/x/exp/slices"
)

// EvaluateDifficulty evaluates the difficulty of a Sudoku puzzle heuristically
// and returns the score on a scale from 1.0 (easiest) to 5.0 hardest. It can
// also return an error if the given board has contradictions, is unsolvable,
// etc. It should be passed a board that didn't have elimination applied to it.
//
// The heuristics are based on 4 factors:
//
// 1. How many hints (filled-in squares) the board has.
// 2. How many hints remain after running a round of elimination (first-order
//    Sudoku solving value deduction).
// 3. How many hints does a row or column with the minimal number of hints have
// 4. How many guesses a backtracking search requires to solve the board
//    (averaged over multiple runs).
//
// This approach was partially inspired by the paper "Sudoku Puzzles Generating:
// from Easy to Evil" by Xiang-Sun ZHANG's research group.
func EvaluateDifficulty(values Values) (float64, error) {
	hintsBeforeElimination := CountHints(values)

	// Count the lower bound (minimal number) of hints in individual rows and
	// cols, pre elimination.
	minHints := 9

	index := func(row, col int) Index {
		return row*9 + col
	}

	// ... first the rows.
	for row := 0; row < 9; row++ {
		rowCount := 0
		for col := 0; col < 9; col++ {
			if values[index(row, col)].Size() == 1 {
				rowCount++
			}
		}
		if rowCount < minHints {
			minHints = rowCount
		}
	}

	// ... then the columns.
	for col := 0; col < 9; col++ {
		colCount := 0
		for row := 0; row < 9; row++ {
			if values[index(row, col)].Size() == 1 {
				colCount++
			}
		}
		if colCount < minHints {
			minHints = colCount
		}
	}

	// Run elimination and count how many hints are on the board after it.
	vcopy := slices.Clone(values)
	if !EliminateAll(vcopy) {
		return 0, fmt.Errorf("contradiction in board")
	}
	hintsAfterElimination := CountHints(vcopy)

	// Run a number of randomized searches and count the average search count.
	EnableStats = true
	var totalSearches uint64 = 0
	iterations := 10
	for i := 0; i < iterations; i++ {
		Stats.Reset()
		_, solved := Solve(vcopy, SolveOptions{Randomize: true})
		if !solved {
			return 0, fmt.Errorf("cannot solve")
		}
		totalSearches += Stats.NumSearches
	}
	EnableStats = false
	averageSearches := float64(totalSearches) / float64(iterations)

	// Assign difficulty scores based on ranges in each category.
	var hintsBeforeDifficulty float64
	if hintsBeforeElimination > 50 {
		hintsBeforeDifficulty = 1.0
	} else if hintsBeforeElimination > 35 {
		hintsBeforeDifficulty = 2.0
	} else if hintsBeforeElimination > 31 {
		hintsBeforeDifficulty = 3.0
	} else if hintsBeforeElimination > 27 {
		hintsBeforeDifficulty = 4.0
	} else {
		hintsBeforeDifficulty = 5.0
	}

	var hintsAfterDifficulty float64
	if hintsAfterElimination > 55 {
		hintsAfterDifficulty = 1.0
	} else if hintsAfterElimination > 42 {
		hintsAfterDifficulty = 2.0
	} else if hintsAfterElimination > 37 {
		hintsAfterDifficulty = 3.0
	} else if hintsAfterElimination > 33 {
		hintsAfterDifficulty = 4.0
	} else {
		hintsAfterDifficulty = 5.0
	}

	var minHintsDifficulty float64
	if minHints >= 5 {
		minHintsDifficulty = 1.0
	} else if minHints == 4 {
		minHintsDifficulty = 2.0
	} else if minHints == 3 {
		minHintsDifficulty = 3.0
	} else if minHints >= 1 {
		minHintsDifficulty = 4.0
	} else {
		minHintsDifficulty = 5.0
	}

	var searchDifficulty float64
	if averageSearches <= 1.0 {
		searchDifficulty = 1.0
	} else if averageSearches < 3.0 {
		searchDifficulty = 2.0
	} else if averageSearches < 10.0 {
		searchDifficulty = 3.0
	} else if averageSearches < 40.0 {
		searchDifficulty = 4.0
	} else {
		searchDifficulty = 5.0
	}

	// Assign final difficulty with weights
	difficulty := 0.5*hintsAfterDifficulty +
		0.3*hintsBeforeDifficulty +
		0.05*minHintsDifficulty +
		0.15*searchDifficulty

	return difficulty, nil
}


text:
package sudoku

import (
	"log"
	"math/rand"
)

// Generate generates a random Sudoku board that has a single solution, with
// at-most hintCount hints remaining on the board. Note that this cannot be
// always reliably done when the count is low (lower than 23 or so), because
// generating a board with a single solution that has a low number of initial
// hints is very hard.
// There are no guarantees made about the difficulty of the generated board,
// though higher hint counts generally correlate with easier boards. It's
// recommended to generate a large number of boards using this function and
// evaluate their difficulty separately using EvaluateDifficulty.
// Notes:
//  * Make sure the default rand source is seeded if you really want to get
//    random boards.
//  * This function may take a while to run when given a low hintCount.
func Generate(hintCount int) Values {
	empty := EmptyBoard()
	board, solved := Solve(empty, SolveOptions{Randomize: true})
	if !solved || !IsSolved(board) {
		log.Fatal("unable to generate solved board from empty")
	}

	removalOrder := rand.Perm(81)
	count := 81

	for _, sq := range removalOrder {
		savedDigit := board[sq]
		// Try to remove the number from square sq.
		board[sq] = FullDigitsSet()

		solutions := SolveAll(board, 2)
		switch len(solutions) {
		case 0:
			// Some sort of bug, because removing a square from a solved board should
			// never result in an unsolvable board.
			log.Fatal("got a board without solutions")
		case 1:
			count--
			if count <= hintCount {
				return board
			}
		default:
			// The board has multiple solutions with this square emptied, so put it
			// back and try again with the next square.
			board[sq] = savedDigit
		}
	}

	return board
}

// GenerateSymmetrical is similar to Generate, but it generates symmetrical
// boards with 180-degree rotational symmetry.
// Because of this additional constraint, it may have more trouble generating
// boards with a small hintCount than Generate, so you'll have to run it more
// times in a loop to find a good low-hint-count board.
func GenerateSymmetrical(hintCount int) Values {
	empty := EmptyBoard()
	board, solved := Solve(empty, SolveOptions{Randomize: true})
	if !solved || !IsSolved(board) {
		log.Fatal("unable to generate solved board from empty")
	}

	// This function works just like Generate, but instead of picking a random
	// square out of all 81, it picks a random square from the first half of the
	// board and then attempts to remove both this square and its reflection.
	removalOrder := rand.Perm(41)
	count := 81

	for _, sq := range removalOrder {
		// Find sq's reflection; note that in the middle row reflectSq could equal
		// sq - we take this into account when counting how many hints remain on
		// the board.
		reflectSq := 80 - sq

		savedDigit := board[sq]
		savedReflect := board[reflectSq]

		board[sq] = FullDigitsSet()
		board[reflectSq] = FullDigitsSet()

		solutions := SolveAll(board, 2)
		switch len(solutions) {
		case 0:
			log.Fatal("got a board without solutions")
		case 1:
			// We may have removed just one or two hints.
			count--
			if sq != reflectSq {
				count--
			}
			if count <= hintCount {
				return board
			}
		default:
			board[sq] = savedDigit
			board[reflectSq] = savedReflect
		}
	}

	return board
}


text:
package sudoku

import (
	"fmt"
	"testing"
)

func TestIsMember(t *testing.T) {
	d := Digits(0b0000000001100010)

	var i uint16
	for i = 1; i <= 9; i++ {
		if i == 1 || i == 5 || i == 6 {
			if !d.IsMember(i) {
				t.Errorf("got isMember(%016b, %d)=false, want true", d, i)
			}
		} else {
			if d.IsMember(i) {
				t.Errorf("got isMember(%016b, %d)=true, want false", d, i)
			}
		}
	}
}

func TestAdd(t *testing.T) {
	d := Digits(0b0000000001100010)

	if d.String() != "156" {
		t.Errorf("got %v, want 156", d.String())
	}

	wantd2 := Digits(0b0000000001100110)
	if d.Add(2) != wantd2 {
		t.Errorf("got wantd2=%v, want=%v", d.Add(2), wantd2)
	}

	wantd5 := Digits(0b0000000001100010)
	if d.Add(5) != wantd5 {
		t.Errorf("got wantd5=%v, want=%v", d.Add(5), wantd5)
	}
}

func TestRemove(t *testing.T) {
	d := Digits(0b0000000001100010)

	wantd5 := Digits(0b0000000001000010)
	if d.Remove(5) != wantd5 {
		t.Errorf("got wantd5=%v, want=%v", d.Remove(5), wantd5)
	}

	wantd8 := Digits(0b0000000001100010)
	if d.Remove(8) != wantd8 {
		t.Errorf("got wantd8=%v, want=%v", d.Remove(8), wantd8)
	}

	wantd6no1 := Digits(0b0000000001000000)
	donly6 := d.Remove(1).Remove(5)
	if donly6 != wantd6no1 {
		t.Errorf("got wantd6no1=%v, want=%v", donly6, wantd6no1)
	}

	if donly6.String() != "6" {
		t.Errorf("got %v, want 6", donly6.String())
	}
}

func TestRemoveAll(t *testing.T) {
	d := Digits(0b0011100110)

	got1 := d.RemoveAll(0b0011000000)
	want1 := Digits(0b0000100110)
	if got1 != want1 {
		t.Errorf("got %v, want %v", got1, want1)
	}

	got2 := d.RemoveAll(0b0001000110)
	want2 := Digits(0b0010100000)
	if got2 != want2 {
		t.Errorf("got %v, want %v", got2, want2)
	}
}

func TestAddRemoveAllSize(t *testing.T) {
	// Exhaustive testing that adds/removes every digits and tests that IsMember
	// also keeps working.

	// Start with zero. Each iteration adds one digit, tests membership, then
	// removes the digit and tests again.
	d := Digits(0)

	testNoMembers := func() {
		if d.Size() != 0 {
			t.Errorf("got size=%v, want 0", d.Size())
		}

		for dig := uint16(1); dig <= 9; dig++ {
			if d.IsMember(dig) {
				t.Errorf("got IsMember=true for %v, want false", dig)
			}
		}
	}
	testNoMembers()

	for dig := uint16(1); dig <= 9; dig++ {
		t.Run(fmt.Sprintf("dig=%v", dig), func(t *testing.T) {
			// Add 'dig' to set
			d = d.Add(dig)

			if d.Size() != 1 {
				t.Errorf("got size=%v, want 1", d.Size())
			}

			off := d.SingleMemberDigit()
			if off != dig {
				t.Errorf("got SingleMemberDigit=%v, want %v", off, dig)
			}

			// For each 'dig2', check set membership
			for dig2 := uint16(1); dig2 <= 9; dig2++ {
				if dig2 == dig {
					if !d.IsMember(dig2) {
						t.Errorf("got IsMember=false for %v, want true", dig2)
					}
				} else {
					if d.IsMember(dig2) {
						t.Errorf("got IsMember=true for %v, want false", dig2)
					}
				}
			}

			d = d.Remove(dig)
			testNoMembers()
		})
	}
}

func TestFullDigitsSet(t *testing.T) {
	d := FullDigitsSet()
	for dig := uint16(1); dig <= 9; dig++ {
		if !d.IsMember(dig) {
			t.Errorf("got IsMember=false for %v, want true", dig)
		}
	}

	if d.String() != "123456789" {
		t.Errorf("got %v, want all digits", d.String())
	}
}

func TestSingleDigitSet(t *testing.T) {
	d := SingleDigitSet(5)
	for dig := uint16(1); dig <= 9; dig++ {
		if dig == 5 {
			if !d.IsMember(dig) {
				t.Errorf("got IsMember=false for 5, want true")
			}
		} else {
			if d.IsMember(dig) {
				t.Errorf("got IsMember=true for %v, want false", dig)
			}
		}
	}
}

func TestTwoMemberDigits(t *testing.T) {
	d := Digits(0b0000000000100100)
	d1, d2 := d.twoMemberDigits()
	if d1 != 2 || d2 != 5 {
		t.Errorf("got %v,%v, want 2 and 5", d1, d2)
	}

	d = Digits(0b0000001000000010)
	d1, d2 = d.twoMemberDigits()
	if d1 != 1 || d2 != 9 {
		t.Errorf("got %v,%v, want 1 and 9", d1, d2)
	}
}


text:
package sudoku

import (
	"fmt"
	"math/bits"
	"strings"
)

// Digits represents a set of possible digits for a Sudoku square. The functions
// in this file perform set operations on Digits, as needed for Sudoku.
type Digits uint16

// Internal representation: digit N is represented by the Nth lowest bit in
// the Digits value, e.g.:
//
//    Digits = 0b0000_0000_0110_0010
//
// has the bits N=1,5,6 are set, so it represents the set of digits {1, 5, 6}

// FullDigitsSet returns a Digits with all possible digits set.
func FullDigitsSet() Digits {
	return 0b0000001111111110
}

// SingleDigitSet returns a Digits with a single digit 'n' set.
func SingleDigitSet(n uint16) Digits {
	return 1 << n
}

// IsMember checks whether digit n is a member of the digit set d.
func (d Digits) IsMember(n uint16) bool {
	return (d & (1 << n)) != 0
}

// Add adds digit n to set d and returns the new set.
func (d Digits) Add(n uint16) Digits {
	return d | (1 << n)
}

// Remove removes digit n from set d and returns the new set.
func (d Digits) Remove(n uint16) Digits {
	return d &^ (1 << n)
}

// RemoveAll removes all digits represented by dn from d and returns
// the new set.
func (d Digits) RemoveAll(dn Digits) Digits {
	return d &^ dn
}

// Size returns the size of the set - the number of digits in it.
func (d Digits) Size() int {
	return bits.OnesCount16(uint16(d))
}

// SingleMemberDigit returns the digit that's a member of a 1-element set; this
// assumes that the set indeed has a single element.
func (d Digits) SingleMemberDigit() uint16 {
	return uint16(bits.TrailingZeros16(uint16(d)))
}

// twoMemberDigits returns the only two digits that are member of a 2-element
// set; this assumes that the set indeed has two elements.
func (d Digits) twoMemberDigits() (uint16, uint16) {
	d1 := uint16(bits.TrailingZeros16(uint16(d)))
	d2 := 16 - uint16(bits.LeadingZeros16(uint16(d))) - 1
	return d1, d2
}

// String implements the fmt.Stringer interface for Digits.
func (d Digits) String() string {
	var parts []string
	for i := uint16(1); i <= uint16(9); i++ {
		if d.IsMember(i) {
			parts = append(parts, fmt.Sprintf("%v", i))
		}
	}
	return strings.Join(parts, "")
}


text:
module github.com/eliben/go-sudoku

go 1.18

require golang.org/x/exp v0.0.0-20220613132600-b0d781184e0d


text:
package sudoku

import (
	"fmt"
	"log"
	"math/rand"
	"strings"
	"testing"
	"time"

	"golang.org/x/exp/slices"
)

func TestInit(t *testing.T) {
	// Smoke testing for the top-level vars initialized in init()
	if len(unitlist) != 27 {
		t.Errorf("got len=%v, want 27", len(unitlist))
	}

	wantUnits := []Unit{
		{18, 19, 20, 21, 22, 23, 24, 25, 26},
		{2, 11, 20, 29, 38, 47, 56, 65, 74},
		{0, 1, 2, 9, 10, 11, 18, 19, 20}}

	if !slices.EqualFunc(wantUnits, units[20], func(a, b Unit) bool {
		return slices.Equal(a, b)
	}) {
		t.Errorf("got units[20]=%v\nwant %v", units[20], wantUnits)
	}

	gotPeers := peers[20]
	slices.Sort(gotPeers)
	wantPeers := []Index{0, 1, 2, 9, 10, 11, 18, 19, 21, 22, 23, 24, 25, 26, 29, 38, 47, 56, 65, 74}
	if !slices.Equal(wantPeers, gotPeers) {
		t.Errorf("got peers[20]=%v\n want %v", peers[20], wantPeers)
	}
}

func TestAssignElimination(t *testing.T) {
	vals := EmptyBoard()

	if IsSolved(vals) {
		t.Errorf("an empty board is solved")
	}

	// Assign a digit to square 20; check that this digit is the only candidate
	// in square 20, and that it was eliminated from all the peers of 20.
	assign(vals, 20, 5)

	if vals[20].Size() != 1 || vals[20].SingleMemberDigit() != 5 {
		t.Errorf("got vals[20]=%v", vals[20])
	}

	for sq := 0; sq <= 80; sq++ {
		if slices.Contains(peers[20], sq) {
			if vals[sq].IsMember(5) {
				t.Errorf("got member 5 in peer square %v", sq)
			}
		} else {
			if !vals[sq].IsMember(5) {
				t.Errorf("got no member 5 in non-peer square %v", sq)
			}
		}
	}
}

// Easy board from Norvig's example that's solved by constraint propagation
// w/o any search.
var easyboard1 string = "003020600900305001001806400008102900700000008006708200002609500800203009005010300"

// Hard board from Norvig's example.
var hardboard1 string = "4.....8.5.3..........7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......"
var hardboard2 string = "..53.....8......2..7..1.5..4....53...1..7...6..32...8..6.5....9..4....3......97.."

// This is the program-generated sudoku board Norvig reported takes the longest
// for his solver to crack. Note that this board has a very large number of
// solutions.
var hardlong string = `
. . . |. . 6 |. . .
. 5 9 |. . . |. . 8
2 . . |. . 8 |. . .
------+------+------
. 4 5 |. . . |. . .
. . 3 |. . . |. . .
. . 6 |. . 3 |. 5 4
------+------+------
. . . |3 2 5 |. . 6
. . . |. . . |. . .
. . . |. . . |. . .`

func TestParseBoard(t *testing.T) {
	// Need to test the "false" variant of ParseBoard too
	v, err := ParseBoard(easyboard1, true)
	if err != nil {
		t.Fatal(err)
	}

	if !IsSolved(v) {
		t.Errorf("expect easy board to be solved w/ elimination")
	}

	// Without runElimination, the easy board will not be immediately solved by
	// parsing.
	ve, err := ParseBoard(easyboard1, false)
	if err != nil {
		t.Fatal(err)
	}

	if IsSolved(ve) {
		t.Errorf("expect easy board to not be solved w/o elimination")
	}

	if !EliminateAll(ve) || !IsSolved(ve) {
		t.Errorf("expect easy board to be solved with explicit elimination")
	}

	// Harder board that isn't fully solved without search.
	v2, err := ParseBoard(hardboard1, true)
	if err != nil {
		t.Fatal(err)
	}

	if IsSolved(v2) {
		t.Errorf("expect hard board to not be solved")
	}

	// Count how many squares are solved immediately in this puzzle and compare
	// to the number Norvig got.
	var solvedSquares int
	for _, d := range v2 {
		if d.Size() == 1 {
			solvedSquares++
		}
	}

	if solvedSquares != 20 {
		t.Errorf("got %v solved squares, want 20", solvedSquares)
	}
}

func TestSolveBoard(t *testing.T) {
	v, err := ParseBoard(hardboard1, true)
	if err != nil {
		log.Fatal(err)
	}
	vcopy := slices.Clone(v)
	vs, success := Solve(v)
	if !slices.Equal(v, vcopy) {
		t.Errorf("Solve modified board; before=%v, after=%v", vcopy, v)
	}

	if !success || !IsSolved(vs) {
		t.Errorf("expect hardboard1 to be solved by search")
	}

	// Should work on the easy board also (even though it's solved with the
	// initial parse)
	v2, err := ParseBoard(easyboard1, true)
	if err != nil {
		log.Fatal(err)
	}
	v2, success2 := Solve(v2)

	if !success2 || !IsSolved(v2) {
		t.Errorf("expect easy board to be solved by search")
	}

	// And the other hard board
	v3, err := ParseBoard(hardboard2, true)
	if err != nil {
		log.Fatal(err)
	}
	v3, success3 := Solve(v3)

	if !success3 || !IsSolved(v3) {
		t.Errorf("expect hardboard2 to be solved by search")
	}
}

func TestSolveWithStats(t *testing.T) {
	// The easy board is solved just by calling ParseBoard, needing no search.
	WithStats(func() {
		_, err := ParseBoard(easyboard1, true)
		if err != nil {
			t.Fatal(err)
		}

		if Stats.NumAssigns == 0 {
			t.Errorf("got NumAssigns==0")
		}
		if Stats.NumSearches != 0 {
			t.Errorf("got NumSearches=%v, want 0", Stats.NumSearches)
		}

		// For the hard board, we'll find both assigns and searches
		Stats.Reset()

		v, err := ParseBoard(hardboard1, true)
		if err != nil {
			t.Fatal(err)
		}
		_, _ = Solve(v)

		if Stats.NumAssigns == 0 {
			t.Errorf("got NumAssigns==0")
		}
		if Stats.NumSearches == 0 {
			t.Errorf("got NumSearches==0")
		}
	})
}

func TestIsSolved(t *testing.T) {
	v, err := ParseBoard(easyboard1, true)
	if err != nil {
		t.Fatal(err)
	}

	if !IsSolved(v) {
		t.Errorf("expect easy board to be solved")
	}

	// Now modify the board and make sure it's not considered "solved" any more.
	// ... modify by trying to add options to each square separately.
	for sq := range v {
		vcopy := slices.Clone(v)
		vcopy[sq] = vcopy[sq].Add(6).Add(8)

		if IsSolved(vcopy) {
			t.Errorf("expect board to not be solved after modification: %v", vcopy)
		}
	}
}

func TestSolveAll(t *testing.T) {
	v, err := ParseBoard(hardboard1, true)
	if err != nil {
		log.Fatal(err)
	}
	vcopy := slices.Clone(v)
	vs := SolveAll(v, -1)

	if len(vs) != 1 {
		t.Errorf("got %v solutions, want 1", len(vs))
	}
	if !IsSolved(vs[0]) {
		t.Errorf("got %v, want solved board", vs[0])
	}

	if !slices.Equal(v, vcopy) {
		t.Errorf("SolveAll modified its input values")
	}

	// Now generate a multiple-solution board, by replacing all instances
	// of 1 and 2 in the solved board by the digits set "12", permitting any
	// combination of them to solve the board.
	board := vs[0]
	for sq, d := range board {
		if d.Size() == 1 && (d.IsMember(1) || d.IsMember(2)) {
			board[sq] = d.Add(1).Add(2)
		}
	}

	vs = SolveAll(board, -1)
	if len(vs) != 2 {
		t.Errorf("got %v solved boards, want 2", len(vs))
	}

	if !IsSolved(vs[0]) || !IsSolved(vs[1]) {
		t.Errorf("got unsolved boards")
	}

	// Now try to limit max=1, see that SolveAll only returns a single solution.
	vs = SolveAll(board, 1)
	if len(vs) != 1 {
		t.Errorf("got %v solved boards, want 1", len(vs))
	}

	if !IsSolved(vs[0]) {
		t.Errorf("got unsolved boards")
	}

	// Create a board with a contradiction on purpose, and verify that SolveAll
	// returns an empty list.
	v, err = ParseBoard(hardboard1, true)
	if err != nil {
		log.Fatal(err)
	}
	v[30] = SingleDigitSet(1)
	v[31] = SingleDigitSet(2)
	v[32] = SingleDigitSet(3)
	vs = SolveAll(v, -1)
	if len(vs) != 0 {
		t.Errorf("expect unsolvable, got %v", vs)
	}
}

func TestHardlong(t *testing.T) {
	v, err := ParseBoard(hardlong, true)
	if err != nil {
		log.Fatal(err)
	}

	// The "hardlong" puzzle has multiple solutions. Norvig says he found 13
	// different solutions, but there are vastly more. Use SolveAll to explore
	// the first 1000 or so.

	// Find the first 1000 solutions
	vs := SolveAll(v, 1000)

	if len(vs) < 1000 {
		t.Errorf("got %v solutions, expected at least 1000", len(vs))
	}

	for _, v := range vs {
		if !IsSolved(v) {
			t.Errorf("got unsolved board %v", v)
		}
	}
}

func TestApplyTwinsStrategy(t *testing.T) {
	// Basic test: on an empty board, leave only candidates 38 for two squares
	// and verify that twin elimination removed what was needed.

	v := EmptyBoard()
	d38 := Digits(0).Add(3).Add(8)
	v[30] = d38
	v[31] = d38

	ApplyTwinsStrategy(v)

	for _, sq := range []Index{27, 28, 39, 32, 33, 34, 35, 39, 40, 41, 48, 49, 50} {
		if v[sq].IsMember(3) || v[sq].IsMember(8) {
			t.Errorf("got board[%v]=%s, expect no 3 or 8", sq, v[sq])
		}
	}
}

// This board is unsolvable, and it takes the search a while to figure this
// out.
var impossible string = `
. . . |. . 5 |. 8 .
. . . |6 . 1 |. 4 3
. . . |. . . |. . .
------+------+------
. 1 . |5 . . |. . .
. . . |1 . 6 |. . .
3 . . |. . . |. . 5
------+------+------
5 3 . |. . . |. 6 1
. . . |. . . |. . 4
. . . |. . . |. . .`

// Run this test but skip in "not short" mode
func TestImpossible(t *testing.T) {
	if testing.Short() {
		t.Skip("skipping test in short mode.")
	}

	WithStats(func() {
		v, err := ParseBoard(impossible, true)
		if err != nil {
			log.Fatal(err)
		}
		v, success := Solve(v)

		if success || IsSolved(v) {
			t.Errorf("got solved board for impossible")
		}
		fmt.Printf("searches=%v, assigns=%v\n", Stats.NumSearches, Stats.NumAssigns)
	})
}

func TestSolveHardest(t *testing.T) {
	// The "hardest" puzzles Norvig found online (taken from
	// https://norvig.com/hardest.txt)
	hardest := `
85...24..72......9..4.........1.7..23.5...9...4...........8..7..17..........36.4.
..53.....8......2..7..1.5..4....53...1..7...6..32...8..6.5....9..4....3......97..
12..4......5.69.1...9...5.........7.7...52.9..3......2.9.6...5.4..9..8.1..3...9.4
...57..3.1......2.7...234......8...4..7..4...49....6.5.42...3.....7..9....18.....
7..1523........92....3.....1....47.8.......6............9...5.6.4.9.7...8....6.1.
1....7.9..3..2...8..96..5....53..9...1..8...26....4...3......1..4......7..7...3..
1...34.8....8..5....4.6..21.18......3..1.2..6......81.52..7.9....6..9....9.64...2
...92......68.3...19..7...623..4.1....1...7....8.3..297...8..91...5.72......64...
.6.5.4.3.1...9...8.........9...5...6.4.6.2.7.7...4...5.........4...8...1.5.2.3.4.
7.....4...2..7..8...3..8.799..5..3...6..2..9...1.97..6...3..9...3..4..6...9..1.35
....7..2.8.......6.1.2.5...9.54....8.........3....85.1...3.2.8.4.......9.7..6....
`
	for _, board := range strings.Split(hardest, "\n") {
		board = strings.TrimSpace(board)
		if len(board) > 0 {
			v, err := ParseBoard(board, true)
			if err != nil {
				log.Fatalf("error for board %v: %v", board, err)
			}

			// Use ApplyTwinsStrategy here for correctness testing
			if !ApplyTwinsStrategy(v) {
				t.Errorf("ApplyTwinsStrategy returns contradiction")
			}

			vs, success := Solve(v)
			if !success || !IsSolved(vs) {
				t.Errorf("not solved board %v", board)
			}

			// Now solve again, with randomization
			vsr, success := Solve(v, SolveOptions{Randomize: true})
			if !success || !IsSolved(vsr) {
				t.Errorf("not solved randomized board %v", board)
			}
		}
	}
}

func TestSolveEmpty(t *testing.T) {
	vals := EmptyBoard()
	vres, solved := Solve(vals)
	if !solved {
		t.Errorf("want Solve(empty) to report success")
	}

	if !IsSolved(vres) {
		t.Errorf("want solved result board; got:\n%v", Display(vres))
	}

	// Try a few randomized solutions
	for i := 0; i < 10; i++ {
		vs, solved := Solve(vals, SolveOptions{Randomize: true})
		if !solved {
			t.Errorf("want Solve(empty) to report success")
		}
		if !IsSolved(vs) {
			t.Errorf("want solved result board; got:\n%v", Display(vs))
		}
	}
}

func BenchmarkParseBoardAssign(b *testing.B) {
	// Benchmark how long it takes to parse a board and run full constraint
	// propagation. We know that for easyboard1 it's fully solved with
	// constraint propagation after parsing.
	for i := 0; i < b.N; i++ {
		_, _ = ParseBoard(easyboard1, true)
	}
}

func BenchmarkSolveBoardHardlong(b *testing.B) {
	for i := 0; i < b.N; i++ {
		v, err := ParseBoard(hardlong, true)
		if err != nil {
			log.Fatal(err)
		}
		_, success := Solve(v)
		if !success {
			log.Fatal("not solved")
		}
	}
}

func BenchmarkSolveBoardHardlongRandomized(b *testing.B) {
	rand.Seed(time.Now().UnixNano())
	for i := 0; i < b.N; i++ {
		v, err := ParseBoard(hardlong, true)
		if err != nil {
			log.Fatal(err)
		}
		_, success := Solve(v, SolveOptions{Randomize: true})
		if !success {
			log.Fatal("not solved")
		}
	}
}

func BenchmarkSolveEmpty(b *testing.B) {
	// Benchmark how long it takes to "solve" an empty board.
	empty := EmptyBoard()
	for i := 0; i < b.N; i++ {
		_, _ = Solve(empty)
	}
}

func BenchmarkApplyTwinsStrategy(b *testing.B) {
	for i := 0; i < b.N; i++ {
		v := EmptyBoard()
		d38 := Digits(0).Add(3).Add(8)
		v[30] = d38
		v[31] = d38

		ApplyTwinsStrategy(v)
	}
}

func BenchmarkSolveEmptyRandomized(b *testing.B) {
	// Benchmark how long it takes to "solve" an empty board,
	// with randomization. Each solution will be different.
	rand.Seed(time.Now().UnixNano())
	empty := EmptyBoard()
	for i := 0; i < b.N; i++ {
		_, _ = Solve(empty)
	}
}


text:
golang.org/x/exp v0.0.0-20220613132600-b0d781184e0d h1:vtUKgx8dahOomfFzLREU8nSv25YHnTgLBn4rDnWZdU0=
golang.org/x/exp v0.0.0-20220613132600-b0d781184e0d/go.mod h1:Kr81I6Kryrl9sr8s2FK3vxD90NdsKWRuOIl2O4CvYbA=


text:
package sudoku

import (
	"fmt"
	"io"
	"math/rand"
	"strings"

	"github.com/eliben/go-sudoku/svg"
	"golang.org/x/exp/slices"
)

// Index represents a square on the Sudoku board; it's a number in the inclusive
// range [0, 80] that stands for row*9+col.
//
// These are the squares designated by an Index:
//
//  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 | 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 64 65 | 66 67 68 | 69 70 71
// 72 73 74 | 75 76 77 | 78 79 80
type Index = int

// Unit is a list of square indices that belong to the same Sudoku
// unit - a row, column or 3x3 block which should contain unique digits.
// There are overlaps between many units - e.g. a single square will be a member
// of a row unit, a column unit and a 3x3 block unit.
type Unit = []Index

// Values represents a Sudoku board in a format that's usable for solving.
// An element at index [i] in Values represents Sudoku square i (see the
// documentation of the Index type), and contains a set of all candidate
// digits for this square.
type Values []Digits

// unitlist is the list of all units that exist on the board.
var unitlist []Unit

// units maps an index to a list of units that contain that square.
// The mapping is a slice, i.e. units[i] is a list of all the units
// that contain the square with index i.
var units [][]Unit

// peers maps an index to a list of unique peers - other indices that share
// some unit with this index (it won't contain the index itself).
var peers [][]Index

func init() {
	index := func(row, col int) Index {
		return row*9 + col
	}

	// row units
	for row := 0; row < 9; row++ {
		var rowUnit []Index
		for col := 0; col < 9; col++ {
			rowUnit = append(rowUnit, index(row, col))
		}
		unitlist = append(unitlist, rowUnit)
	}

	// column units
	for col := 0; col < 9; col++ {
		var colUnit []Index
		for row := 0; row < 9; row++ {
			colUnit = append(colUnit, index(row, col))
		}
		unitlist = append(unitlist, colUnit)
	}

	// 3x3 block units
	for blockRow := 0; blockRow < 3; blockRow++ {
		for blockCol := 0; blockCol < 3; blockCol++ {
			var blockUnit []Index

			for row := 0; row < 3; row++ {
				for col := 0; col < 3; col++ {
					blockUnit = append(blockUnit, index(blockRow*3+row, blockCol*3+col))
				}
			}
			unitlist = append(unitlist, blockUnit)
		}
	}

	// For each index i, units[i] is a list of all units that contain i.
	units = make([][]Unit, 81)
	for i := 0; i < 81; i++ {
		for _, unit := range unitlist {
			if slices.Index(unit, i) >= 0 {
				units[i] = append(units[i], slices.Clone(unit))
			}
		}
	}

	// For each index i, peers[i] is a list of unique indices that share some
	// unit with i.
	peers = make([][]Index, 81)
	for i := 0; i < 81; i++ {
		for _, unit := range units[i] {
			for _, candidate := range unit {
				// This uses linear search to ensure uniqueness, but this calculation is
				// only done once at solver creation so we don't particularly care about
				// its speed.
				if candidate != i && slices.Index(peers[i], candidate) < 0 {
					peers[i] = append(peers[i], candidate)
				}
			}
		}
	}
}

// ParseBoard parses a Sudoku board given in textual representation, and returns
// it as Values. The textual representation is as described in
// http://norvig.com/sudoku.html: a string with a sequence of 81 runes in the
// set [0123456789.], where 0 or . mean "unassigned". All other runes in the
// string are ignored, making the input format flexible (it can be just 81
// runes on a single line, or it can be a nice multi-line |-separated board
// representation - since all runes outside the aforementioned set are ignored).
//
// If runElimination is false, the board is returned immediately after parsing.
// If runElimination is true, ParseBoard will invoke EliminateAll on the board
// and return the result. This is recommended when the board is then passed to
// a solver.
// It returns an error if there was an issue parsing the board, or if the board
// isn't a valid Sudoku board (e.g. contradictions exist).
func ParseBoard(str string, runElimination bool) (Values, error) {
	var dgs []uint16

	// Iterate and grab only the supported runes; ignore all others.
	for _, r := range str {
		if r >= '0' && r <= '9' {
			dgs = append(dgs, uint16(r)-uint16('0'))
		} else if r == '.' {
			dgs = append(dgs, 0)
		}
	}

	if len(dgs) != 81 {
		return nil, fmt.Errorf("got only %v digits in board, want 81", len(dgs))
	}

	// Start with an empty board.
	values := EmptyBoard()

	// Assign square digits based on the parsed board. Note that this runs
	// constraint propagation and may discover contradictions.
	for sq, d := range dgs {
		if d != 0 {
			values[sq] = SingleDigitSet(d)
		}
	}

	if runElimination && !EliminateAll(values) {
		return nil, fmt.Errorf("contradiction when eliminating board")
	}

	return values, nil
}

// EliminateAll runs elimination on all assigned squares in values. It applies
// first-order Sudoku heuristics on the entire board. Returns true if the
// elimination is successful, and false if the board has a contradiction.
func EliminateAll(values Values) bool {
	for sq, d := range values {
		if d.Size() == 1 {
			// Because of how eliminate() works, we prepare for it by remembering
			// which digit this square has assigned, setting the square to the full
			// set of digits and then calling eliminate on all digits except the
			// assigned one.
			digit := d.SingleMemberDigit()
			values[sq] = FullDigitsSet()
			for dn := uint16(1); dn <= 9; dn++ {
				if dn != digit {
					if !eliminate(values, sq, dn) {
						return false
					}
				}
			}
		}
	}
	return true
}

// CountHints counts the total number of hints on the board.
func CountHints(values Values) int {
	hintcount := 0
	for _, d := range values {
		if d.Size() == 1 {
			hintcount++
		}
	}
	return hintcount
}

// assign attempts to assign digit to values[square], propagating
// constraints from the assignment. values is modified.
// It returns true if the assignment succeeded, and false if the assignment
// fails resulting in an invalid Sudoku board.
func assign(values Values, square Index, digit uint16) bool {
	if EnableStats {
		Stats.NumAssigns++
	}

	for d := uint16(1); d <= 9; d++ {
		// For each d 1..9 that's != digit, if d is set in
		// values[square], try to eliminate it.
		if values[square].IsMember(d) && d != digit {
			if !eliminate(values, square, d) {
				return false
			}
		}
	}
	return true
}

// eliminate removes digit from the candidates in values[square], propagating
// constraints. values is modified.
// It returns false if this results in an invalid Sudoku board; otherwise
// returns true.
func eliminate(values Values, square Index, digit uint16) bool {
	if !values[square].IsMember(digit) {
		// Already eliminated
		return true
	}

	// Remove digit from the candidates in square.
	values[square] = values[square].Remove(digit)

	switch values[square].Size() {
	case 0:
		// No remaining options for square -- this is a contradiction.
		return false
	case 1:
		// A single digit candidate remaining in the square -- this creates a new
		// constraint. Eliminate this digit from all peer squares.
		remaining := values[square].SingleMemberDigit()
		for _, peer := range peers[square] {
			if !eliminate(values, peer, remaining) {
				return false
			}
		}
	}

	// Since digit was eliminated from square, it's possible that we'll find a
	// position for this digit in one of the units the square belongs to.
UnitLoop:
	for _, unit := range units[square] {
		// Looking for a single square in this unit that has 'digit' as one of its
		// candidates. sqd marks the square, or -1 if no such square was found.
		sqd := -1
		for _, sq := range unit {
			if values[sq].IsMember(digit) {
				if sqd == -1 {
					sqd = sq
				} else {
					// More than one square has 'digit' as a candidate, so we won't be
					// able to simplify things.
					continue UnitLoop
				}
			}
		}
		if sqd == -1 {
			// Contradiction: no places left in this unit for 'digit'
			return false
		}

		// There's only a single place left in the unit for 'digit' to go, so
		// assign it.
		if !assign(values, sqd, digit) {
			return false
		}
	}

	return true
}

// Display returns a visual representation of values, with all the digit
// candidates as a string in each cell.
func Display(values Values) string {
	// Find maximum length of one square.
	var maxlen int = 0
	for _, d := range values {
		if d.Size() > maxlen {
			maxlen = d.Size()
		}
	}
	width := maxlen + 1

	line := strings.Join([]string{
		strings.Repeat("-", width*3),
		strings.Repeat("-", width*3),
		strings.Repeat("-", width*3)}, "+")

	var sb strings.Builder
	for sq, d := range values {
		fmt.Fprintf(&sb, "%[1]*s", -width, fmt.Sprintf("%[1]*s", (width+d.Size())/2, d))
		if sq%9 == 2 || sq%9 == 5 {
			sb.WriteString("|")
		}
		if sq%9 == 8 {
			sb.WriteRune('\n')
		}
		if sq == 26 || sq == 53 {
			sb.WriteString(line + "\n")
		}
	}
	return sb.String()
}

// DisplayAsInput returns a 2D Sudoku input board corresponding to values.
// It treats solved squares (with one candidate) as hints that are filled into
// the board, and unsolved squares (with more than one candidate) as empty.
func DisplayAsInput(values Values) string {
	line := strings.Join([]string{
		strings.Repeat("-", 6),
		strings.Repeat("-", 6),
		strings.Repeat("-", 6)}, "+")

	var sb strings.Builder
	for sq, d := range values {
		ds := d.String()
		if d.Size() > 1 {
			ds = "."
		}
		fmt.Fprintf(&sb, "%s ", ds)
		if sq%9 == 2 || sq%9 == 5 {
			sb.WriteString("|")
		}
		if sq%9 == 8 {
			sb.WriteRune('\n')
		}
		if sq == 26 || sq == 53 {
			sb.WriteString(line + "\n")
		}
	}
	return sb.String()
}

// DisplayAsSVG write the board's visual representation in SVG format into w.
// The difficulty is emitted too.
func DisplayAsSVG(w io.Writer, values Values, difficulty float64) {
	startX := 50
	startY := 50
	width := 800
	height := 900
	cellsize := 80
	canvas := svg.New(w, width, height)

	for sq, d := range values {
		col := sq % 9
		x := startX + col*cellsize

		row := sq / 9
		y := startY + row*cellsize

		canvas.Rect(x, y, cellsize, cellsize, "stroke:black; stroke-width:2; fill:white")
		if d.Size() == 1 {
			canvas.Text(x+cellsize/2, y+cellsize/2, d.String(), "text-anchor:middle; dominant-baseline:middle; font-family:Helvetica; font-size:32px; fill:black")
		}
	}

	// Wider squares around 3x3 blocks
	for br := 0; br < 3; br++ {
		for bc := 0; bc < 3; bc++ {
			canvas.Rect(startX+bc*cellsize*3, startY+br*cellsize*3, cellsize*3, cellsize*3, "stroke:black; stroke-width:5; fill-opacity:0.0")
		}
	}

	difficultyText := fmt.Sprintf("Difficulty: %.2f out of 5", difficulty)
	canvas.Text(startX, startY+9*cellsize+cellsize/2, difficultyText, "font-family:Helvetica; font-size:16px; fill:black")

	canvas.End()
}

// EmptyBoard creates an "empty" Sudoku board, where each square can potentially
// contain any digit.
func EmptyBoard() Values {
	vals := make(Values, 81)
	for sq := range vals {
		vals[sq] = FullDigitsSet()
	}
	return vals
}

// IsSolved checks whether values is a properly solved Sudoku board, with all
// the constraints satisfied.
func IsSolved(values Values) bool {
	for _, unit := range unitlist {
		var dset Digits
		for _, sq := range unit {
			// Some squares have more than a single candidate? Not solved.
			if values[sq].Size() != 1 {
				return false
			}
			dset = dset.Add(values[sq].SingleMemberDigit())
		}
		// Not all digits covered by this unit? Not solved.
		if dset != FullDigitsSet() {
			return false
		}
	}
	return true
}

// findSquareWithFewestCandidates finds a square in values with more than one
// digit candidate, but the smallest number of such candidates.
func findSquareWithFewestCandidates(values Values) Index {
	var squareToTry Index = -1
	var minSize int = 10
	for sq, d := range values {
		if d.Size() > 1 && d.Size() < minSize {
			minSize = d.Size()
			squareToTry = sq
		}
	}
	return squareToTry
}

// SolveOptions is a container of options for the Solve function.
type SolveOptions struct {
	// Randomize tells the solver to randomly shuffle its digit selection when
	// attempting to guess a value for a square. For actual randomness, the
	// rand package's default randomness source should be properly seeded before
	// invoking Solve.
	Randomize bool
}

// Solve runs a backtracking search to solve the board given in values.
// It returns true and the solved values if the search succeeded and we ended up
// with a board with only a single candidate per square; otherwise, it returns
// false. The input values is not modified.
// The solution process can be configured by providing SolveOptions.
func Solve(values Values, options ...SolveOptions) (Values, bool) {
	if len(options) > 1 {
		panic("Solve cannot accept more than a single SolveOptions")
	}

	squareToTry := findSquareWithFewestCandidates(values)

	// If we didn't find any square with more than one candidate, the board is
	// solved!
	if squareToTry == -1 {
		return values, true
	}

	if EnableStats {
		Stats.NumSearches++
	}

	var candidates = []uint16{1, 2, 3, 4, 5, 6, 7, 8, 9}
	if len(options) > 0 && options[0].Randomize {
		rand.Shuffle(len(candidates), func(i, j int) {
			candidates[i], candidates[j] = candidates[j], candidates[i]
		})
	}

	for _, d := range candidates {
		// Try to assign sq with each one of its candidate digits. If this results
		// in a successful Solve() - we've solved the board!
		if values[squareToTry].IsMember(d) {
			vcopy := slices.Clone(values)
			if assign(vcopy, squareToTry, d) {
				if vresult, solved := Solve(vcopy, options...); solved {
					return vresult, true
				}
			}
		}
	}
	return values, false
}

// SolveAll finds all solutions to the given board and returns them. If no
// solutions were found, an empty list is returned. max can specify the
// (approximate) maximal number of solutions to find; a value <= 0 means "all of
// them". Often more solutions than max will be returned, but not a lot more
// (maybe 2-3x as many). values is not modified.
// Warning: this function can take a LONG time to run for boards with multiple
// solutions, and it can consume enormous amounts of memory because it has to
// remember each solution it finds. For some boards it will run forever (e.g.
// finding all solutions on an empty board). If in doubt, use the max parameter
// to restrict the number.
func SolveAll(values Values, max int) []Values {
	squareToTry := findSquareWithFewestCandidates(values)

	// If we didn't find any square with more than one candidate, the board is
	// solved!
	if squareToTry == -1 {
		return []Values{values}
	}

	var allSolved []Values

	for d := uint16(1); d <= 9; d++ {
		// Try to assign sq with each one of its candidate digits. If this results
		// in a successful Solve() - we've solved the board!
		if values[squareToTry].IsMember(d) {
			vcopy := slices.Clone(values)
			if assign(vcopy, squareToTry, d) {
				if vsolved := SolveAll(vcopy, max); len(vsolved) > 0 {
					allSolved = append(allSolved, vsolved...)
					if max > 0 && len(allSolved) >= max {
						return allSolved
					}
				}
			}
		}
	}
	return allSolved
}

// EnableStats enables statistics collection during the processes of solving.
// When stats are enabled, solving will be slightly slower.
//
// Note: statistics collection is NOT SAFE FOR CONCURRENT ACCESS.
var EnableStats bool = false

type StatsCollector struct {
	NumSearches uint64
	NumAssigns  uint64
}

// Stats is the global variable for accessing statistics from this package.
// It's recommended to call Stats.Reset() before solving a board, and access
// the Stats fields after it's done.
var Stats StatsCollector

func (s *StatsCollector) Reset() {
	s.NumSearches = 0
	s.NumAssigns = 0
}

// WithStats helps run any block of code with stats enabled.
func WithStats(f func()) {
	EnableStats = true
	defer func() {
		EnableStats = false
	}()
	Stats.Reset()

	f()
}

// ApplyTwinsStrategy applies the "naked twins" Sudoku strategy to the given
// board and updates it. It returns false if there was a contradiction
// discovered while applying the strategy.
// The "naked twins" strategy is to find two digits that are the only candidates
// in two squares in the same unit. This helps eliminate these digits from
// other squares in the same unit.
// For example, two squares in the same row cell may have 38 as their only
// candidates, which means that 3 and 8 must occupy these squares (though we
// don't know which goes where), and that no other square in the unit may have
// either 3 or 8.
func ApplyTwinsStrategy(values Values) bool {
	// The strategy is repeated to a "fixed point" where further runs don't end
	// up changing the board in any way.
RepeatStrategy:
	for {
		for _, unit := range unitlist {
			// dcount will map Digits->count, counting how many times a certain
			// combination of digit candidates appears in this unit.
			dcount := make(map[Digits]int)
			for _, sq := range unit {
				dcount[values[sq]]++
			}

			var removed bool
			for d, count := range dcount {
				if d.Size() == 2 && count == 2 {
					// Found it! Now go over all the squares in this unit other than the
					// ones with these exact two candidates, and eliminate both candidates
					// from them.
					d1, d2 := d.twoMemberDigits()
					for _, sq := range unit {
						if values[sq].Size() >= 2 && values[sq] != d {
							if values[sq].IsMember(d1) {
								if !eliminate(values, sq, d1) {
									return false
								}
								removed = true
							}
							if values[sq].IsMember(d2) {
								if !eliminate(values, sq, d2) {
									return false
								}
								removed = true
							}
						}
					}
				}
				// We've eliminated some values from the board, so it's not safe to
				// proceed with this loop since dcount may hold stale values. We
				// instead repeat the strategy on the whole board.
				if removed {
					continue RepeatStrategy
				}
			}
		}

		// We went over all units and didn't remove anything in this round, so
		// no point repeating.
		return true
	}
}


text:
package sudoku

import (
	"math/rand"
	"testing"
	"time"
)

func TestGenerate(t *testing.T) {
	board := Generate(30)
	vs := SolveAll(board, -1)
	if len(vs) != 1 {
		t.Errorf("got %v solutions, want 1", len(vs))
	}

	if !IsSolved(vs[0]) {
		t.Errorf("got unsolved board")
	}
}

func TestGenerateSymmetrical(t *testing.T) {
	rand.Seed(time.Now().UnixNano())
	//for {
	board := GenerateSymmetrical(30)
	vs := SolveAll(board, -1)
	if len(vs) != 1 {
		t.Errorf("got %v solutions, want 1", len(vs))
	}

	if !IsSolved(vs[0]) {
		t.Errorf("got unsolved board")
	}

	// Check symmetry
	for sq := 0; sq < 41; sq++ {
		if board[sq].Size() != board[80-sq].Size() {
			t.Errorf("squares %v != %v on board, expected symmetry", sq, 80-sq)
		}
	}
}


text:
# go-sudoku

This is a toolkit for solving and generating [Sudoku
puzzles](https://en.wikipedia.org/wiki/Sudoku) in Go.

## What's inside

All the code is in a single top-level package: `sudoku`. It's broadly separated
into three parts:

* `sudoku.go`: board representation and functions for parsing boards from
  strings, emitting boards back to output and solving Sudoku puzzles. The
  basic solver uses constraint propagation and recursive search and is based on
  [Peter Norvig's old post](https://norvig.com/sudoku.html), although the Go
  code is about 100x faster than Norvig's Python (faster compiled language but
  also an optimized board representation).

  Contains additional functionality like finding _all_ the solutions of a given
  puzzle and not just a single solution.

* `generator.go`: generate valid Sudoku puzzles that have a single solution.
  The algorithm is based on a mish-mash of information found online and tweaked
  by me. Contains additional functionality like generating _symmetrical_
  Sudoku boards.

  Note: generating hard-to-solve boards with a single solution is fairly
  difficult. The best way to do this in practice seems to be to generate a
  _large_ number of boards, sifting through them to find the hardest ones.
  These boards can then be transformed in a myriad ways to retain the same
  difficulty but look and feel very different (through swapping rows and
  columns, rotations, and permuting the existing hint digits). Therefore,
  a single genuienly hard board can be replayed in many different ways.

* `difficulty.go`: code to evaluate the difficulty of a given Sudoku puzzle;
  the approach was partially inspired by the paper "Sudoku Puzzles Generating:
  from Easy to Evil" by Xiang-Sun ZHANG's research group.

The `cmd` directory has two command-line tools: `generator` and `solver` that
demonstrate the use of the package.

## Testing

Some tests take a while to run, so they are excluded if the `-short` testing
flag is provided:

    $ go test -v -short ./...

## Generating printable boards

`go-sudoku` includes some rudimentary functionality to emit a Sudoku board into
a printable SVG format, like this:

![SVG board sample](/doc/genb.svg)

You can invoke the `cmd/generator` command with the `-svgout` flag to see this
in action, or use the web interface.

## Web interface

This repository includes a web interface for generating Sudoku puzzles, by
compiling the Go code to WebAssembly. To run it locally:

		$ cd cmd/wasm
		$ make build
		$ make serve

This will run a local webserver; open http://localhost:8899 in your browser
to generate puzzles!

The repository also has a GitHub actions setup to automatically deploy the web
interface to https://eliben.github.io/go-sudoku on each commit.


text:
# Simple workflow for deploying static content to GitHub Pages
name: Deploy static content to Pages

on:
  # Runs on pushes targeting the default branch
  push:
    branches: ["main"]

  # Allows you to run this workflow manually from the Actions tab
  workflow_dispatch:

# Sets permissions of the GITHUB_TOKEN to allow deployment to GitHub Pages
permissions:
  contents: read
  pages: write
  id-token: write

# Allow one concurrent deployment
concurrency:
  group: "pages"
  cancel-in-progress: true

jobs:
  # Single deploy job since we're just deploying
  deploy:
    environment:
      name: github-pages
      url: ${{ steps.deployment.outputs.page_url }}
    runs-on: ubuntu-latest
    steps:
      - name: Checkout
        uses: actions/checkout@v3

      - name: Set up Go
        uses: actions/setup-go@v3
        with:
          go-version: 1.18

      - name: Setup Pages
        uses: actions/configure-pages@v2

      - name: Build wasm
        run: |
          cd cmd/wasm
          make build

      - name: Upload artifact
        uses: actions/upload-pages-artifact@v1
        with:
          # Upload
          path: 'cmd/wasm/assets'

      - name: Deploy to GitHub Pages
        id: deployment
        uses: actions/deploy-pages@v1


text:
name: Run Go tests

on:
  push:
    branches: [ "main" ]
  pull_request:
    branches: [ "main" ]

jobs:

  build:
    runs-on: ubuntu-latest
    steps:
    - uses: actions/checkout@v3

    - name: Set up Go
      uses: actions/setup-go@v3
      with:
        go-version: 1.18

    - name: Test
      # Running ./... won't work since cmd/wasm needs js tag settings.
      run: go test -v . ./svg


text:
# These are taken from Norvig's collection of hard puzzles found online; a
# combination of http://norvig.com/top95.txt and https://norvig.com/hardest.txt
4.....8.5.3..........7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......
52...6.........7.13...........4..8..6......5...........418.........3..2...87.....
6.....8.3.4.7.................5.4.7.3..2.....1.6.......2.....5.....8.6......1....
48.3............71.2.......7.5....6....2..8.............1.76...3.....4......5....
....14....3....2...7..........9...3.6.1.............8.2.....1.4....5.6.....7.8...
......52..8.4......3...9...5.1...6..2..7........3.....6...1..........7.4.......3.
6.2.5.........3.4..........43...8....1....2........7..5..27...........81...6.....
.524.........7.1..............8.2...3.....6...9.5.....1.6.3...........897........
6.2.5.........4.3..........43...8....1....2........7..5..27...........81...6.....
.923.........8.1...........1.7.4...........658.........6.5.2...4.....7.....9.....
6..3.2....5.....1..........7.26............543.........8.15........4.2........7..
.6.5.1.9.1...9..539....7....4.8...7.......5.8.817.5.3.....5.2............76..8...
..5...987.4..5...1..7......2...48....9.1.....6..2.....3..6..2.......9.7.......5..
3.6.7...........518.........1.4.5...7.....6.....2......2.....4.....8.3.....5.....
1.....3.8.7.4..............2.3.1...........958.........5.6...7.....8.2...4.......
6..3.2....4.....1..........7.26............543.........8.15........4.2........7..
....3..9....2....1.5.9..............1.2.8.4.6.8.5...2..75......4.1..6..3.....4.6.
45.....3....8.1....9...........5..9.2..7.....8.........1..4..........7.2...6..8..
.237....68...6.59.9.....7......4.97.3.7.96..2.........5..47.........2....8.......
..84...3....3.....9....157479...8........7..514.....2...9.6...2.5....4......9..56
.98.1....2......6.............3.2.5..84.........6.........4.8.93..5...........1..
..247..58..............1.4.....2...9528.9.4....9...1.........3.3....75..685..2...
4.....8.5.3..........7......2.....6.....5.4......1.......6.3.7.5..2.....1.9......
.2.3......63.....58.......15....9.3....7........1....8.879..26......6.7...6..7..4
1.....7.9.4...72..8.........7..1..6.3.......5.6..4..2.........8..53...7.7.2....46
4.....3.....8.2......7........1...8734.......6........5...6........1.4...82......
.......71.2.8........4.3...7...6..5....2..3..9........6...7.....8....4......5....
6..3.2....4.....8..........7.26............543.........8.15........8.2........7..
.47.8...1............6..7..6....357......5....1..6....28..4.....9.1...4.....2.69.
......8.17..2........5.6......7...5..1....3...8.......5......2..4..8....6...3....
38.6.......9.......2..3.51......5....3..1..6....4......17.5..8.......9.......7.32
...5...........5.697.....2...48.2...25.1...3..8..3.........4.7..13.5..9..2...31..
.2.......3.5.62..9.68...3...5..........64.8.2..47..9....3.....1.....6...17.43....
.8..4....3......1........2...5...4.69..1..8..2...........3.9....6....5.....2.....
..8.9.1...6.5...2......6....3.1.7.5.........9..4...3...5....2...7...3.8.2..7....4
4.....5.8.3..........7......2.....6.....5.8......1.......6.3.7.5..2.....1.8......
1.....3.8.6.4..............2.3.1...........958.........5.6...7.....8.2...4.......
1....6.8..64..........4...7....9.6...7.4..5..5...7.1...5....32.3....8...4........
249.6...3.3....2..8.......5.....6......2......1..4.82..9.5..7....4.....1.7...3...
...8....9.873...4.6..7.......85..97...........43..75.......3....3...145.4....2..1
...5.1....9....8...6.......4.1..........7..9........3.8.....1.5...2..4.....36....
......8.16..2........7.5......6...2..1....3...8.......2......7..3..8....5...4....
.476...5.8.3.....2.....9......8.5..6...1.....6.24......78...51...6....4..9...4..7
.....7.95.....1...86..2.....2..73..85......6...3..49..3.5...41724................
.4.5.....8...9..3..76.2.....146..........9..7.....36....1..4.5..6......3..71..2..
.834.........7..5...........4.1.8..........27...3.....2.6.5....5.....8........1..
..9.....3.....9...7.....5.6..65..4.....3......28......3..75.6..6...........12.3.8
.26.39......6....19.....7.......4..9.5....2....85.....3..2..9..4....762.........4
2.3.8....8..7...........1...6.5.7...4......3....1............82.5....6...1.......
6..3.2....1.....5..........7.26............843.........8.15........8.2........7..
1.....9...64..1.7..7..4.......3.....3.89..5....7....2.....6.7.9.....4.1....129.3.
.........9......84.623...5....6...453...1...6...9...7....1.....4.5..2....3.8....9
.2....5938..5..46.94..6...8..2.3.....6..8.73.7..2.........4.38..7....6..........5
9.4..5...25.6..1..31......8.7...9...4..26......147....7.......2...3..8.6.4.....9.
...52.....9...3..4......7...1.....4..8..453..6...1...87.2........8....32.4..8..1.
53..2.9...24.3..5...9..........1.827...7.........981.............64....91.2.5.43.
1....786...7..8.1.8..2....9........24...1......9..5...6.8..........5.9.......93.4
....5...11......7..6.....8......4.....9.1.3.....596.2..8..62..7..7......3.5.7.2..
.47.2....8....1....3....9.2.....5...6..81..5.....4.....7....3.4...9...1.4..27.8..
......94.....9...53....5.7..8.4..1..463...........7.8.8..7.....7......28.5.26....
.2......6....41.....78....1......7....37.....6..412....1..74..5..8.5..7......39..
1.....3.8.6.4..............2.3.1...........758.........7.5...6.....8.2...4.......
2....1.9..1..3.7..9..8...2.......85..6.4.........7...3.2.3...6....5.....1.9...2.5
..7..8.....6.2.3...3......9.1..5..6.....1.....7.9....2........4.83..4...26....51.
...36....85.......9.4..8........68.........17..9..45...1.5...6.4....9..2.....3...
34.6.......7.......2..8.57......5....7..1..2....4......36.2..1.......9.......7.82
......4.18..2........6.7......8...6..4....3...1.......6......2..5..1....7...3....
.4..5..67...1...4....2.....1..8..3........2...6...........4..5.3.....8..2........
.......4...2..4..1.7..5..9...3..7....4..6....6..1..8...2....1..85.9...6.....8...3
8..7....4.5....6............3.97...8....43..5....2.9....6......2...6...7.71..83.2
.8...4.5....7..3............1..85...6.....2......4....3.26............417........
....7..8...6...5...2...3.61.1...7..2..8..534.2..9.......2......58...6.3.4...1....
......8.16..2........7.5......6...2..1....3...8.......2......7..4..8....5...3....
.2..........6....3.74.8.........3..2.8..4..1.6..5.........1.78.5....9..........4.
.52..68.......7.2.......6....48..9..2..41......1.....8..61..38.....9...63..6..1.9
....1.78.5....9..........4..2..........6....3.74.8.........3..2.8..4..1.6..5.....
1.......3.6.3..7...7...5..121.7...9...7........8.1..2....8.64....9.2..6....4.....
4...7.1....19.46.5.....1......7....2..2.3....847..6....14...8.6.2....3..6...9....
......8.17..2........5.6......7...5..1....3...8.......5......2..3..8....6...4....
963......1....8......2.5....4.8......1....7......3..257......3...9.2.4.7......9..
15.3......7..4.2....4.72.....8.........9..1.8.1..8.79......38...........6....7423
..........5724...98....947...9..3...5..9..12...3.1.9...6....25....56.....7......6
....75....1..2.....4...3...5.....3.2...8...1.......6.....1..48.2........7........
6.....7.3.4.8.................5.4.8.7..2.....1.3.......2.....5.....7.9......1....
....6...4..6.3....1..4..5.77.....8.5...8.....6.8....9...2.9....4....32....97..1..
.32.....58..3.....9.428...1...4...39...6...5.....1.....2...67.8.....4....95....6.
...5.3.......6.7..5.8....1636..2.......4.1.......3...567....2.8..4.7.......2..5..
.5.3.7.4.1.........3.......5.8.3.61....8..5.9.6..1........4...6...6927....2...9..
..5..8..18......9.......78....4.....64....9......53..2.6.........138..5....9.714.
..........72.6.1....51...82.8...13..4.........37.9..1.....238..5.4..9.........79.
...658.....4......12............96.7...3..5....2.8...3..19..8..3.6.....4....473..
.2.3.......6..8.9.83.5........2...8.7.9..5........6..4.......1...1...4.22..7..8.9
.5..9....1.....6.....3.8.....8.4...9514.......3....2..........4.8...6..77..15..6.
.....2.......7...17..3...9.8..7......2.89.6...13..6....9..5.824.....891..........
3...8.......7....51..............36...2..4....7...........6.13..452...........8..
85...24..72......9..4.........1.7..23.5...9...4...........8..7..17..........36.4.
..53.....8......2..7..1.5..4....53...1..7...6..32...8..6.5....9..4....3......97..
12..4......5.69.1...9...5.........7.7...52.9..3......2.9.6...5.4..9..8.1..3...9.4
...57..3.1......2.7...234......8...4..7..4...49....6.5.42...3.....7..9....18.....
7..1523........92....3.....1....47.8.......6............9...5.6.4.9.7...8....6.1.
1....7.9..3..2...8..96..5....53..9...1..8...26....4...3......1..4......7..7...3..
1...34.8....8..5....4.6..21.18......3..1.2..6......81.52..7.9....6..9....9.64...2
...92......68.3...19..7...623..4.1....1...7....8.3..297...8..91...5.72......64...
.6.5.4.3.1...9...8.........9...5...6.4.6.2.7.7...4...5.........4...8...1.5.2.3.4.
7.....4...2..7..8...3..8.799..5..3...6..2..9...1.97..6...3..9...3..4..6...9..1.35
....7..2.8.......6.1.2.5...9.54....8.........3....85.1...3.2.8.4.......9.7..6....


text:
# The hardest board Norvig generated in a million random puzzles. Note: it has
# more than one solution.
.....6....59.....82....8....45........3........6..3.54...325..6..................


text:
package svg

import (
	"encoding/xml"
	"fmt"
	"io"
	"text/template"
)

type Canvas struct {
	writer io.Writer
}

var headerText = `<?xml version="1.0"?>
<svg width="{{.Width}}" height="{{.Height}}"
     xmlns="http://www.w3.org/2000/svg"
     xmlns:xlink="http://www.w3.org/1999/xlink">
`

var headerTemplate = template.Must(template.New("header").Parse(headerText))

func New(writer io.Writer, width, height int) *Canvas {
	c := &Canvas{writer: writer}

	headerTemplate.Execute(writer, struct {
		Width, Height int
	}{Width: width, Height: height})

	return c
}

func (c *Canvas) End() {
	fmt.Fprintf(c.writer, "</svg>\n")
}

func (c *Canvas) Rect(x, y, width, height int, style string) {
	fmt.Fprintf(c.writer, `<rect x="%v" y="%v" width="%v" height="%v"`, x, y, width, height)
	if len(style) > 0 {
		fmt.Fprintf(c.writer, ` style="%s"`, style)
	}
	fmt.Fprintf(c.writer, "/>\n")
}

func (c *Canvas) Text(x, y int, text string, style string) {
	fmt.Fprintf(c.writer, `<text x="%v" y="%v"`, x, y)
	if len(style) > 0 {
		fmt.Fprintf(c.writer, ` style="%s"`, style)
	}
	fmt.Fprintf(c.writer, ">")
	xml.Escape(c.writer, []byte(text))
	fmt.Fprintf(c.writer, "</text>\n")
}


text:
package svg

import (
	"bytes"
	"strings"
	"testing"
)

// Sanity testing.
func TestSvgWrite(t *testing.T) {
	var buf bytes.Buffer

	width := 120
	height := 233
	x := 77
	y := 88

	canvas := New(&buf, width, height)
	canvas.Rect(x, y, 100, 200, "my style")
	canvas.Text(x+10, y+1, "hello", "")
	canvas.End()

	result := buf.String()

	want := `<?xml version="1.0"?>
<svg width="120" height="233"
     xmlns="http://www.w3.org/2000/svg"
     xmlns:xlink="http://www.w3.org/1999/xlink">
<rect x="77" y="88" width="100" height="200" style="my style"/>
<text x="87" y="89">hello</text>
</svg>`

	if strings.TrimSpace(result) != strings.TrimSpace(want) {
		t.Errorf("got:\n %s\n\nwant:\n %v\n", result, want)
	}
}


text:
package main

import (
	"bufio"
	"flag"
	"fmt"
	"log"
	"math/rand"
	"os"
	"strings"
	"time"

	"github.com/eliben/go-sudoku"
)

var statsFlag = flag.Bool("stats", false, "enable stats for solving")
var randomizeFlag = flag.Bool("randomize", false, "randomize solving order")
var actionFlag = flag.String("action", "solve", "action to perform: solve, count")

func main() {
	flag.Usage = func() {
		out := flag.CommandLine.Output()
		fmt.Fprintln(out, "usage: solver [options] <input in stdin>")
		fmt.Fprintln(out, "Options:")
		flag.PrintDefaults()
	}
	flag.Parse()

	switch *actionFlag {
	case "solve":
		solveAndReport()
	case "count":
		countHints()
	default:
		flag.Usage()
		log.Fatal("Please select one of the supported actions.")
	}
}

func solveAndReport() {
	var totalDuration time.Duration = 0
	var maxDuration time.Duration = 0
	var totalSearches uint64 = 0
	var totalDifficulty float64
	var maxSearches uint64 = 0
	var numBoards int = 0
	var numSolved int = 0

	if *statsFlag {
		sudoku.EnableStats = true
	}

	if *randomizeFlag {
		rand.Seed(time.Now().UnixNano())
	}

	boards := getInputBoards()
	for _, board := range boards {
		numBoards++

		v, err := sudoku.ParseBoard(board, false)
		if err != nil {
			log.Fatal(err)
		}
		d, err := sudoku.EvaluateDifficulty(v)
		if err != nil {
			log.Fatal(err)
		}
		totalDifficulty += d

		tStart := time.Now()
		sudoku.EliminateAll(v)
		v, _ = sudoku.Solve(v, sudoku.SolveOptions{Randomize: *randomizeFlag})
		if err != nil {
			log.Fatal(err)
		}
		tElapsed := time.Now().Sub(tStart)

		totalDuration += tElapsed
		if tElapsed > maxDuration {
			maxDuration = tElapsed
		}

		if sudoku.IsSolved(v) {
			numSolved++
		}

		if *statsFlag {
			totalSearches += sudoku.Stats.NumSearches
			if sudoku.Stats.NumSearches > maxSearches {
				maxSearches = sudoku.Stats.NumSearches
			}
			sudoku.Stats.Reset()
		}
	}

	fmt.Printf("Solved %v/%v boards\n", numSolved, numBoards)
	fmt.Printf("Average difficulty: %.2v\n", totalDifficulty/float64(numBoards))
	fmt.Printf("Duration average=%-15v max=%v\n", totalDuration/time.Duration(numBoards), maxDuration)
	if *statsFlag {
		fmt.Printf("Searches average=%-15.2f max=%v\n", float64(totalSearches)/float64(numBoards), maxSearches)
	}
}

func countHints() {
	boards := getInputBoards()
	for _, board := range boards {
		fmt.Println("board:", board)
		v, err := sudoku.ParseBoard(board, false)
		if err != nil {
			log.Fatal(err)
		}
		fmt.Println("|")

		initialNumHints := sudoku.CountHints(v)
		fmt.Printf("\\ initial num hints:           %v\n", initialNumHints)

		sudoku.EliminateAll(v)
		afterElimNumHints := sudoku.CountHints(v)
		fmt.Printf("  num hints after elimination: %v\n", afterElimNumHints)

		sudoku.ApplyTwinsStrategy(v)
		afterTwinsNumHints := sudoku.CountHints(v)
		fmt.Printf("  num hints after twins:       %v\n", afterTwinsNumHints)
		fmt.Println("")
	}
}

// getInputBoards reads input boards from stdin, ignores comments and empty
// lines and returns them.
func getInputBoards() []string {
	var boards []string

	scanner := bufio.NewScanner(os.Stdin)
	// Expect one board per line, ignoring whitespace and lines starting with '#'.
	for scanner.Scan() {
		board := strings.TrimSpace(scanner.Text())
		if len(board) == 0 || strings.HasPrefix(board, "#") {
			continue
		}
		boards = append(boards, board)
	}

	if err := scanner.Err(); err != nil {
		log.Fatal(err)
	}
	return boards
}


text:
package main

import (
	"bytes"
	"fmt"
	"log"
	"math/rand"
	"syscall/js"
	"time"

	"github.com/eliben/go-sudoku"
)

func main() {
	rand.Seed(time.Now().UnixNano())
	fmt.Println("go-sudoku wasm")

	// Export the jsGenerateBoard function to JS.
	js.Global().Set("generateBoard", jsGenerateBoard)

	// For the Go code to be usable from JS, the main function has to run forever.
	<-make(chan bool)
}

// jsGenerateBoard wraps the functionality we need from this package, for use
// in the web interface. It creates a function that takes two parameters:
// an integer hint count, and a boolean "is symmetrical" flag. It returns
// the SVG generated for the board as a string.
var jsGenerateBoard = js.FuncOf(func(this js.Value, args []js.Value) interface{} {
	if len(args) != 2 {
		return fmt.Sprintf("got %v args, want 2", len(args))
	}
	hintCount := args[0].Int()
	symmetrical := args[1].Bool()

	var board sudoku.Values
	if symmetrical {
		board = sudoku.GenerateSymmetrical(hintCount)
	} else {
		board = sudoku.Generate(hintCount)
	}

	d, err := sudoku.EvaluateDifficulty(board)
	if err != nil {
		log.Fatal(err)
	}

	var buf bytes.Buffer
	sudoku.DisplayAsSVG(&buf, board, d)
	return buf.String()
})


text:
.PHONY: build serve

build:
	GOOS=js GOARCH=wasm go build -o assets/gosudoku.wasm main.go

serve:
	go run ./server