Language: C++

#include <stdio.h>
#include <iostream>
#include <cstdio>
#include <algorithm> 
#include <stdio.h>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <map>
#include <list>
#include <utility>
#include <set>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <numeric>
#include <bitset>
#include <stdio.h>
#include <string.h>
#include <assert.h>
#include <random>

using namespace std;

static const double EPS = 1e-9;
template<class T> bool INRANGE(T x, T a, T b) { return a <= x && x <= b; }
template<class T> void amin(T &a, T v) { if (a > v) a = v; }
template<class T> void amax(T &a, T v) { if (a < v) a = v; }
int ROUND(double x) { return (int)(x + 0.5); }
bool ISINT(double x) { return fabs(ROUND(x) - x) <= EPS; }
bool ISEQUAL(double x, double y) { return fabs(x - y) <= EPS * max(1.0, max(fabs(x), fabs(y))); }
double SQSUM(double x, double y) { return x * x + y * y; }
#define PI  (acos(-1))
#define ARRAY_NUM(a) (sizeof(a)/sizeof(a[0])) 
#define NG (-1)
#define BIG ((int)1e9+10)
#define BIGLL ((ll)4e18)
#define SZ(a) ((int)(a).size())
#define SQ(a) ((a)*(a))
typedef unsigned long long ull;
typedef long long ll;

using namespace std;
using Int = long long;
template<typename T1, typename T2> inline void chmin(T1 &a, T2 b) { if (a>b) a = b; }
template<typename T1, typename T2> inline void chmax(T1 &a, T2 b) { if (a<b) a = b; }

enum AlgorithmMode
{
	SA,
	HillClimb,
};

// 焼きなまし法の状態をIStateを実装する必要がある（ポリモーフィズムは使ってない。テンプレートのかわり）
class IState
{
public:
	virtual	void	saInit() = 0;
	virtual	double	saScore() = 0;			// ペナルティなども含んだ、焼きなまし用のスコア。
	virtual	double	saMove() = 0;
	virtual	void	saForward() = 0;
	virtual	void	saResume() = 0;
	virtual void	saSetProgressRatio(double progressRatio, double temp) = 0;
#ifdef SA_STATS
	virtual int 	saGetNeighborID() = 0;
#endif // SA_STATS
};

/*
SAState state(this);
state.mCommands = saParams->prevBestCommands;
SimulatedAnnealing<SAState> sa;

double          maxTemp = 2.0;
double          minTemp = 0.001;
double          tempExponential = 1;
int             numMultiStart = 1;

sa.init(state, SA, 5000, 1000, maxTemp, maxTemp, minTemp, tempExponential, numMultiStart, saParams->handleID);

#if 1
ll startTime = getTime();
sa.runByNumLoops(static_cast<ll>(1e6));
cerr << "Duration= " << getTime() - startTime << endl;
#else
ll startTime = getTime();
ll endTime = startTime + saParams->duration;
sa.runByClock(endTime - getTime());
#endif

#ifdef SA_STATS
sa.showStats();
#endif // SA_STATS

saParams->retBestCommands = sa.getBestState().mCommands;
saParams->retBestScore = sa.getBestScore();
*/

class RandXor
{
public:
	RandXor()
	{
		init();
	}

	void init()
	{
		x = 123456789;
		y = 362436069;
		z = 521288629;
		w = 88675123;
	}

	void init(ll seed)
	{
		x = 123456789;
		y = static_cast<unsigned int>(362436069 + seed * 674890201);
		z = static_cast<unsigned int>(521288629 - seed * seed * 164907679);
		w = static_cast<unsigned int>(88675123 * seed);
	}


	inline unsigned int rand()
	{
		unsigned int t;
		t = (x ^ (x << 11)); x = y; y = z; z = w; return(w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)));
	}

	unsigned int next() {
		return rand();
	}

	int nextInt(int num) {
		return rand() % num;
	}

	// [a,b)
	int next(int a, int b) {
		return a + (rand() % (b - a));
	}

	// [-d,d]
	int nextDelta(int d) {
		return next(-d, d + 1);
	}

	double nextDouble() {
		return (rand() + 0.5) * (1.0 / 4294967296.0);
	}

	template <class T>
	void randomShuffle(vector <T> &a)
	{
		const int n = SZ(a);
		for (int i = n - 1; i > 0; --i) {
			swap(a[i], a[nextInt(i + 1)]);
		}
	}

