for (int32_t r = 0; r < NX; r++)
          { for (int32_t s = 0; s < NX; s++) 
              { fprintf(stderr, " %+18.9f", A[r*NX + s]); }
            fprintf(stderr, "  ");
            fprintf(stderr, " %+18.9f", B[r]);
            fprintf(stderr, "\n");
          }
        /* Print the Least Squares system: */
        for (int32_t r = 0; r < NX; r++)
          { for (int32_t s = 0; s < NX; s++) 
              { fprintf(stderr, " %+18.9f", A[r*NX + s]); }
            fprintf(stderr, "  ");
            fprintf(stderr, " %+18.9f", B[r]);
            fprintf(stderr, "\n");
          }
    
    if (verbose) { fprintf(stderr, "inverting the moment matrix...\n"); }
    double *M = rmxn_alloc(NX,NX);
    (void)rmxn_inv(NX, A, M);
    if (verbose) { rmxn_gen_print(stderr, NX, NX, M, "%+18.9f", "  ","\n  ","\n", "[ "," "," ]"); }

    if (verbose) { fprintf(stderr, "solving system...\n"); }
    rmxn_map_col(NX, NX, M, B, C);
    if (verbose) { rmxn_gen_print(stderr, NX, 1, C, "%+18.9f", "  ","\n  ","\n", "[ "," "," ]"); }

 
void lif_remove_averages
  ( int32_t NZ, 
    double Z[], 
    double W[], 
    int32_t NX, 
    double X[], 
    double *Zave_P,
    double Xave[],
    bool_t verbose
  );
  /*  Computes the weighted average {Zave} of the values {Z[0..NZ-1]},
    and the weighted average {Xave[k]} of {X[0..NZ-1][k]}, for each {k} in {0..NX-1]},
    all with weights {W[0..NZ-1]}. Then subtracts those averages 
    from {Z[0..NZ-1]} and {X[0..NZ-1][k]}.  Returns {Zave} in {*Zave_P}. */

    double *Z_save = NULL; /* {Z[0..NZ-1]} before shifting to zero average. */  
    double Zave = NAN; /* Weighted average value of {Z[0..NZ-1]}. */
    double Xave[NX]; /* {Xave[kx]} is weighted average value of {X[0..NZ-1][kx]}. */
    if (unitTerm)
      { /* Save original {Z[0..NZ-1]}: */
        Z_save = rn_alloc(NZ); /* {Z[0..NZ-1]} shifted to zero average. */
        for (int32_t i = 0; i < NZ; i++) { Z_save[i] = Z[i]; }
        lif_remove_averages(NZ, Z, W, NX, X, &Zave, Xave, verbose);
      }
  
void lif_remove_averages
  ( int32_t NZ, 
    double Z[], 
    double W[], 
    int32_t NX, 
    double X[], 
    double *Zave_P,
    double Xave[],
    bool_t verbose
  )
  {
    if (verbose) { fprintf(stderr, "removing variable averages...\n"); }
    /* Compute the weighted sums of {X[k],Z} in {Xave[k],Zave}: */
    for (int32_t k = 0; k < NX; k++) { Xave[k] = 0; }
    double Zave = 0;
    double sumW = 0;
    for (int32_t i = 0; i < NZ; i++)
      { double Wi = W[i];
        demand(Wi >= 0, "invalid weight {W[i]}");
        double *Xi = &(X[i*NX]); 
        for (int32_t k = 0; k < NX; k++) { Xave[k] += Wi*Xi[k]; }
        Zave += Wi*Z[i];
        sumW += Wi;
      }
    /* Convert weighted sums to weighted averages: */
    if (sumW > 0)
      { for (int32_t k = 0; k < NX; k++) 
          { Xave[k] /= sumW;
            if (verbose) { fprintf(stderr, "    Xave[%2d] = %14.8f\n", k, Xave[k]); }
          }
        Zave /= sumW;
        if (verbose) { fprintf(stderr, "    Zave     = %14.8f\n", Zave); }
      }
    /* Subtract averages from all variables: */
    for (int32_t i = 0; i < NZ; i++)
      { double *Xi = &(X[i*NX]); 
        for (int32_t k = 0; k < NX; k++) { Xi[k] -= Xave[k]; }
        Z[i] -= Zave;
      }
      
    (*Zave_P) = Zave;
  }

double lif_compute_constant_term(double Zave, int32_t NX, double Xave[], double C[], bool_t verbose);
  /* Assumes that {C[0..NX+1]} are the coefficients of the linear regression of 
    dependent variable {Z} as a function of independent variables {X[0..NX-1]},
    after all variables had theit averages {Zave} and {Xave[0..NX-1]} subtracted. 
    Computes the constant term of the affine regression of the original variables. */
    
    if (unitTerm)
      { /* Compure the constant term {C[NX]}: */
        assert(NC == NX + 1);
        C[NX] = lif_compute_constant_term(Zave, NX, Xave, C, verbose);
        /* Restore the original {Z[0..NZ-1]} data: */
        for (int32_t i = 0; i < NZ; i++) { Z[i] = Z_save[i]; }
        free(Z_save); 
        Z_save = NULL;
      }

double lif_compute_constant_term(double Zave, int32_t NX, double Xave[], double C[], bool_t verbose)
  { if (verbose) { fprintf(stderr, "Computing constant term...\n"); }
    double Cunit = Zave;
    for (int32_t k = 0; k < NX; k++) { Cunit -= C[k]*Xave[k]; }
    return Cunit;
  }