# HG changeset patch
# User Jordi Gutiérrez Hermoso <jordigh@octave.org>
# Date 1421631722 18000
# Node ID cf8c8c825c6d8e67c2ac75974c9de1ec9e9ac807
# Parent  65f6fad65f4918b3f12867ab5bb7702c3f6bbf1b
jmedcouple: change all ints to long, to avoid overflow

diff --git a/jmedcouple.c++ b/jmedcouple.c++
--- a/jmedcouple.c++
+++ b/jmedcouple.c++
@@ -24,26 +24,25 @@
   weights W.
 */
 double wmedian(const vector<double>& A,
-               const vector<int>& W)
+               const vector<long>& W)
 {
-  typedef pair<double,int> aw_t;
+  typedef pair<double, long> aw_t;
 
-  int n = A.size();
+  long n = A.size();
   vector<aw_t> AW(n);
-  for (int i = 0; i < n; i ++) 
+  for (long i = 0; i < n; i ++)
     AW[i] = make_pair(A[i], W[i]);
 
   long wtot = sum(W);
 
-  int beg = 0;
-  int end = n - 1;
+  long beg = 0;
+  long end = n - 1;
 
   while (true) {
-    int mid = (beg + end)/2;
+    long mid = (beg + end)/2;
     nth_element(AW.begin(), AW.begin() + mid, AW.end(),
                 [](const aw_t& l, const aw_t& r){return l.first > r.first;});
     double trial = AW[mid].first;
-
     long wleft = 0, wright = 0;
 
     for (aw_t& aw: AW ) {
@@ -69,7 +68,7 @@
                  double eps1 = numeric_limits<double>::epsilon(), 
                  double eps2 = numeric_limits<double>::min())
 {
-  size_t n = X.size(), n2 = (n - 1)/2;
+  long n = X.size(), n2 = (n - 1)/2;
   
   if (n < 3)
     return 0;
@@ -106,7 +105,7 @@
   copy_if(Z.begin(), Z.end(), back_inserter(Zminus), 
           [=](double z){return Zeps >= z;});
 
-  size_t 
+  long
     n_plus = Zplus.size(), 
     n_minus = Zminus.size();
 
@@ -117,7 +116,7 @@
     In case a and be are within epsilon of the median, the kernel
     is the signum of their position.
   */
-  auto h_kern = [&](int i, int j) {
+  auto h_kern = [&](long i, long j) {
     double a = Zplus[i];
     double b = Zminus[j];
 
@@ -131,8 +130,8 @@
   };
  
   // Init left and right borders
-  vector<int> L(n_plus, 0);
-  vector<int> R(n_plus, n_minus - 1);
+  vector<long> L(n_plus, 0);
+  vector<long> R(n_plus, n_minus - 1);
 
   long Ltot = 0;
   long Rtot = n_minus*n_plus;
@@ -140,10 +139,11 @@
   
   // kth pair algorithm (Johnson & Mizoguchi)
   while (Rtot - Ltot > n_plus) {
+
     // First, compute the median inside the given bounds
     vector<double> A;
-    vector<int> W;
-    for (int i = 0; i < n_plus; i++) {
+    vector<long> W;
+    for (long i = 0; i < n_plus; i++) {
       if (L[i] <= R[i]) {
         A.push_back(h_kern(i, (L[i] + R[i])/2));
         W.push_back(R[i] - L[i] + 1);
@@ -155,11 +155,11 @@
 
     //Compute new left and right boundaries, based on the weighted
     //median
-    vector<int> P(n_plus), Q(n_plus);
+    vector<long> P(n_plus), Q(n_plus);
 
     {
-      int j = 0;
-      for (int i = n_plus - 1; i >= 0; i--) {
+      long j = 0;
+      for (long i = n_plus - 1; i >= 0; i--) {
         while (j < n_minus and h_kern(i, j) - Am > Am_eps)
           j++;
         P[i] = j - 1;
@@ -167,8 +167,8 @@
     }
 
     {
-      int j = n_minus - 1;
-      for (int i = 0; i < n_plus; i++) {
+      long j = n_minus - 1;
+      for (long i = 0; i < n_plus; i++) {
         while (j >= 0 and h_kern(i, j) - Am < -Am_eps)
           j--;
         Q[i] = j + 1;
@@ -197,8 +197,8 @@
   // to find it in, just between the left and right boundaries. This
   // space is of size Rtot - Ltot which is <= n_plus
   vector<double> A;
-  for (int i = 0; i < n_plus; i++) {
-    for (int j = L[i]; j <= R[i]; j++)
+  for (long i = 0; i < n_plus; i++) {
+    for (long j = L[i]; j <= R[i]; j++)
       A.push_back(h_kern(i, j));
   }
   nth_element(A.begin(), A.begin() + medc_idx - Ltot, A.end(),