void SortLines(long *lin, long N, double **d)
  {
    if (N <= 2) return;
    /* Find a diametral pair of points, put it in lin[0], lin[1]: */
    { long i, j, iMax, jMax, k;
      double dMax = 0;
      for (i=0; i<N; i++)
        for (j=0; j<N; j++)
          { double dd = d[lin[i]][lin[j]]; 
            if (dd > dMax) { dMax = dd; iMax = i; jMax = j; }
          }
      k = lin[0]; lin[0] = lin[iMax]; lin[iMax] = k;
      k = lin[1]; lin[1] = lin[jMax]; lin[jMax] = k;
    }
    /* Now insert all other points along the path: */
    { long i, r;
      r = 0;
      for (i=2; i<N; i++)
        { long ki, ka, kb, j, jMin;
          double da, db, c, cMin;
          /* Advance "r" by "N/i", cycling in the range "[0..N-i-1]": */
          r += N/i;
          while (r>=N-i) { r -= N-i; }
          /* Bring "lin[i+r]" to "lin[i]": */
          if (r!=0) { long k = lin[i+r]; lin[i+r] = lin[i]; lin[i] = k; }
          /* Insert "lin[i]" in the path "lin[0..i-1]". */
          /* Find the pair "lin[jMin-1], lin[jMin]" such that */
          /* inserting "lin[i]" between them will do least damage: */
          ki = lin[i];
          kb = lin[0]; 
          db = d[ki][kb];
          /* Consider inserting "lin[i]" before "lin[0]": */
          cMin = db; jMin = 0;
          for (j=1; j<i; j++)
            { /* Consider inserting "lin[i]" between "lin[j-1]" and "lin[j]": */
              ka = kb; da = db; 
              kb = lin[j]; db = d[ki][kb];
              c = da + db - d[ka][kb];
              if (c < cMin) { cMin = c; jMin = j; }
            }
          /* Consider inserting "lin[i]" after "lin[i]": */
          if (db < cMin) { cMin = db; jMin = i; }
          /* Now insert "lin[i]" at "lin[jMin]": */
          if (jMin < i)
            { /* Open space for "lin[i]": */
              for (j=i; j>jMin; j--) { lin[j] = lin[j-1]; }
              lin[jMin] = ki;
            }
        }
    }
  }