fork(1) download 
  1. /* C++11 implementation for "Selection in X + Y and matrices with sorted
  2.    rows and columns" by A. Mirzaian and E. Arjomandi.
  3.    http://w...content-available-to-author-only...u.ca/~andy/pubs/X%2BY.pdf
  4.  */
  5.  
  6. #include <algorithm>
  7. #include <vector>
  8. #include <iostream>
  9. #include <assert.h>
  10.  
  11. int qceil(int x)
  12. {
  13.   return (x + 3) / 4;
  14. }
  15.  
  16. int qhigh(int x, int n)
  17. {
  18.   if (n & 1) // odd
  19.     return qceil(x + 2 * n + 1);
  20.   else // even
  21.     return n + 1 + qceil(x);
  22. }
  23.  
  24. typedef std::vector<int> Vec;
  25.  
  26. struct Ranks
  27. {
  28.   int ra_less;
  29.   int rb_more;
  30.   Vec l;
  31.  
  32.   Ranks(const Vec& a1, const Vec& a2, int a, int b)
  33.   {
  34.     int n = a1.size() - 1;
  35.     ra_less = 0;
  36.     rb_more = n * n;
  37.     l.reserve(12 * n + 1);
  38.     l.push_back(0);
  39.  
  40.     int j = n;
  41.     for (int i = 1; i <= n; ++i)
  42.     {
  43.       while (j && a1[i] + a2[j] <= a)
  44. 	--j;
  45.  
  46.       ra_less += j;
  47.       int jj = j;
  48.  
  49.       while (jj && a1[i] + a2[jj] < b)
  50.       {
  51. 	l.push_back(a1[i] + a2[jj]);
  52. 	--jj;
  53.       }
  54.  
  55.       rb_more -= jj;
  56.     }
  57.   }
  58. };
  59.  
  60. Vec half(const Vec& a)
  61. {
  62.   Vec res;
  63.   res.reserve(2 + (a.size() - 1) / 2);
  64.   res.push_back(0);
  65.  
  66.   for (int i = 1; i < a.size(); i += 2)
  67.     res.push_back(a[i]);
  68.  
  69.   if (a.size() & 1)
  70.     res.push_back(a[a.size() - 1]);
  71.  
  72.   return res;
  73. }
  74.  
  75. struct AB
  76. {
  77.   int a;
  78.   int b;
  79. };
  80.  
  81. AB biselect(const Vec& a1, const Vec& a2, int k1, int k2)
  82. {
  83.   AB res;
  84.   int n = a1.size() - 1;
  85.   assert(n > 0);
  86.  
  87.   if (n == 1)
  88.   {
  89.     res.a = res.b = a1[1] + a2[1];
  90.   }
  91.   else if (n == 2)
  92.   {
  93.     Vec l {a1[1] + a2[1], a1[1] + a2[2], a1[2] + a2[1], a1[2] + a2[2]};
  94.     std::nth_element(l.begin(),
  95. 		     l.end() - k1,
  96. 		     l.end());
  97.     res.a = l[4 - k1];
  98.     std::nth_element(l.begin(),
  99. 		     l.end() - k2,
  100. 		     l.end());
  101.     res.b = l[4 - k2];
  102.   }
  103.   else
  104.   {
  105.     int k1_half = qhigh(k1, n);
  106.     int k2_half = qceil(k2);
  107.     AB ab = biselect(half(a1), half(a2), k1_half, k2_half);
  108.     Ranks ranks(a1, a2, ab.a, ab.b);
  109.     int r1 = k1 + ranks.rb_more - n * n;
  110.     int r2 = k2 + ranks.rb_more - n * n;
  111.  
  112.     if (ranks.ra_less <= k1 - 1)
  113.       res.a = ab.a;
  114.     else if (r1 <= 0)
  115.       res.a = ab.b;
  116.     else
  117.     {
  118.       std::nth_element(ranks.l.begin() + 1,
  119. 		       ranks.l.end() - r1,
  120. 		       ranks.l.end());
  121.       res.a = *(ranks.l.end() - r1);
  122.     }
  123.  
  124.     if (ranks.ra_less <= k2 - 1)
  125.       res.b = ab.a;
  126.     else if (r2 <= 0)
  127.       res.b = ab.b;
  128.     else
  129.     {
  130.       std::nth_element(ranks.l.begin() + 1,
  131. 		       ranks.l.end() - r2,
  132. 		       ranks.l.end());
  133.       res.b = *(ranks.l.end() - r2);
  134.     }
  135.  }
  136.  
  137.   return res;
  138. }
  139.  
  140. int select(const Vec& a1, const Vec& a2, int k)
  141. {
  142.   assert(a1.size() == a2.size());
  143.   AB ab = biselect(a1, a2, k, k);
  144.   assert(ab.a == ab.b);
  145.   return ab.a;
  146. }
  147.  
  148. int check(const Vec& a1, const Vec& a2, int k)
  149. {
  150.   int n = a1.size() - 1;
  151.   Vec sum;
  152.   sum.reserve(n * n);
  153.  
  154.   for (int i = 1; i <= n; ++i)
  155.     for (int j = 1; j <= n; ++j)
  156.       sum.push_back(a1[i] + a2[j]);
  157.  
  158.   std::nth_element(sum.begin(),
  159. 		   sum.end() - k,
  160. 		   sum.end());
  161.   return *(sum.end() - k);
  162. }
  163.  
  164. #include <random>
  165.  
  166. int main()
  167. {
  168.   std::random_device rd;
  169.   std::default_random_engine e1(rd());
  170.   std::uniform_int_distribution<int> uniform_dist(1, 6);
  171.  
  172.   int n = 100000;
  173.   Vec a1, a2;
  174.   a1.reserve(n+1);
  175.   a2.reserve(n+1);
  176.   int x = 1000;
  177.  
  178.   a1.push_back(0);
  179.   for (int i = 0; i < n; ++i)
  180.   {
  181.     x -= uniform_dist(e1);
  182.     a1.push_back(x);
  183.   }
  184.  
  185.   x = 1000;
  186.  
  187.   a2.push_back(0);
  188.   for (int i = 0; i < n; ++i)
  189.   {
  190.     x -= uniform_dist(e1);
  191.     a2.push_back(x);
  192.   }
  193.  
  194.   std::cout << select(a1, a2, 1 + n * n / 7) << '\n'
  195.     /*<< check(a1, a2, 1 + n * n / 7) << '\n'*/;
  196.  
  197.   return 0;
  198. }
  199.  
Success #stdin #stdout 0.03s 3776KB
comments (?)
stdin 
Standard input is empty
stdout 
-68290
https://ideone.com/qe1YHA
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!