private:
	unsigned int x;
	unsigned int y;
	unsigned int z;
	unsigned int w;
};





template <class STATE>
class SimulatedAnnealing
{
public:
	const STATE& getBestState() const { return mBestState; }
	const STATE& getCurrentState() const { return mCurrentState; }
	const double getBestScore() const { return mBestScore; }
	const ll getLoopCount() const { return mLoopCount; }

	void init(
		const STATE&	state,			// STATEに、mProbのthisポインタを渡すなどの自由度を持たす狙い。こうしないとTemplateにもう１つクラスを足すしかない。
		AlgorithmMode	algorithmMode,
		ll				numSubLoops,
		ll				numTimerCheckLoops,
		double			maxStartTemp = 1.0,
		double			maxEndTemp = 1.0,
		double			minTemp = 0.001,
		double			tempExponential = 0.5,
		int				numMultiStart = 1,
		int				randomSeed = 1)
	{
		mCurrentState = state;
		mCurrentState.saInit();
		mCurrentScore = mCurrentState.saScore();

		mStartState = mBestState = mCurrentState;
		mStartScore = mBestScore = mCurrentScore;

#ifdef SA_STATS
		fprintf(stderr, "Start Score   =%8.2f time=%9lld\n", mBestScore, mLoopCount);
#endif // SA_STATS


		mAlgorithmMode = algorithmMode;

		mNumSubLoops = numSubLoops;
		mNumTimerCheckLoops = numTimerCheckLoops;

		mBestTime = 0;
		mLoopCount = 0;

		// For SA
		mRandXor.init(randomSeed);
		mMaxStartTemp = maxStartTemp;
		mMaxEndTemp = maxEndTemp;
		mMinTemp = minTemp;
		mTempExponential = tempExponential;
		mNumMultiStart = numMultiStart;

#ifdef SA_STATS
		mPhases.clear();
#endif // SA_STATS
	}

	void runByClock(const ll durationMillSecond)
	{
		updateParameters(0.0);
		mSubLoopDownCount = mNumSubLoops;
		mTimeCheckDownCount = mNumTimerCheckLoops;

		ll startTime = getTime();
		ll ellapsedTime;
		while ((ellapsedTime = getTime() - startTime) < durationMillSecond)
		{
			runOnce();
			if (mSubLoopDownCount <= 0)
			{
				updateParameters(static_cast<double>(ellapsedTime) / durationMillSecond);
				mSubLoopDownCount = mNumSubLoops;
			}

			if (mTimeCheckDownCount <= 0)
			{
				if ((ellapsedTime = getTime() - startTime) < durationMillSecond)
				{
					mTimeCheckDownCount = mNumTimerCheckLoops;
				}
				else
				{
					break;
				}
			}
		}
	}

	void runByNumLoops(const ll numLoops)
	{
		updateParameters(0.0);
		mSubLoopDownCount = mNumSubLoops;

		while (mLoopCount < numLoops)
		{
			runOnce();
			if (mSubLoopDownCount <= 0)
			{
				updateParameters(static_cast<double>(mLoopCount) / numLoops);
				mSubLoopDownCount = mNumSubLoops;
			}
		}
	}

	double computeTemp(double progressRatio) const
	{
		const double x = progressRatio;
		const double xp = pow(mNumMultiStart * x - floor(mNumMultiStart * x), mTempExponential);	// [0,1]
		const double stepped_max = (mMaxEndTemp - mMaxStartTemp) * floor(mNumMultiStart * x) / mNumMultiStart + mMaxStartTemp; // 最大値の階段関数
		const double y = mMinTemp * xp + stepped_max * (1.0 - xp);
		return y;
	}

