Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Toy problemsは役立たずか(Floydの問題・ネタバレ) http://blog.unfindable.net/archives/6103
//Are Toy Problems Useful?
//in Selected Papers on Computer Science by Knuth (p.172)
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <set>
#include <unordered_map> // faster than map
#include <algorithm>
#include <ctime>
//#define SQRT sqrtl
//typedef long double real;
#define SQRT sqrt
typedef double real;
typedef unsigned int uint;
using namespace std;
const real DELTA = 1000.;
const uint N = 50;
const uint sizeA = 20;
int getKey(real x)
{
return static_cast<int>((x - floor(x)) * 1e8);
}
real getSum(uint code, vector<real>* group, bool ignoreLast)
{
real sum = 0;
uint size = group->size();
if (ignoreLast) --size;
for (uint i = 0; i < size; ++i) if (code >> i & 1) sum += (*group)[i];
return sum;
}
void printSolution(set<int>* solution)
{
real sum = 0;
stringstream ss;
for (auto i : *solution) {
sum += SQRT(i);
ss << i << ",";
}
string str = ss.str();
cout << setprecision(30) << sum << ": {" << str.substr(0, str.size()-1) << "}" << endl;
}
int main()
{
auto start = clock();
real dummy;
real goal = 0;
vector<int> all, a, b, c;
vector<real> groupA, groupB, groupC; //平方根を計算しておく
for (uint i = 1; i <= N; ++i) {
all.push_back(i); //{1, 2, ..., 50}
real t = SQRT(i);
goal += t;
if (t == floor(t)) {
c.push_back(i); //{1, 4, 9, 16, 25, 36, 49}
groupC.push_back(t);
}
else if (a.size() < sizeA) {
a.push_back(i); //{2, 3, 5, 6, 7, 8, 10, 11, 12, 13,...
groupA.push_back(t);
}
else {
b.push_back(i);
groupB.push_back(t);
}
}
goal /= 2.;
real fracPartOfGoal = modf(goal, &dummy);
//aのsubset sumを計算し、記憶する
unordered_multimap<int, pair<int, real>* > subsetSumOfA; //キー:概数、値:codeAと和のペア
for (uint codeA = 0; codeA < (1u << groupA.size()); ++codeA) {
real sumA = getSum(codeA, &groupA, false);
int key = getKey(sumA);
auto p = new pair<int, real>(codeA, sumA);
subsetSumOfA.insert(pair<int, pair<int, real>* >(key, p));
}
//bのsubset sumを計算する
real minDelta = DELTA;
int codeAforMinDelta = 0, codeBforMinDelta = 0;
real sumAforMinDelta = 0, sumBforMinDelta = 0;
//#pragma omp parallel
for (uint codeB = 0; codeB < (1u << (groupB.size() - 1)); ++codeB) { //bの1つの要素は決めておいてよい
real sumB = getSum(codeB, &groupB, true);
real t = fabs(fracPartOfGoal - modf(sumB, &dummy)); //興味があるのは非常に近い場合だけだからこれでいい
auto range = subsetSumOfA.equal_range(getKey(t));
for (auto i = range.first; i != range.second; ++i) {
int codeA = i->second->first;
real sumA = i->second->second;
real delta = fabs(fracPartOfGoal - modf(sumA, &dummy) - modf(sumB, &dummy));
if (delta < minDelta) {
//#pragma omp critical
{
minDelta = delta;
codeAforMinDelta = codeA;
codeBforMinDelta = codeB;
sumAforMinDelta = sumA;
sumBforMinDelta = sumB;
}
}
}
}
//整数部分を揃える
real sum = sumAforMinDelta + sumBforMinDelta;
minDelta = DELTA;
int codeCforMinDelta = 0;
for (uint codeC = 0; codeC < (1u << groupC.size()); ++codeC) {
real delta = fabs(goal - sum - getSum(codeC, &groupC, false));
if (delta < minDelta) {
minDelta = delta;
codeCforMinDelta = codeC;
}
}
std::cout << setprecision(30) << goal << ": goal" << endl;
auto solution1 = set<int>(); //ソート済みになる
auto solution2 = set<int>();
for (uint i = 0; i < groupA.size(); ++i) if (codeAforMinDelta >> i & 1) solution1.insert(a[i]);
for (uint i = 0; i < groupB.size() - 1; ++i) if (codeBforMinDelta >> i & 1) solution1.insert(b[i]);
for (uint i = 0; i < groupC.size(); ++i) if (codeCforMinDelta >> i & 1) solution1.insert(c[i]);
set_difference(all.begin(), all.end(), solution1.begin(), solution1.end(), inserter(solution2, solution2.end()));
printSolution(&solution1);
printSolution(&solution2);
clock_t finish = clock();
real duration = real(finish - start) / CLOCKS_PER_SEC;
cout << "Elapsed time: " << setprecision(3) << duration << " sec." << endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment