#ifndef TF_INDICATORS_CURRENCY_MATH_MQH
#define TF_INDICATORS_CURRENCY_MATH_MQH

bool TF_Solve7(double &coefficients[],double &rhs[],double &solution[])
  {
   ArrayResize(solution,7);
   for(int pivot=0;pivot<7;pivot++)
     {
      int best=pivot;
      double best_value=MathAbs(coefficients[pivot*7+pivot]);
      for(int row=pivot+1;row<7;row++)
        {
         double value=MathAbs(coefficients[row*7+pivot]);
         if(value>best_value) { best_value=value; best=row; }
        }
      if(best_value<1.0e-10) return(false);
      if(best!=pivot)
        {
         for(int col=pivot;col<7;col++)
           {
            double tmp=coefficients[pivot*7+col];
            coefficients[pivot*7+col]=coefficients[best*7+col];
            coefficients[best*7+col]=tmp;
           }
         double tmp_rhs=rhs[pivot]; rhs[pivot]=rhs[best]; rhs[best]=tmp_rhs;
        }
      double divisor=coefficients[pivot*7+pivot];
      for(int col=pivot;col<7;col++) coefficients[pivot*7+col]/=divisor;
      rhs[pivot]/=divisor;
      for(int row=0;row<7;row++)
        {
         if(row==pivot) continue;
         double factor=coefficients[row*7+pivot];
         if(MathAbs(factor)<DBL_EPSILON) continue;
         for(int col=pivot;col<7;col++) coefficients[row*7+col]-=factor*coefficients[pivot*7+col];
         rhs[row]-=factor*rhs[pivot];
        }
     }
   for(int i=0;i<7;i++) solution[i]=rhs[i];
   return(true);
  }

bool TF_CurrencyLeastSquares(const int &base_index[],const int &quote_index[],
                             const double &observations[],const double &weights[],
                             double &strengths[],double &residual_rms)
  {
   int count=ArraySize(observations);
   if(count<7 || ArraySize(base_index)!=count || ArraySize(quote_index)!=count || ArraySize(weights)!=count)
      return(false);
   double normal[];
   double rhs[];
   ArrayResize(normal,49);
   ArrayResize(rhs,7);
   ArrayInitialize(normal,0.0);
   ArrayInitialize(rhs,0.0);

   for(int equation=0;equation<count;equation++)
     {
      double x[7];
      ArrayInitialize(x,0.0);
      if(base_index[equation]>=0 && base_index[equation]<7) x[base_index[equation]]+=1.0;
      if(quote_index[equation]>=0 && quote_index[equation]<7) x[quote_index[equation]]-=1.0;
      double weight=MathMax(weights[equation],0.0);
      for(int row=0;row<7;row++)
        {
         rhs[row]+=weight*x[row]*observations[equation];
         for(int col=0;col<7;col++) normal[row*7+col]+=weight*x[row]*x[col];
        }
     }

   double reduced[];
   if(!TF_Solve7(normal,rhs,reduced)) return(false);
   ArrayResize(strengths,8);
   for(int i=0;i<7;i++) strengths[i]=reduced[i];
   strengths[7]=0.0;
   double mean=0.0;
   for(int i=0;i<8;i++) mean+=strengths[i];
   mean/=8.0;
   for(int i=0;i<8;i++) strengths[i]-=mean;

   double weighted_error=0.0,weight_sum=0.0;
   for(int equation=0;equation<count;equation++)
     {
      double fitted=strengths[base_index[equation]]-strengths[quote_index[equation]];
      double error=observations[equation]-fitted;
      double weight=MathMax(weights[equation],0.0);
      weighted_error+=weight*error*error;
      weight_sum+=weight;
     }
   residual_rms=(weight_sum>0.0 ? MathSqrt(weighted_error/weight_sum) : DBL_MAX);
   return(true);
  }

#endif