	// 時刻・進行割合・温度の更新
	void updateParameters(double progressRatio)
	{
		mProgressRatio = progressRatio;
		// mTemp = mStartTemp + (mEndTemp - mStartTemp) * sqrt(mProgressRatio);	// TODO : 等差等比mixなど、ありうる。
		mTemp = computeTemp(mProgressRatio);
		mInvTemp = 1.0 / mTemp;
		mCurrentState.saSetProgressRatio(mProgressRatio, mTemp);
		mCurrentScore = mCurrentState.saScore();

#ifdef SA_STATS
		PHASE	phase;
		phase.neighbors.resize(STATE::NUM_NEIGHBORS); // これだと複数のSAをやるとき、まずいので、新しいneighborがでてきたとき初期化resizeという手もある。
		mPhases.push_back(phase);
#endif // SA_STATS

	}

	void runOnce()
	{
		const double deltaScore = mCurrentState.saMove();
		bool isAccepted = deltaScore > 0;	// スコアが良くなった

#ifdef SA_STATS
		const int	neighborID = mCurrentState.saGetNeighborID();
		PHASE&		phase = mPhases.back();
		NEIGHBOR&	neighbor = phase.neighbors[neighborID];
		phase.countTurns++;
#endif // SA_STATS

		// 焼きなましのときはスコアが悪くなっても
		// 最初progressRatio=0.0のときは、スコアが良くなろうが悪くなろうが、常にnextを使用
		// 最後progressRatio=1.0のときは、スコアが改善したときのみ、nextを使用
		if (mAlgorithmMode == SA && !isAccepted)
		{
			const double expCut = -6;
			double left;
			double v = deltaScore * mInvTemp;
			if (v < expCut)
			{
				isAccepted = false;
				left = 0.0;
			}
			else
			{
				left = exp(v);
				const double r = mRandXor.nextDouble();
				isAccepted = left > r;	// スコアが悪くなっても強制遷移
			}


#ifdef SA_STATS
			if (isAccepted)
			{
				neighbor.countAcceptedMinus++;
			}
			else
			{
				if (left > 0.0001)
				{
					neighbor.countRejectedMinus++;
				}
				else
				{
					neighbor.countImpossible++;
				}
			}
#endif // SA_STATS
			//			fprintf(stderr, "temp=%f deltaScore=%f left=%f ac=%d\n", mTemp, deltaScore, left, (int)isAccepted);
		}
		else
		{
#ifdef SA_STATS
			neighbor.countAcceptedPlus++;
#endif // SA_STATS
		}

		if (isAccepted)
		{
			mCurrentScore += deltaScore;
			mCurrentState.saForward();
		}
		else
		{
			mCurrentState.saResume();
		}

		// ベストスコアは
		if (mCurrentScore > mBestScore)
		{
			mBestScore = mCurrentScore;
			mBestState = mCurrentState;
			mBestTime = mLoopCount;
			fprintf(stderr, "New best Score=%8.2f time=%9lld (%.3f)\n", mBestScore, mLoopCount, mProgressRatio);
		}

		mLoopCount++;
		mSubLoopDownCount--;
		mNumTimerCheckLoops--;

#ifdef SA_STATS
		phase.totalCurrentScore += mCurrentScore;
		phase.totalBestScore += mBestScore;
#endif // SA_STATS
	}

#ifdef SA_STATS
	struct NEIGHBOR
	{
		ll getCountAll() const
		{
			return countAcceptedPlus + countAcceptedMinus + countRejectedMinus + countImpossible;
		}

		double getRatioAcceptedPlus() const
		{
			return static_cast<double>(countAcceptedPlus) / getCountAll();
		}

		double getRatioAcceptedMinus() const
		{
			return static_cast<double>(countAcceptedMinus) / getCountAll();
		}

		double getRatioRejectedMinus() const
		{
			return static_cast<double>(countRejectedMinus) / getCountAll();
		}

		double getRatioImpossible() const
		{
			return static_cast<double>(countImpossible) / getCountAll();
		}

		ll countAcceptedPlus = 0;	// deltaScoreプラスでaccepted
		ll countAcceptedMinus = 0;	// deltaScoreマイナスでaccepted
		ll countRejectedMinus = 0;	// deltaScroreマイナスでrejected
		ll countImpossible = 0;	// deltaScoreが低すぎてrejected
	};

