//+------------------------------------------------------------------+
//|                                                        cmath.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
//+------------------------------------------------------------------+
//| utilities for working with complex numbers                       |
//+------------------------------------------------------------------+

//+------------------------------------------------------------------+
//| get real part from container of complex numbers                  |
//+------------------------------------------------------------------+
void real(complex& data[],double& out[])
  {
   // Resize the output array to match the source array size
   ArrayResize(out,data.Size());
   // Extract the real component element-by-element
   for(ulong i = 0; i<out.Size(); i++)
      out[i] = data[i].real;
  }

//+------------------------------------------------------------------+
//| get imaginary part from container of complex numbers             |
//+------------------------------------------------------------------+
void imaginary(complex& data[],double& out[])
  {
   // Resize the output array to match the source array size
   ArrayResize(out,data.Size());
   // Extract the imaginary component element-by-element
   for(ulong i = 0; i<out.Size(); i++)
      out[i] = data[i].imag;
  }

//+------------------------------------------------------------------+
//| get real part from container of complex numbers                  |
//+------------------------------------------------------------------+
vector real(vectorc& data)
  {
   // Initialize an output vector of the same size with zeros
   vector out=vector::Zeros(data.Size());
   for(ulong i = 0; i<out.Size(); i++)
      out[i] = data[i].real;
   return out;
  }

//+------------------------------------------------------------------+
//| get imaginary part from container of complex numbers             |
//+------------------------------------------------------------------+
vector imaginary(vectorc& data)
  {
   // Initialize an output vector of the same size with zeros
   vector out=vector::Zeros(data.Size());
   for(ulong i = 0; i<out.Size(); i++)
      out[i] = data[i].imag;
   return out;
  }

//+------------------------------------------------------------------+
//| get real part from container of complex numbers                  |
//+------------------------------------------------------------------+
matrix real(matrixc& data)
  {
   // Initialize an output matrix matching dimensions with zeros
   matrix out=matrix::Zeros(data.Rows(),data.Cols());
   // Iterate through rows and columns to copy real values
   for(ulong i = 0; i<out.Rows(); i++)
      for(ulong j = 0; j<out.Cols(); j++)
         out[i,j] = data[i,j].real;
   return out;
  }

//+------------------------------------------------------------------+
//| get imaginary part from container of complex numbers             |
//+------------------------------------------------------------------+
matrix imaginary(matrixc& data)
  {
   // Initialize an output matrix matching dimensions with zeros
   matrix out=matrix::Zeros(data.Rows(),data.Cols());
   // Iterate through rows and columns to copy imaginary values
   for(ulong i = 0; i<out.Rows(); i++)
      for(ulong j = 0; j<out.Cols(); j++)
         out[i,j] = data[i,j].imag;
   return out;
  }

//+------------------------------------------------------------------+
//| create complex containers from real numbers                      |
//+------------------------------------------------------------------+
vectorc as_complex(vector& re,vector& im)
  {
   // Size output based on the largest input vector to prevent truncation
   vectorc out = vectorc::Zeros(MathMax(re.Size(),im.Size()));
   for(ulong i = 0; i<out.Size(); i++)
     {
      // Fallback to 0.0 if index exceeds the size of an uneven input vector
      out[i].real = i<re.Size()?re[i]:0.0;
      out[i].imag = i<im.Size()?im[i]:0.0;
     }
   return out;
  }

//+------------------------------------------------------------------+
//| create complex containers from real numbers                      |
//+------------------------------------------------------------------+
matrixc as_complex(matrix& re,matrix& im)
  {
   // Size output based on the maximum bounds of both matrices
   matrixc out = matrixc::Zeros((ulong)MathMax(re.Rows(),im.Rows()),(ulong)MathMax(re.Cols(),im.Cols()));
   for(ulong i = 0; i<out.Rows(); i++)
     {
      for(ulong j = 0; j<out.Cols(); j++)
        {
         // Safely extract real and imaginary values within bound checks
         out[i,j].real = (i<re.Rows() && j<re.Cols())?re[i,j]:0.0;
         out[i,j].imag = (i<im.Rows() && j<im.Cols())?im[i,j]:0.0;
        }
     }
   return out;
  }

