• It finds the smallest n pair sums.
  
• The solution is returned as a pivot value, less than which there are n pair sums.
  
• The pivot is found by binary search between A[0]+B[0] and A[sqrt(n)] +B[sqrt(n)].  Each possible pivot is evaluated by a counting function.
  
• The counting function returns as soon as more than n pairs are found below the pivot.
  
• The counting function starts on the diagonal and works left and right columnwise.  Since there are only sqrt(n) plausible starting points there, this option seemed better than starting on an axis, but I haven't tested this assertion.
  

#include <iostream>

#include "stdlib.h"
#include "math.h"

// 

int binsearch_le( int* X, int n, int v )
{
int l = 0, u = n-1, tmp = 0;

if( X[u] < v ) return u;

while( X[l] != v && l < u-1 )

{
//
std::cout << "l: " << l << " u:" << u << "\n";

tmp = (l+u)>>1;

if( X[tmp] > v )


u = tmp;

else if( X[tmp] < v )


l = tmp;

else


return tmp;

};

return l;
}

// counter for pair tests in cnt; for testing 
int paircount;
#define COUNTIT {paircount++;}
#define RESETCOUNT {paircount = 0;}
#define OUTPUTCOUNT {std::cout << "\npaircount: " << paircount << "\n"; }

// the main operation to count:  add 2 elements and compare to pivot
int pairtest( int a, int b, int v )
{
COUNTIT;
return a+b <= v;
}

// count pair sums <= val
// d is start pos on diagonal, max <= val
// return count as soon as > n
int cnt( int val, int *A, int *B, int n, int d )
{
int i,j;
int c = 0;
// left

for( i = d, j = d; i >=0 && c <= n; i-- )

{

while( j+1 < n && pairtest(A[i], B[j+1], val)  )  // peek ahead


j++ ;

c += j+1;
// j is highest <= in column

};

for( i = d+1, j = d; j >= 0 && i < n && c <= n; i++ )

{

while( j >= 0 && !pairtest(A[i], B[j], val) )


j--;

c += j+1;

}

return c;
}

// returns pivot for n pair sums < pivot
int bottomn( int *A, int *B, int n  )
{

int i,j, a, b;

int l, u, lsum, usum, newsum, sqrtn, count =0;
int * diags = 0;

RESETCOUNT;

l = 0;
u = sqrtn = lround(ceil(sqrt(n)));

// pre-calc sums along diagonal up to sqrtn
diags = new int[sqrtn];
for( i = 0; i < sqrtn; i ++ )

diags[i] = A[i]+B[i];

lsum = A[ l ] + B[ l ];
usum = A[ u ] + B[ u ];

// loop:
// interpolate new pivot value in (lsum, usum) (note quadratic growth, d-count/d-pivot = C*pivot)
// find nearest (<=) diagonal pivot position
// find count outward from diag. pivot
// set lsum or usum to pivot based on count </> n
// stop when count = n (or smaller pivot would count < n)
while( lsum < usum -1 )

{

newsum = (lsum +usum) >> 1;

j = binsearch_le( diags, sqrtn, newsum );

count = cnt( newsum, A, B, n, j );
//
std::cout << "lsum: " << lsum << " usum: " << usum <<" newsum:" << newsum << " count: " << count << " j: " << j << "\n";
//
std::cout << "diags[j]: " << diags[j] << "\n";

if( count > n )


usum = newsum;

else if( count < n )


lsum = newsum;

else


break;

};

OUTPUTCOUNT;

return newsum;
}

int qsort_le(const void *a, const void *b)
{
 return *(int *)a - *(int *)b ;
}

int main (int argc, char * const argv[]) 
{
int pivot;
    // insert code here...

int n;
int i=0;
int j=0;
int *A, *B;
int ntrials = 10;

int pairtot=0;

if( argc < 2) return -1;

n = atoi(argv[1]);

A = new int[n];
B = new int[n];

for( j = 0; j < ntrials; j++ )

{

for(i = 0; i < n; i++ )


{


A[i] = random() % 1000000000;


B[i] = random() % 1000000000;


};


qsort( A, n, sizeof( *A), &qsort_le );

qsort( B, n, sizeof( *B), &qsort_le );


pivot = bottomn( A, B, n );

pairtot += paircount; // op counter


int ii, jj;

int checkcnt = 0;

// check result

for( ii = 0; ii < n; ii++ )


for( jj = 0; jj < n && A[ii]+B[jj] <= pivot; jj++ )



checkcnt++;


std::cout << "n: " << n <<" Pivot:" << pivot << " checkcnt: " << checkcnt << "\n";

};

std::cout << "avg paircount: " << pairtot/ntrials << " avg/n: "<< pairtot/(n*ntrials*1.0) << "\n";

return 0;
}

     10    9    5    4    3    2
12   22   21   17   ..   ..   ..
11   21   20   16   ..  
 9   19   18   14   ..
 6   16   15   ..
 1   ..   ..
 0   ..