	struct PHASE
	{
		ll countTurns = 0;
		double totalCurrentScore = 0;
		double totalBestScore = 0;
		vector <NEIGHBOR> neighbors;
		double getAverageCurrentScore() const
		{
			return totalCurrentScore / countTurns;
		}
		double getAverageBestScore() const
		{
			return totalBestScore / countTurns;
		}
	};

	vector <PHASE>	mPhases;

	void	showStats() const
	{
		for (int phaseID = 0; phaseID < SZ(mPhases); phaseID += 5)
		{
			const PHASE& phase = mPhases[phaseID];
			fprintf(stderr, "%5d %lld current=%8.2f best=%8.2f"
				, phaseID
				, phase.countTurns
				, phase.getAverageCurrentScore()
				, phase.getAverageBestScore()
			);

			for (int neighborID = 0; neighborID < SZ(phase.neighbors); ++neighborID)
			{
				const NEIGHBOR& neighbor = phase.neighbors[neighborID];

				fprintf(stderr, "[%1d]=%5lld(%4lld/%4lld/%4lld/%4lld)"
					, neighborID
					, neighbor.getCountAll()
					, neighbor.countAcceptedPlus
					, neighbor.countAcceptedMinus
					, neighbor.countRejectedMinus
					, neighbor.countImpossible
				);

			}

			fprintf(stderr, "\n");
		}
		fprintf(stderr, "LoopCount = %lld\n", mLoopCount);
	}
#endif // SA_STATS

	double	mCurrentScore = 0.0;
	double	mBestScore = 0.0;
	double	mStartScore = 0.0;

private:
	AlgorithmMode	mAlgorithmMode = SA;

	RandXor	mRandXor;

	double	mMaxStartTemp = 2.0;
	double	mMaxEndTemp = 0.001;
	double	mMinTemp = 0.001;
	double  mTempExponential = 1.0;
	int		mNumMultiStart = 1;

	double	mProgressRatio = 0.0;
	double	mTemp = 0.0;
	double  mInvTemp = 0.0;

	STATE	mCurrentState;
	STATE	mBestState;
	STATE	mStartState;

	ll		mLoopCount = 0;
	ll		mSubLoopDownCount = 0;
	ll		mNumSubLoops = 0;
	ll		mNumTimerCheckLoops = 0;
	ll		mTimeCheckDownCount = 0;
	ll		mBestTime = 0;
};


class SAState
{
public:
	SAState()
	{
	}

	SAState(const vector < vector < int > >& input_a)
	{
		a = input_a;
		H = SZ(a);
		W = SZ(a[0]);
	}

	void    saInit()
	{
		freq.resize(W);
		for (int x = 0; x < W; ++x)
		{
			freq[x].resize(H);
			for (int y = 0; y < H; ++y)
			{
				int c = a[y][x] / W;
				freq[x][c]++;
			}
		}
	}

	double  saScore()
	{
		return saRealScore();
	}

	double  saRealScore()
	{
		int ret = 0;

		for (int x = 0; x < W; ++x)
		{
			for (int c = 0; c < H; ++c)
			{
				ret -= SQ(freq[x][c]-1);
			}
		}
		return ret;
	}


	double  saMove()
	{
		mY = mRand.nextInt(H);
		mX1 = mRand.nextInt(W);
		mX2 = mRand.nextInt(W - 1);
		if (mX2 >= mX1) mX1++;

		const int c1 = a[mY][mX1] / W;
		const int c2 = a[mY][mX2] / W;

		double deltaScore = 0.0;
		deltaScore -= SQ(freq[mX1][c1] - 1);
		deltaScore -= SQ(freq[mX1][c2] - 1);
		deltaScore += SQ(freq[mX1][c1]-1 - 1);
		deltaScore += SQ(freq[mX1][c2]+1 - 1);

		deltaScore -= SQ(freq[mX2][c1] - 1);
		deltaScore -= SQ(freq[mX2][c2] - 1);
		deltaScore += SQ(freq[mX2][c1] + 1 - 1);
		deltaScore += SQ(freq[mX2][c2] - 1 - 1);

		return -deltaScore;
	}

