fork download
copy 
  1. //Are Toy Problems Useful?
  2. //in Selected Papers on Computer Science by Knuth (p.172)
  3.  
  4. #include <iostream>
  5. #include <sstream>
  6. #include <iomanip>
  7. #include <vector>
  8. #include <set>
  9. #include <unordered_map> // faster than map
  10. #include <algorithm>
  11. #include <ctime>
  12. //#define SQRT sqrtl
  13. //typedef long double real;
  14. #define SQRT sqrt
  15. typedef double real;
  16. typedef unsigned int uint;
  17. using namespace std;
  18.  
  19. const real DELTA = 1000.;
  20. const uint N = 50;
  21. const uint sizeA = 20;
  22.  
  23. int getKey(real x)
  24. {
  25.   return static_cast<int>((x - floor(x)) * 1e8);
  26. }
  27.  
  28. real getSum(uint code, vector<real>* group, bool ignoreLast)
  29. {
  30.   real sum = 0;
  31.   uint size = group->size();
  32.   if (ignoreLast) --size;
  33.   for (uint i = 0; i < size; ++i) if (code >> i & 1) sum += (*group)[i];
  34.   return sum;
  35. }
  36.  
  37. void printSolution(set<int>* solution)
  38. {
  39.   real sum = 0;
  40.   stringstream ss;
  41.   for (auto i : *solution) {
  42.     sum += SQRT(i);
  43.     ss << i << ",";
  44.   }
  45.   string str = ss.str();
  46.   cout << setprecision(30) << sum << ": {" << str.substr(0, str.size()-1) << "}" << endl;
  47. }
  48.  
  49. int main()
  50. {
  51.   auto start = clock();
  52.  
  53.   real dummy;
  54.   real goal = 0;
  55.   vector<int> all, a, b, c;
  56.   vector<real> groupA, groupB, groupC; //平方根を計算しておく
  57.   for (uint i = 1; i <= N; ++i) {
  58.     all.push_back(i); //{1, 2, ..., 50}
  59.     real t = SQRT(i);
  60.     goal += t;
  61.     if (t == floor(t)) {
  62.       c.push_back(i); //{1, 4, 9, 16, 25, 36, 49}
  63.       groupC.push_back(t);
  64.     }
  65.     else if (a.size() < sizeA) {
  66.       a.push_back(i); //{2, 3, 5, 6, 7, 8, 10, 11, 12, 13,...
  67.       groupA.push_back(t);
  68.     }
  69.     else {
  70.       b.push_back(i);
  71.       groupB.push_back(t);
  72.     }
  73.   }
  74.   goal /= 2.;
  75.   real fracPartOfGoal = modf(goal, &dummy);
  76.  
  77.   //aのsubset sumを計算し、記憶する
  78.   unordered_multimap<int, pair<int, real>* > subsetSumOfA; //キー：概数、値：codeAと和のペア
  79.   for (uint codeA = 0; codeA < (1u << groupA.size()); ++codeA) {
  80.     real sumA = getSum(codeA, &groupA, false);
  81.     int key = getKey(sumA);
  82.     auto p = new pair<int, real>(codeA, sumA);
  83.     subsetSumOfA.insert(pair<int, pair<int, real>* >(key, p));
  84.   }
  85.  
  86.   //bのsubset sumを計算する
  87.   real minDelta = DELTA;
  88.   int codeAforMinDelta = 0, codeBforMinDelta = 0;
  89.   real sumAforMinDelta = 0, sumBforMinDelta = 0;
  90. //#pragma omp parallel
  91.   for (uint codeB = 0; codeB < (1u << (groupB.size() - 1)); ++codeB) { //bの1つの要素は決めておいてよい
  92.     real sumB = getSum(codeB, &groupB, true);
  93.     real t = fabs(fracPartOfGoal - modf(sumB, &dummy)); //興味があるのは非常に近い場合だけだからこれでいい
  94.     auto range = subsetSumOfA.equal_range(getKey(t));
  95.     for (auto i = range.first; i != range.second; ++i) {
  96.       int codeA = i->second->first;
  97.       real sumA = i->second->second;
  98.       real delta = fabs(fracPartOfGoal - modf(sumA, &dummy) - modf(sumB, &dummy));
  99.       if (delta < minDelta) {
  100. //#pragma omp critical
  101.         {
  102.           minDelta = delta;
  103.           codeAforMinDelta = codeA;
  104.           codeBforMinDelta = codeB;
  105.           sumAforMinDelta = sumA;
  106.           sumBforMinDelta = sumB;
  107.         }
  108.       }
  109.     }
  110.   }
  111.  
  112.   //整数部分を揃える
  113.   real sum = sumAforMinDelta + sumBforMinDelta;
  114.   minDelta = DELTA;
  115.   int codeCforMinDelta = 0;
  116.   for (uint codeC = 0; codeC < (1u << groupC.size()); ++codeC) {
  117.     real delta = fabs(goal - sum - getSum(codeC, &groupC, false));
  118.     if (delta < minDelta) {
  119.       minDelta = delta;
  120.       codeCforMinDelta = codeC;
  121.     }
  122.   }
  123.  
  124.   std::cout << setprecision(30) << goal << ": goal" << endl;
  125.   auto solution1 = set<int>(); //ソート済みになる
  126.   auto solution2 = set<int>();
  127.   for (uint i = 0; i < groupA.size(); ++i) if (codeAforMinDelta >> i & 1) solution1.insert(a[i]);
  128.   for (uint i = 0; i < groupB.size() - 1; ++i) if (codeBforMinDelta >> i & 1) solution1.insert(b[i]);
  129.   for (uint i = 0; i < groupC.size(); ++i) if (codeCforMinDelta >> i & 1) solution1.insert(c[i]);
  130.   set_difference(all.begin(), all.end(), solution1.begin(), solution1.end(), inserter(solution2, solution2.end()));
  131.   printSolution(&solution1);
  132.   printSolution(&solution2);
  133.  
  134.   clock_t finish = clock();
  135.   real duration = real(finish - start) / CLOCKS_PER_SEC;
  136.   cout << "Elapsed time: " << setprecision(3) << duration << " sec." << endl;
  137.   return 0;
  138. }
Success #stdin #stdout 2.56s 35776KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
119.517900301760391812422312796: goal
119.517900301760320758148736786: {1,4,5,7,8,9,12,13,15,16,20,25,27,30,31,33,37,38,41,43,44,46,47,48,49}
119.517900301760462866695888806: {2,3,6,10,11,14,17,18,19,21,22,23,24,26,28,29,32,34,35,36,39,40,42,45,50}
Elapsed time: 2.47 sec.
https://ideone.com/HwbkTo
language:
C++14 (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!