Timc 2021 新年題 Q3

readme.md

這是整理在這篇文章裡的程式碼： https://blog.cruciferslab.net/?p=14

factsqrt.cpp

#include <iostream>
#include <cmath>
#include <map>
#include <queue>
#include <tuple>
using namespace std;

// switch the comment below to see log printed to stderr
//#define LOGERR(x) x
#define LOGERR(x) 

double MAX_LOG = 35; // e^35 < 2^52
const double PI = 3.141592653589793238;

struct _source
{
	int64_t factnum;
	int nsqrt;
};

struct _qitem
{
	int steps;
	int64_t number;
	bool operator < (const _qitem& other) const
	{
		return steps > other.steps || (steps == other.steps && number > other.number);
	}
	operator pair<int, int64_t> () const
	{
		return {steps, number};
	}
};

map<int64_t, int> steps;
map<int64_t, _source> source;
priority_queue<_qitem> q;

double logfactorial(int n)
{
	if(n <= 17) // 17! < e^35 < 18!
	{
		double v = 1;
		for(int i = 2; i <= n; i++)
		{
			v *= i;
		}
		return log(v);
	}
	else
	{
		// Use Ramanujan's approximation
		// https://math.stackexchange.com/a/138326
		return (n * log(n * 1.) - n + log(n * (1. + 4.*n*(1. + 2.*n)) + 1./30)/6 + log(PI)/2);
	}
}

int main()
{
	q.push({0, 2021});
	source[2021] = {-1, 0};
	steps[2021] = 0;
	double v = sqrt(2021);
	for(int i = 1; v >= 3; ++i, v = sqrt(v))
	{
		int64_t iv = (int64_t)(floor(v));
		q.push({i+1, iv});
		source[iv] = {-1, i};
		steps[iv] = i+1;
		LOGERR(cerr << "Update " << iv << " to " << (i+1) << " steps" << endl);
	}
	
	int s;
	int64_t	num;
	while(tie(s, num) = (pair<int,int64_t>)q.top(), q.pop(), num != 110)
	{
		if(steps.find(num) != steps.end() && steps[num] < s) continue;
		LOGERR(cerr << "Doing " << num << "! with " << s << " steps" << endl);
		steps[num] = s;
		double logfact = logfactorial(num);
		int sqcounter = 0;
		while(logfact > MAX_LOG)
		{
			++sqcounter;
			logfact /= 2;
			LOGERR(cerr << "  ln(" << num << "!) / 2^" << sqcounter << " ~ " << logfact << endl);
		}
		double value = exp(logfact);
		while(value >= 3)
		{
			LOGERR(cerr << "  " << num << "! with sqrt^" << sqcounter << " ~ " << value << endl);
			int64_t iv = (int64_t)(floor(value));
			int newstep = s+sqcounter+1;
			if(sqcounter > 0) newstep++;
			if(steps.find(iv) == steps.end() || steps[iv] > newstep)
			{
				q.push({newstep, iv});
				source[iv] = {num, sqcounter};
				steps[iv] = newstep;
				LOGERR(cerr << "  Update " << iv << " to " << newstep << " steps" << endl);
			}
			++sqcounter;
			value = sqrt(value);
		}
	}
	
	cout << "Found 110 at " << s << " steps:" << endl;
	int64_t x = 110;
	while(x != -1)
	{
		auto src = source[x];
		cout << x << " <- ";
		if(src.factnum == -1)
			cout << "2021";
		else
			cout << src.factnum << "!";
		cout << ", sqrt * " << src.nsqrt << endl;
		x = src.factnum;
	}
	
	return 0;
}

factsqrt2.cpp

#include <iostream>
#include <cmath>
#include <map>
#include <queue>
#include <tuple>
using namespace std;

// switch the comment below to see log printed to stderr
#define LOGERR(x) x
//#define LOGERR(x) 

const double MAX_LOG = 35; // e^35 < 2^52
const int MAX_FACT = 17; // 17! < e^35 < 18!
const double PI = 3.141592653589793238;

struct _source
{
	int64_t factnum;
	int nsqrt;
};

struct _qitem
{
	int steps;
	int64_t number;
	bool operator < (const _qitem& other) const
	{
		return steps > other.steps || (steps == other.steps && number > other.number);
	}
	operator pair<int, int64_t> () const
	{
		return {steps, number};
	}
};

map<int64_t, int> steps;
map<int64_t, _source> source;
priority_queue<_qitem> q;

// Assume called with n > MAX_FACT
double logfactorial(int64_t n)
{
	// Use Ramanujan's approximation
	// https://math.stackexchange.com/a/138326
	return (n * log(n * 1.) - n + log(n * (1. + 4.*n*(1. + 2.*n)) + 1./30)/6 + log(PI)/2);
}

int main()
{
	q.push({0, 2021});
	source[2021] = {-1, 0};
	steps[2021] = 0;
	double v = sqrt(2021);
	for(int i = 1; v >= 3; ++i, v = sqrt(v))
	{
		int64_t iv = (int64_t)(floor(v));
		q.push({i+1, iv});
		source[iv] = {-1, i};
		steps[iv] = i+1;
		LOGERR(cerr << "Update " << iv << " to " << (i+1) << " steps" << endl);
	}
	
	int s;
	int64_t	num;
	while(tie(s, num) = (pair<int,int64_t>)q.top(), q.pop(), num != 110)
	{
		if(steps.find(num) != steps.end() && steps[num] < s) continue;
		LOGERR(cerr << "Doing " << num << "! with " << s << " steps" << endl);
		steps[num] = s;
		double value;
		int sqcounter;
		if(num > MAX_FACT)
		{
			double logfact = logfactorial(num);
			if(isnan(logfact))
			{
				LOGERR(cerr << "  WARNING: logfactorial returns NAN" << endl);
			}
			sqcounter = 0;
			while(logfact > MAX_LOG)
			{
				++sqcounter;
				logfact /= 2;
				LOGERR(cerr << "  ln(" << num << "!) / 2^" << sqcounter << " ~ " << logfact << endl);
			}
			value = exp(logfact);
		}
		else
		{
			value = 1;
			for(int i = 2; i <= num; i++) value *= i;
			sqcounter = 0;
		}
		while(value >= 3)
		{
			LOGERR(cerr << "  " << num << "! with sqrt^" << sqcounter << " ~ " << value << endl);
			int64_t iv = (int64_t)(floor(value));
			int newstep = s+sqcounter+1;
			if(sqcounter > 0) newstep++;
			if(steps.find(iv) == steps.end() || steps[iv] > newstep)
			{
				q.push({newstep, iv});
				source[iv] = {num, sqcounter};
				steps[iv] = newstep;
				LOGERR(cerr << "  Update " << iv << " to " << newstep << " steps" << endl);
			}
			++sqcounter;
			value = sqrt(value);
		}
	}
	
	cout << "Found 110 at " << s << " steps:" << endl;
	int64_t x = 110;
	while(x != -1)
	{
		auto src = source[x];
		cout << x << " <- ";
		if(src.factnum == -1)
			cout << "2021";
		else
			cout << src.factnum << "!";
		cout << ", sqrt * " << src.nsqrt << endl;
		x = src.factnum;
	}
	
	return 0;
}