	void    saForward()
	{
		const int c1 = a[mY][mX1] / W;
		const int c2 = a[mY][mX2] / W;
		freq[mX1][c1]--;
		freq[mX2][c2]--;
		freq[mX1][c2]++;
		freq[mX2][c1]++;

		swap(a[mY][mX1], a[mY][mX2]);
	}

	void    saResume()
	{
	}

	void    saSetProgressRatio(double progressRatio, double temp)
	{
	}

public:
	vector <vector <int> > a;
	vector <vector <int> > freq;
	int H;
	int W;
	int mY;
	int mX1;
	int mX2;

	RandXor mRand;

};

void tenchi(vector <vector <int>>& vs)
{
	const int H = SZ(vs);
	const int W = SZ(vs[0]);

	vector <vector <int>> ret(W, vector <int>(H));
	for (int y = 0; y < W; y++)
	{
		for (int x = 0; x < H; x++)
		{
			ret[y][x] = vs[x][y];
		}
	}

	vs = ret;
}

void printMatrix(const vector <vector <int> > & a)
{
	for (int y = 0; y < SZ(a); ++y)
	{
		for (int x = 0; x < SZ(a[0]); ++x)
		{
			if (x != 0) printf(" ");
			printf("%d", a[y][x]);
		}
		printf("\n");
	}
}

int main() 
{
	int H, W;
	cin >> H >> W;

#if 1 // 本番
	vector <vector <int> > a;
	a.resize(H);

	for (int y = 0; y < H; ++y)
	{
		for (int x = 0; x < W; ++x)
		{
			int tmp;
			cin >> tmp;
			tmp--;
			a[y].push_back(tmp);
		}
	}

#else // テスト用
	vector <vector <int> > a;
	a.resize(H);
	for (int y = 0; y < H; ++y)
	{
		a[y].resize(W);
	}

	vector <int> input;
	for (int i = 0; i < H*W; ++i)
	{
		input.push_back(i);
	}

	random_shuffle(input.begin(), input.end());

	for (int y = 0; y < H; ++y)
	{
		for (int x = 0; x < W; ++x)
		{
			a[y][x] = input[y*W + x];
		}
	}
#endif


	double          maxTemp = 0.5;
	double          minTemp = 0.001;
	double          tempExponential = 1.0;
	int             numMultiStart = 1;

	SAState state(a);
	SimulatedAnnealing<SAState> sa;

	sa.init(state, SA, 5000, 1000, maxTemp, maxTemp, minTemp, tempExponential, numMultiStart);
//	ll startTime = getTime();
	sa.runByNumLoops(static_cast<ll>(1e8));

	vector <vector <int> > ret = sa.getBestState().a;
	for (int y = 0; y < H; ++y)
	{
		for (int x = 0; x < W; ++x)
		{
			ret[y][x]++;
		}
	}

	printMatrix(ret);

	tenchi(ret);
	for (int y = 0; y < SZ(ret); ++y)
	{
		sort(ret[y].begin(), ret[y].end());
	}
	tenchi(ret);

	printMatrix(ret);
	for (int y = 0; y < H; ++y)
	{
		sort(ret[y].begin(), ret[y].end());
	}

//	printMatrix(ret);


//	cerr << "Duration= " << getTime() - startTime << endl;


	return 0;
}


Compiler output:
a.cc: In member function 'void SimulatedAnnealing<STATE>::runByClock(ll)':
a.cc:239:32: error: there are no arguments to 'getTime' that depend on a template parameter, so a declaration of 'getTime' must be available [-fpermissive]
  239 |                 ll startTime = getTime();
      |                                ^~~~~~~
a.cc:239:32: note: (if you use '-fpermissive', G++ will accept your code, but allowing the use of an undeclared name is deprecated)
a.cc:241:40: error: there are no arguments to 'getTime' that depend on a template parameter, so a declaration of 'getTime' must be available [-fpermissive]
  241 |                 while ((ellapsedTime = getTime() - startTime) < durationMillSecond)
      |                                        ^~~~~~~
a.cc:252:53: error: there are no arguments to 'getTime' that depend on a template parameter, so a declaration of 'getTime' must be available [-fpermissive]
  252 |                                 if ((ellapsedTime = getTime() - startTime) < durationMillSecond)
      |                                                     ^~~~~~~