//+------------------------------------------------------------------+
//| create complex containers from real numbers                      |
//+------------------------------------------------------------------+
vectorc as_complex(vector& re)
  {
   vectorc out = vectorc::Zeros(re.Size());
   for(ulong i = 0; i<out.Size(); i++)
     {
      out[i].real = re[i];
      out[i].imag = 0.0; // Imaginary part defaults to zero
     }
   return out;
  }

//+------------------------------------------------------------------+
//| create complex containers from real numbers                      |
//+------------------------------------------------------------------+
matrixc as_complex(matrix& re)
  {
   matrixc out = matrixc::Zeros(re.Rows(),re.Cols());
   for(ulong i = 0; i<out.Rows(); i++)
     {
      for(ulong j = 0; j<out.Cols(); j++)
        {
         out[i,j].real = re[i,j];
         out[i,j].imag = 0.0; // Imaginary part defaults to zero
        }
     }
   return out;
  }

//+------------------------------------------------------------------+
//|  struct for complex vectors                                      |
//+------------------------------------------------------------------+
struct CComplexVector
  {
private:
   // Unpack vectorc structured data into separate real and imaginary vectors
   void              to_real(vectorc& data)
     {
      im = re = vector::Zeros(data.Size());
      for(ulong i = 0; i<data.Size(); i++)
        {
         re[i] = data[i].real;
         im[i] = data[i].imag;
        }
     }

   // Convert purely real vector input by initializing the imaginary vector to zero
   void              to_complex(vector& in)
     {
      im = vector::Zeros(in.Size());
      re = in;
     }

   // Unpack array of complex structures into internal real/imaginary vectors
   void              to_real(complex & data[])
     {
      im = re = vector::Zeros(data.Size());
      for(uint i = 0; i<data.Size(); i++)
        {
         re[i] = data[i].real;
         im[i] = data[i].imag;
        }
     }

   // Convert double array input by initializing internal vectors and assigning values
   void              to_complex(double& in[])
     {
      re = im = vector::Zeros(in.Size());
      re.Assign(in);
     }

   // Reset vector allocations to zero size
   void              reset(void)
     {
      re = vector::Zeros(0);
      im = vector::Zeros(0);
     }
   
public:
   vector            re; // Holds real parts
   vector            im; // Holds imaginary parts

   // Default Constructor
   CComplexVector(void)
     {
      re = vector::Zeros(0);
      im = vector::Zeros(0);
     }

   // Overloaded Constructors
   CComplexVector(vectorc& complx) { to_real(complx); }
   CComplexVector(vector& realv)    { to_complex(realv); }
   CComplexVector(CComplexVector& other)
     {
      re = other.re;
      im = other.im;
     }

   // Overloaded index operator to retrieve elements as a singular complex structure
   complex           operator[](ulong index)
     {
      complex out=0;
      if(index<re.Size())
         {
          out.real = re[index];
          out.imag = im[index];
          return out;
         }
      else
        {
         Print(__FUNCTION__, " invalid index ", index, " : container size ", re.Size());
         return out;
        }
     }

   // Assignment Operators
   void              operator=(vectorc& complex_vector)  { to_real(complex_vector); }
   void              operator=(vector& real_vector)      { to_complex(real_vector); }
   void              operator=(complex& complex_array[]) { to_real(complex_array); }
   void              operator=(double& real_array[])     { to_complex(real_array); }

   // Convert back to native vectorc format
   vectorc           to_vector(void)                     { return as_complex(re,im); }
  };

