//---------------------------------------------------------------------------
// N.V.Shokhirev
// created: 20040809
// modified: 20040809
// modified: 20060103 - switched to real
//---------------------------------------------------------------------------
#ifndef dynarrutilsH
#define dynarrutilsH
#include <math.h>
#include <stdlib.h>
#include "MatUtils.h"
#include "dynarrays.h"
//---------------------------------------------------------------------------

bool SameLim(Lim1D& a1, Lim1D& a2);

bool SameLim(Lim2D& a1, Lim2D& a2);

bool IsSquare(Lim2D& a);

// Diff = max(abs(a1[i] - a2[i]))
real Diff(FArr1D& a1, FArr1D& a2);

// Diff = max(abs(a1[i,j] - a2[i,j]))
real Diff(FArr2D& a1, FArr2D& a2);

// minval = min(a[i,j])
real minval(FArr1D& a, int& i);

// maxval = max(a[i,j])
real maxval(FArr1D& a, int& i);

// maxval = max(abs(a[i,j]))
real maxabs(FArr1D& a, int& i);

// if reset then amin = min(a[i]); amax = max(a[i])
//          else amin = min(amin, a[i]); amax = max(amax, a[i])
void minmax(FArr1D& a, real& amin, real& amax, bool reset);

bool Equal(real x1, real x2, real eps = eps0);

bool Equal(IArr1D& a1, IArr1D& a2);

bool Equal(FArr1D& a1, FArr1D& a2, real eps = eps0); // return (Diff(a1, a2) < eps);

bool Equal(FArr2D& a1, FArr2D& a2, real eps = eps0); // return (Diff(a1, a2) < eps);

// random vector: v1 <= a[i,j] <= v2
void RandArr(IArr1D& a, int v1, int v2);

// random vector: v1 <= a[i,j] <= v2
void RandArr(FArr1D& a, real v1, real v2);

/* random matrix: v1 <= a[i,j] <= v2
sym = 0 - no symmetry
for square a: if sym > 0 then a is symmetric if sym < 0 then a is antisymmetric */
void RandArr(FArr2D& a, real v1, real v2, int sym = 0);

// a[i,i] = a[i,i] + v
void ShiftDiag(FArr2D& a, real v);

// a = v*a
void MultArr(FArr1D& a, real v);

// a = v*a
void MultArr(FArr2D& a, real v);

// Norm2 = (a1*a1)
real Norm2(FArr1D& a1);

// Dot = (a1*a2) - dot or scalar product
real Dot(FArr1D& a1, FArr1D& a2);

// Combine = v1*a1 + v2*a2
FArr1D Combine(real v1, FArr1D& a1, real v2, FArr1D& a2);

// Combine = v1*a1 + v2*a2
FArr2D Combine(real v1, FArr2D& a1, real v2, FArr2D& a2);

// MT = Transpose( a )
FArr2D MT(FArr2D& a);

// MxV = a1 * a2 - vector
FArr1D MxV(FArr2D& a1, FArr1D& a2);

// VxM = a1 * a2 - vector
FArr1D VxM(FArr1D& a1, FArr2D& a2);

// MxM = a1 * a2 - matrix
FArr2D MxM(FArr2D& a1, FArr2D& a2);

/* Fourier transform T0[j] corresponds to non-centered frequency        <br>
                         Q[j] = j*2*Pi/(dT*N)   j = 0, 1, 2, .. , N-1   <br>
   Centered frequency: F[j] = j*2*Pi/(dT*N)  Mmin <= j <= Mmax          <br>
   where Mmin = -((N-1) div 2), Mmax = (M div 2)                        <br>
   The transform corresponding to the centered frequency                <br>
                                    C[m] = T0[ m ]  0 <= m <= Mmax      <br>
                                    C[m] = T0[N-m]  Mmin <= m < 0       <br> */
void CenterToTrans(FArr1D& T, FArr1D& C);

void TransToCenter(FArr1D& T, FArr1D& C);


//  if(i < L) i += D; else if (i > H) i -= D;   return i;
int cycle(int i, int L, int H, int D);

//  if (cyclic) // a(i+D) = a(i-D) = a(i)
//  else 0 outside limits
real extend(FArr1D& a, int i, bool cyclic);

//  if (cyclic) // a[0..N] a(0) = a(N);  a[0] is not used
//  else 0 outside limits
real extend0(FArr1D& a, int i, bool cyclic);

//  5-point 2nd degree smoothing; p1 = p0 - OK
void Smooth5(FArr1D& p0, FArr1D& p1);


#endif