slatenails

find_rings_3.cc

Oct 16th, 2013
155
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <cstdio>
  2. #include <cstdlib>
  3. #include <cmath>
  4. #include <cstring>
  5. #include <vector>
  6.  
  7. int g_stack = 0;
  8.  
  9. // One ring component within a solution
  10. struct Component
  11. {   int num;
  12.     double scale;
  13. };
  14.  
  15. // =============================================================================
  16. // Solution class
  17. // -----------------------------------------------------------------------------
  18. class Solution
  19. {   public:
  20.         // Components of this solution
  21.         inline const std::vector<Component>& components() const
  22.         {   return m_components;
  23.         }
  24.  
  25.         // Add a component to this solution
  26.         void addComponent (const Component& a)
  27.         {   m_components.push_back (a);
  28.         }
  29.  
  30.         // Compare solutions
  31.         bool operator> (const Solution& other) const
  32.         {   // If this solution has less components than the other one, this one
  33.             // is definitely better.
  34.             if (components().size() < other.components().size())
  35.                 return true;
  36.  
  37.             // vice versa
  38.             if (other.components().size() < components().size())
  39.                 return false;
  40.  
  41.             // Calculate the maximum ring number. Since the solutions have equal
  42.             // ring counts, the solutions with lesser maximum rings should result
  43.             // in cleaner code and less new primitives, right?
  44.             int maxA = 0,
  45.                 maxB = 0;
  46.  
  47.             for (size_t i = 0; i < components().size(); ++i)
  48.             {   if (components()[i].num > maxA)
  49.                     maxA = components()[i].num;
  50.  
  51.                 if (other.components()[i].num > maxB)
  52.                     maxB = other.components()[i].num;
  53.             }
  54.  
  55.             if (maxA < maxB)
  56.                 return true;
  57.  
  58.             if (maxB < maxA)
  59.                 return false;
  60.  
  61.             // Solutions have equal rings and equal maximum ring numbers. Let's
  62.             // just say this one is better, at this point it does not matter which
  63.             // one is chosen.
  64.             return true;
  65.         }
  66.  
  67.     private:
  68.         std::vector<Component> m_components;
  69. };
  70.  
  71. std::vector<Solution> solutions;
  72. Solution bestSolution;
  73.  
  74. template<class T> inline T abs (T a)
  75. {   return (a > 0) ? a : -a;
  76. }
  77.  
  78. inline bool isZero (double a)
  79. {   return abs<double> (a) < 0.0001;
  80. }
  81.  
  82. inline bool isInteger (double a)
  83. {   return isZero (a - floor (a));
  84. }
  85.  
  86. // =============================================================================
  87. // ringFinder
  88. //
  89. // This is the main algorithm of the ring finder. It tries to use math to find
  90. // the one ring between r0 and r1. If it fails (the ring number is non-integral),
  91. // it finds an intermediate radius (ceil of the ring number times scale) and
  92. // splits the radius at this point, calling this function again to try find the
  93. // rings between r0 - r and r - r1.
  94. //
  95. // This does not always yield into usable results. If at some point r == r0 or
  96. // r == r1, there is no hope of finding the rings, at least with this algorithm,
  97. // as it would fall into an infinite recursion.
  98. // -----------------------------------------------------------------------------
  99. bool ringFinderRecursor (double r0, double r1, Solution& currentSolution)
  100. {   char tabs[64];
  101.     memset (tabs, '\t', g_stack);
  102.     tabs[g_stack] = '\0';
  103.  
  104.     // Don't recurse too deep.
  105.     if (g_stack >= 5)
  106.         return false;
  107.  
  108.     // Find the scale and number of a ring between r1 and r0.
  109.     double scale = r1 - r0;
  110.     double num = r0 / scale;
  111.  
  112.     // If the ring number is integral, we have found a fitting ring to r0 -> r1!
  113.     if (isInteger (num))
  114.     {   Component cmp;
  115.         cmp.scale = scale;
  116.         cmp.num = (int) round (num);
  117.         currentSolution.addComponent (cmp);
  118.  
  119.         // If we're still at the first recursion, this is the only
  120.         // ring and there's nothing left to do. Guess we found the winner.
  121.         if (g_stack == 0)
  122.         {   solutions.push_back (currentSolution);
  123.             return true;
  124.         }
  125.     }
  126.     else
  127.     {   // Try find solutions by splitting the ring in various positions.
  128.         if (isZero (r1 - r0))
  129.             return false;
  130.  
  131.         double interval;
  132.  
  133.         // Determine interval. The smaller delta between radii, the more precise
  134.         // interval should be used. We can't really use a 0.5 increment when
  135.         // calculating rings to 10 -> 105... that would take ages to process!
  136.         if (r1 - r0 < 0.5)
  137.             interval = 0.1;
  138.         else if (r1 - r0 < 10)
  139.             interval = 0.5;
  140.         else if (r1 - r0 < 50)
  141.             interval = 1;
  142.         else
  143.             interval = 5;
  144.  
  145.         // Now go through possible splits and try find rings for both segments.
  146.         std::vector<double> radii;
  147.  
  148.         // r1 / 2 is a particularly good split point, as r1 / 2 - r1 can be
  149.         // filled with just one ring, possibly covering a lot of space.s
  150.         if (r1 / 2 > r0)
  151.             radii.push_back (r1 / 2);
  152.  
  153.         for (double r = r0 + interval; r < r1; r += interval)
  154.             radii.push_back (r);
  155.  
  156.         for (std::vector<double>::iterator it = radii.begin(); it != radii.end(); ++it)
  157.         {   double r = *it;
  158.             Solution sol = currentSolution;
  159.  
  160.             g_stack++;
  161.             bool res = ringFinderRecursor (r0, r, sol) && ringFinderRecursor (r, r1, sol);
  162.             g_stack--;
  163.  
  164.             if (res)
  165.             {   // We succeeded in finding radii for this segment. If the stack is 0, this
  166.                 // is the first recursion to this function. Thus there are no more ring segments
  167.                 // to process and we can add the solution.
  168.                 //
  169.                 // If not, when this function ends, it will be called again with more arguments.
  170.                 // Accept the solution to this segment by setting currentSolution to sol, and
  171.                 // return true to continue processing.
  172.                 if (g_stack == 0)
  173.                     solutions.push_back (sol);
  174.                 else
  175.                 {   currentSolution = sol;
  176.                     return true;
  177.                 }
  178.             }
  179.         }
  180.  
  181.         return false;
  182.     }
  183.  
  184.     return true;
  185. }
  186.  
  187. // =============================================================================
  188. // Main function. Call this with r0 and r1. If this returns true, use bestSolution
  189. // for the solution that was presented.
  190. // -----------------------------------------------------------------------------
  191. bool findRings (double r0, double r1)
  192. {   Solution sol;
  193.  
  194.     // Recurse in and try find solutions.
  195.     ringFinderRecursor (r0, r1, sol);
  196.  
  197.     // Compare the solutions and find the best one. The solution class has an operator>
  198.     // overload to compare two solutions.
  199.     const Solution* best = NULL;
  200.     for (std::vector<Solution>::iterator solp = solutions.begin(); solp != solutions.end(); ++solp)
  201.     {   const Solution& sol = *solp;
  202.  
  203.         if (best == NULL || sol > *best)
  204.             best = &sol;
  205.     }
  206.  
  207.     // If we found one best solution, return it now.
  208.     if (best)
  209.     {   bestSolution = *best;
  210.         return true;
  211.     }
  212.  
  213.     // This should only happen if no solutions were found in the first place.
  214.     return false;
  215. }
  216.  
  217. int main (int argc, char* argv[])
  218. {   char* endptr0, *endptr1;
  219.     double r0, r1;
  220.  
  221.     // Validate input
  222.     if (argc != 3)
  223.     {   fprintf (stderr, "usage: %s <r0> <r1>\n", argv[0]);
  224.         return 1;
  225.     }
  226.  
  227.     r0 = strtod (argv[1], &endptr0);
  228.     r1 = strtod (argv[2], &endptr1);
  229.  
  230.     if (*endptr0 != '\0' || *endptr1 != '\0')
  231.     {   fprintf (stderr, "error: arguments are not numbers\n");
  232.         return 1;
  233.     }
  234.  
  235.     if (r0 > r1)
  236.         std::swap (r0, r1);
  237.  
  238.     // Use the algorithm to find a solution
  239.     if (!findRings (r0, r1))
  240.     {   printf ("can't find rings between %g and %g\n", r0, r1);
  241.         return 1;
  242.     }
  243.  
  244.     // Iterate through the solution and print the rings involved.
  245.     for (std::vector<Component>::const_iterator it = bestSolution.components().begin(); it != bestSolution.components().end(); ++it)
  246.         printf ("Ring %d, scale by %g (%g -> %g)\n", it->num, it->scale, it->num * it->scale, (it->num + 1) * it->scale);
  247. }
RAW Paste Data