//+------------------------------------------------------------------+
//|  struct for complex matrices                                     |
//+------------------------------------------------------------------+
struct CComplexMatrix
  {
private:
   // Unpack matrixc structure into separate internal real and imaginary matrices
   void              to_real(matrixc& data)
     {
      im = re = matrix::Zeros(data.Rows(),data.Cols());
      for(ulong i = 0; i<data.Rows(); i++)
        {
         for(ulong j = 0; j<data.Cols(); j++)
           {
            re[i,j] = data[i][j].real;
            im[i,j] = data[i][j].imag;
           }
        }
     }

   // Convert purely real matrix input by zeroing out the imaginary matrix
   void              to_complex(matrix& in)
     {
      im = matrix::Zeros(in.Rows(),in.Cols());
      re = in;
     }

   // Clear memory footprint of internal matrices
   void              reset(void)
     {
      re = matrix::Zeros(0,0);
      im = matrix::Zeros(0,0);
     }
   
public:
   matrix            re; // Holds real matrix space
   matrix            im; // Holds imaginary matrix space

   // Default Constructor
   CComplexMatrix(void)
     {
      re = matrix::Zeros(0,0);
      im = matrix::Zeros(0,0);
     }

   // Overloaded Constructors
   CComplexMatrix(matrixc& complx)       { to_real(complx); }
   CComplexMatrix(matrix& realm)         { to_complex(realm); }
   CComplexMatrix(CComplexMatrix& other)
     {
      re = other.re;
      im = other.im;
     }

   // Safe element retrieval at explicit row/column coordinates
   complex           At(ulong row_index,ulong column_index)
     {
      complex out = 0;
      if(row_index<re.Rows() && column_index<re.Cols())
         {
          out.real = re[row_index,column_index];
          out.imag = im[row_index,column_index];
          return out;
         }
      else
        {
         Print(__FUNCTION__, " invalid row or column index ", " container rows ", re.Rows(), " container columns ", re.Cols());
         return complex(0.0);
        }
     }

   // Assignment Operators
   void              operator=(matrixc& complex_matrix) { to_real(complex_matrix); }
   void              operator=(matrix& real_matrix)     { to_complex(real_matrix); }

   // Convert back to native matrixc format
   matrixc           to_matrix(void)                    { return as_complex(re,im); }
  };

//+------------------------------------------------------------------+
//|         power of complex number                                  |
//+------------------------------------------------------------------+
complex c_pow(complex base, complex exponent)
  {
   // Calculated via identity: base^exponent = exp(exponent * ln(base))
   return c_exp(exponent*c_log(base));
  }

//+------------------------------------------------------------------+
//| simple power of complex number                                   |
//+------------------------------------------------------------------+
vectorc c_pow(vectorc &vect, complex exponent)
  {
   vectorc out(vect.Size());
   for(ulong i = 0; i<vect.Size(); i++)
      out[i] = c_pow(vect[i],exponent);
   return out;
  }

//+------------------------------------------------------------------+
//| simple power of complex number                                   |
//+------------------------------------------------------------------+
matrixc c_pow(matrixc &matrx, complex exponent)
  {
   matrixc out(matrx.Rows(),matrx.Cols());
   // Process column by column using vector-based power calculations
   for(ulong i = 0; i<out.Cols(); i++)
      out.Col(c_pow(matrx.Col(i),exponent),i);
   return out;
  }

//+------------------------------------------------------------------+
//|         power of complex number                                  |
//+------------------------------------------------------------------+
complex c_pow(complex value, double alpha)
  {
   complex out = value;
   double r = c_abs(value);                         // Magnitude
   double theta = atan2(value.imag,value.real);     // Phase Angle
   double n_r = pow(r,alpha);                       // New Magnitude
   double n_theta = alpha*theta;                    // New Phase Angle

   // Convert back from polar form to rectangular form
   out.real = n_r*cos(n_theta);
   out.imag = n_r*sin(n_theta);
   return out;
  }

//+------------------------------------------------------------------+
//| simple power of complex number                                   |
//+------------------------------------------------------------------+
vectorc c_pow(vectorc &vect, double alpha)
  {
   vectorc out(vect.Size());
   for(ulong i = 0; i<vect.Size(); i++)
      out[i] = c_pow(vect[i],alpha);
   return out;
  }

//+------------------------------------------------------------------+
//| simple power of complex number                                   |
//+------------------------------------------------------------------+
matrixc c_pow(matrixc &matrx, double alpha)
  {
   matrixc out(matrx.Rows(),matrx.Cols());
   // Process matrix columns individually
   for(ulong i = 0; i<out.Cols(); i++)
      out.Col(c_pow(matrx.Col(i),alpha),i);
   return out;
  }

//+------------------------------------------------------------------+
//|  Complex Signum Function                                         |
//+------------------------------------------------------------------+
complex signum(complex z)
  {
   double magnitude = sqrt((z.real*z.real) + (z.imag*z.imag));
   if(!magnitude)
      return 0+0i; // Handle division by zero for origin coordinates
   
   complex magz;
   magz.real = magnitude;
   magz.imag = 0.0;

   // Returns unit vector directed in the same phase angle (z / |z|)
   return z/magz;
  }

//+------------------------------------------------------------------+
//| absolute value of a complex number (magnitude)                   |
//+------------------------------------------------------------------+
double c_abs(complex z)
  {
   return sqrt((z.real*z.real) + (z.imag*z.imag));
  }

//+------------------------------------------------------------------+
//| absolute value of a complex number (magnitude)                   |
//+------------------------------------------------------------------+
vector c_abs(vectorc &z)
  {
   CComplexVector v(z);
   // Vectorized calculation of Euclidean length
   return sqrt((v.re*v.re) + (v.im*v.im));
  }

//+------------------------------------------------------------------+
//| absolute value of a complex number (magnitude)                   |
//+------------------------------------------------------------------+
matrix abs(matrixc &z)
  {
   CComplexMatrix v(z);
   // Matrix-wide vectorized calculation of absolute values
   return sqrt((v.re*v.re) + (v.im*v.im));
  }

//+------------------------------------------------------------------+
//| sqrt of complex number                                           |
//+------------------------------------------------------------------+
complex c_sqrt(complex z)
  {
   double mag_z = c_abs(z);
   mag_z=sqrt(mag_z);                           // Square root of the magnitude
   double theta = atan2(z.imag,z.real);         // Phase Angle
   complex out;
   // Half-angle principles for complex square root calculation
   out.real = mag_z*cos(theta/2.0);
   out.imag = mag_z*sin(theta/2.0);
   return out;
  }

//+------------------------------------------------------------------+
//| sqrt of complex vector                                           |
//+------------------------------------------------------------------+
vectorc c_sqrt(vectorc &z)
  {
   vectorc out = z;
   for(ulong i = 0; i<out.Size(); i++)
      out[i] = c_sqrt(z[i]);
   return out;
  }

//+------------------------------------------------------------------+
//| sqrt of complex vector                                           |
//+------------------------------------------------------------------+
matrixc c_sqrt(matrixc &z)
  {
   matrixc out = z;
   for(ulong i = 0; i<out.Cols(); i++)
      out.Col(c_sqrt(z.Col(i)),i);
   return out;
  }

//+------------------------------------------------------------------+
//|  log of complex number                                           |
//+------------------------------------------------------------------+
complex c_log(complex z)
  {
   complex out;
   // Natural logarithm of complex number: ln(z) = ln(|z|) + i*theta
   out.real = c_abs(z);
   out.imag = atan2(z.imag,z.real);
   return out;
  }

//+------------------------------------------------------------------+
//|  log of complex number                                           |
//+------------------------------------------------------------------+
vectorc c_log(vectorc &z)
  {
   vectorc out = z;
   for(ulong i = 0; i<out.Size(); i++)
      out[i] = c_log(z[i]);
   return out;
  }

//+------------------------------------------------------------------+
//|  log of complex number                                           |
//+------------------------------------------------------------------+
matrixc c_log(matrixc &z)
  {
   matrixc out = z;
   for(ulong i = 0; i<out.Cols(); i++)
      out.Col(c_log(z.Col(i)),i);
   return out;
  }

//+------------------------------------------------------------------+
//|  exp function of complex number                                  |
//+------------------------------------------------------------------+
complex c_exp(complex z)
  {
   complex out;
   // Euler's Formula conversion: exp(x + iy) = exp(x) * (cos(y) + i*sin(y))
   out.real = exp(z.real)*cos(z.imag);
   out.imag = exp(z.real)*sin(z.imag);
   return out;
  }

//+------------------------------------------------------------------+
//|  exp function of complex number                                  |
//+------------------------------------------------------------------+
vectorc c_exp(vectorc &z)
  {
   vectorc out = z;
   for(ulong i = 0; i<out.Size(); i++)
      out[i] = c_exp(z[i]);
   return out;
  }

//+------------------------------------------------------------------+
//|  exp function of complex number                                  |
//+------------------------------------------------------------------+
matrixc c_exp(matrixc &z)
  {
   matrixc out = z;
   for(ulong i = 0; i<out.Cols(); i++)
      out.Col(c_exp(z.Col(i)),i);
   return out;
  }

//+-----------------------------------------------------------------------+
//| Generate a complex Vandermonde matrix                                 |
//| in : vectorc                                                          |
//| num_columns : ulong,  default: size of in input vector.               |
//| ascending : bool, (default: false) Order of the powers of the columns.|
//| If true, the powers increase                                          |
//| from left to right, if false (the default) they are reversed          |
//+-----------------------------------------------------------------------+
matrixc vanderc(vectorc &in, ulong num_columns = 0,bool ascending = false)
  {
   ulong n = num_columns?num_columns:in.Size();

   matrixc out(in.Size(),n);

   // Populate columns by geometric power progression rules specified by 'ascending'
   for(ulong i = 0; i<n; i++)
      out.Col(c_pow(in,ascending?double(i):double(n-1-i)),i);

   return out;
  }

//+------------------------------------------------------------------+
//|process vector representing complex numbers in compact form       |
//+------------------------------------------------------------------+
vectorc compact_to_complex(vector& v,vectorc& bp)
  {
   vectorc v_out = vectorc::Zeros(v.Size());

   for(ulong i = 0; i<v.Size() && i<bp.Size();)
     {
      // If the blueprint tracking parameter has an imaginary element, parse conjugate pairs
      if(bp[i].imag)
        {
         v_out[i].real = v[i];
         if(i+1 < v.Size())
           {
            // Maps the pair out as real + i*imag and real - i*imag sequential indexes
            v_out[i].imag = v[i+1];
            v_out[i+1].real = v_out[i].real;
            v_out[i+1].imag = -1.0*v_out[i].imag; // Conjugate inversion
            i+=2; // Step over both elements processed
           }
        }
      else
        {
         // Purely real number fallback step
         v_out[i].real = v[i];
         i+=1;
        }
     }

   return v_out;
  }

//+------------------------------------------------------------------+
//| process matrix representing complex numbers in compact form      |
//+------------------------------------------------------------------+
matrixc compact_to_complex(matrix& m, vectorc& bp)
  {
   matrixc m_out = matrixc::Zeros(m.Rows(), m.Cols());

   // Row-by-row matrix translation leveraging the underlying vector unpacker
   for(ulong i = 0; i<m.Rows(); i++)
     {
      vector row = m.Row(i);
      vectorc row_c = compact_to_complex(row,bp);
      m_out.Row(row_c,i);
     }
   return m_out;
  }
//+------------------------------------------------------------------+