Quantum Information and Computation library
(QIClib)

• Every functions, classes, or constants in QIClib belongs to the namespace qic.
• First time users may want to have a look on sample programs.
• First time users may want to have a look on frequently asked questions.
• If you discover any bugs or regressions, please report them
• Notes on API additions

• Classes and constants
• Functions
• Classes for quantum discord like features
• Preprocessor macros to tune the library in compile time

• Note: These features depend on NLopt.
• Note: If you want to use old API's for discord like function, define QICLIB_USE_OLD_DISCORD macro before including QIClib header. See preprocessor macros in QIClib for details. Old API's for discord like functions are not maintained anymore.

Classes and constants

• Class for library initialization. Automatically loaded if you use this library.
• Prints current time and date upon program start. After program completion, prints current time and date along with total runtime.
• To disable this feature, define QICLIB_NO_INIT_MESSAGE macro before including QIClib header. See preprocessor macros in QIClib for details.

• Several predefined special matrices and vectors like, basis, projectors, operators etc. POD_TYPE can be either double or float.
• This is a singleton class, can be constructed only once. You can only take const reference of this class using SPM<POD_TYPE>::get_instance() function (see below). For the ease of use, QIClib library automatically instantiates this class for double and float as spm and spmf respectively. Use of spm or spmf is recommended rather than explicit SPM<POD_TYPE>::get_instance().
• spm.S gives field of Pauli matrices as constant complex matrices (const cx_mat). spm.S(0), spm.S(1), spm.S(2), and spm.S(3) are respectively I_2x2, σ_x, σ_y, and σ_z.
  (Similarly for spmf.S).
• spm.basis2 gives field of 2-dimensional basis vectors as constant complex column vectors (const cx_vec). spm.basis2(0, i) and spm.basis2(1, i) are respectively |0> and |1> in i direction (i = 0, 1, 2, 3).
  spm.proj2 gives corresponding projectors as constant complex matrices (const cx_mat).
  (Similarly for spmf.basis2 and spmf.proj2).
• spm.basis3 gives field of 3-dimensional basis vectors constant complex column vectors (const cx_vec). spm.basis3(0, i), spm.basis3(1, i) and spm.basis3(2, i) are respectively |0>, |1> and |2> in i direction (i = 0, 1, 2, 3).
  spm.proj3 gives corresponding projectors constant complex matrices (const cx_mat).
  (Similarly for spmf.basis3 and spmf.proj3).
• spm.bell.phim, spm.bell.phip, spm.bell.psim, and spm.bell.psip give four Bell states as constant real vectors (const vec).
  (Similarly for spmf.bell).
• Note: You cannot take a copy of the whole SPM<POD_TYPE> class, but you can take reference to const. For this use const SPM<POD_TYPE>& SOMENAME = SPM<POD_TYPE>::get_instance();. However, for each member variables, like S or proj3 etc., you can either take a copy or reference to const (See the following example).
  Note: Taking copy or reference to const is not advisable for the member variables. Direct use of member variables is recommended.
• Example:

   
  	  	const SPM<double>& myspm 
     	 	   = SPM<double>::get_instance(); // take const ref of the whole class
  	   
  		const cx_mat& S1 = spm.S(1);   // \sigma_x matrix, take reference
  		cx_mat S2 = spm.S(2);          // \sigma_y matrix, take a copy
  
  		const cx_vec& plus = spm.basis2(0, 1);     // |+> state
  		const cx_mat& proj_plus = spm.proj2(0, 1); // |+><+| projector
  
  		const cx_vec& U = spm.basis3(0, 3); // 1st eigenvector of 
                  //Spin S_z operator in 3-dim
  
  		const cx_mat& M = spm.basis3(1, 3); // 2nd eigenvector of 
                  //Spin S_z operator in 3-dim
  	
  		const cx_mat& M = spm.basis3(2, 3); // 3rd eigenvector of 
                  //Spin S_z operator in 3-dim
  	      

• Several predefined constant quantum gates and generators of quantum gates. POD_TYPE can be either double or float.
• This is a singleton class, can be constructed only once. You can only take const reference of this class using GATES<POD_TYPE>::get_instance() function (see below). For the ease of use, QIClib library automatically instantiates this class for double and float as gates and gatesf respectively. Use of gates or gatesf is recommended rather than explicit GATES<POD_TYPE>::get_instance().
• gates.X, gates.Y, gates.Z, and gates.Had give Pauli X,Y, and Z and Hadamard gates respectively. gates.Y is a constant complex matrix (const cx_mat), where as others are constant real matrices (const mat).
  (Similarly for gatesf.X, gatesf.Y, gatesf.Z, and gatesf.Had).
• gates.CNOT, gates.CZ and gates.swap give controlled NOT, controlled Pauli Z, and swap gates respectively as constant real matrices (const mat).
  (Similarly for gatesf.CNOT, gatesf.CZ, and gatesf.swap).
• gates.Tof and gates.Fred give Toffoli and Fredkin gates respectively as constant real matrices (const mat).
  (Similarly for gatesf.Tof and gatesf.Fred).
• gates.U2(theta, unitV): Generates a member of U(2) group. theta is a double and uniV is a 3-dimensional unit column vector (vec). Output is a fixed size 2x2 complex unitary matrix (cx_mat22) as U = cos(theta/2) I + i sin(theta/2)(unitV.σ).
  (Similarly for gatesf.U2(theta, unitV)).
• gates.PS(phi): Generates a phase-shift gate according to the angle phi (double). Output is a fixed size 2x2 complex unitary matrix (cx_mat22) as U = diag(1, exp(i*phi)).
  (Similarly for gatesf.PS(phi)).
• gates.qft(dim): Generates unitary matrix for quantum Fourier transform of dimension specified by dim, as complex matrix (cx_mat).
  (Similarly for gatesf.qft(dim)).
• Note: You cannot take a copy of the whole GATES<POD_TYPE> class, but you can take reference to const. For this use const GATES<POD_TYPE>& SOMENAME = GATES<POD_TYPE>::get_instance();. However, for each member variables, like Had or CNOT etc., you can either take a copy or reference to const. Member functions like U2 or qft etc. behave like normal functions. (See the following example).
  Note: Taking copy or reference to const is not advisable for the member variables. Direct use of member variables is recommended.
• Example:

   
  	  	const GATES<double>& mygates 
     	 	   = GATES<double>::get_instance(); // take const ref of the whole class
  	   
  		const mat& XX = gates.X;   // \sigma_x matrix, take reference
  		cx_mat YY = gates.Y;       // \sigma_y matrix, take a copy
  
  		cx_mat Q = gates.qft(8); // quantum Fourier transform 
  		                         // for 8 dimension 
  	  

• Class with several predefined static member functions to be used with funcm_sym or funcm_gen.
• func is an alias for Func<double> and funcf is an alias for Func<float>.
• Every member function takes one complex type and returns one complex type.
• List of member functions:
  • func::sin
    func::cos
    func::tan
  • func::asin
    func::acos
    func::atan
  • func::sinh
    func::cosh
    func::tanh
  • func::asinh
    func::acosh
    func::atanh
  • func::sqrt
    func::log
    func::log2
  • func::real
    func::imag
    func::norm

• Predefined random number engines/generators, used in functions like randU etc. This class is thread safe.
• This is a singleton class, can be constructed only once. You can only take reference of this class using RandomDevices::get_thread_local_instance() or RandomDevices::get_instance(), but this is not advisable at all. QIClib library automatically instantiates this class as redvs. Use rdevs rather than explicit RandomDevices::get_thread_local_instance() or RandomDevices::get_instance().
• There are two member variables, rd and rng. rd is an instance of std::random_device, and rng is an instance of Mersenne twister engine, std::mt19937_64 or std::mt19937 depending on whether Armadillo uses 64bit integers or not. See this.
• By default, rng is seeded with the random return value of rd. To seed rng with user defined value, use rdevs.set_seed(YOUR_SEED). To seed rng with the random return value of rd again, use rdevs.set_seed_random().
• Example:

  		std::cauchy_distribution<double> dis(0,1);
  		// Cauchy distribution with parameters 0 and 1
  
  		double a[10], b[10], c[10];
  
  		// seed is random, by default
  		for(int i = 0; i < 10; ++i) {
  		  a[i] = dis(rdevs.rng); // generate 10 random number
  		                         // according to Cauchy distribution 
  		                         // with parameters 0 and 1	   
  		}
  				   
  		rdevs.set_seed(10); // set seed to 10
  		
  		for(int i = 0; i < 10; ++i) {
  		  b[i] = dis(redvs.rng); // generate 10 random number
  		                         // according to Cauchy distribution 
  		                         // with parameters 0 and 1	   
  		}
  						      
  		rdevs.set_seed_random(); // again seed is random
  		
  		for(int i = 0; i < 10; ++i) {
  		  c[i] = dis(rdevs.rng); // generate 10 random number
  		                         // according to Cauchy distribution 
  		                         // with parameters 0 and 1	   
  		}
  			

Functions

• Hermiticity check. Returns true if A is real symmetric (mat) or Hermitian (cx_mat), else false.
• atol and rtol are optional. By default, atol is floating point precision used in QIClib, and rtol is 10 times of that. See preprocessor macros in QIClib for details.

• Unitarity check. Returns true if A is real orthogonal (mat) or unitary (cx_mat), else false.
• atol and rtol are optional. By default, atol is floating point precision used in QIClib, and rtol is 10 times of that. See preprocessor macros in QIClib for details.

• To check whether the the matrix is normal or not. Returns true if A is a normal matrix, else false.
• atol and rtol are optional. By default, atol is floating point precision used in QIClib, and rtol is 10 times of that. See preprocessor macros in QIClib for details.

• To check purity of a row/column vector or density matrix. Returns true if A is a pure state, else false.
• check_norm and tol are optional. By default, check_norm = true and normalization of the vector or matrix will be checked, and tol is floating point precision used in QIClib. See preprocessor macros in QIClib for details.

• To check whether the matrix or vector represents a valid quantum state. Returns true if A is a valid quantum state, else false.
• tol is optional, by default, tol is floating point precision used in QIClib. See preprocessor macros in QIClib for details.
• Note: Normalization will always be checked.

• To check whether the matrix is diagonalizable. Returns true if A is diagonalizable, else false.

• To check whether two matrices are equal in values. Returns true if A and B are element-wise (approximately) equal, else false.
• typecheck is optional. By default, typecheck=false, and element types (complex or real) of A and B need not to be same.
• atol and rtol are optional. By default, atol is floating point precision used in QIClib, and rtol is 10 times of that. See preprocessor macros in QIClib for details.

• To Convert a dense matrix A (Mat) to a sparse (SpMat) matrix with the same element type.
• tol is optional. Any element having absolute value less than tol will be treated as zero. By default, tol is floating point precision used in QIClib. See preprocessor macros in QIClib for details.

• To Convert a sparse matrix A (SpMat) to a dense (Mat) matrix with the same element type.

• Similar to Python's range() function. Useful in range-based for loops. Returns std::vector of arithmetic (floating point or integral types) values, from (and including) start to (and excluding) stop with step-size equal to step.
• start and step are optional, by default, start = 0 and step = 1.
• Example:

  		for(auto&& i : range(10))
  		{
  		// i from 0 to 9
  		//do something
  		}
  		
  		std::vector<int> R1 = range(2, 10, 2); 
  		std::vector<double> R2 = range(10, -1.5, -0.5);
  	      

• Partial trace of a quantum state, where subsys is a uvec containing the indices of to be traced out subsystems.
• A has be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		mat A = randn<mat>(8, 8); // 8x8 matrix
  		A *= A.t(); // make Hermitian
  		A /= trace(A); //normalize
  
  		cx_mat B = randRho(12); // random 12x12 density matrix, see randRho()
  		
  		// trace out 2nd and 3rd party
  		mat A1 = TrX(A, {2, 3}); // all dim are 2, total 3 parties
  		mat A2 = TrX(A, {2, 3}, 2); // same as before
  		
  		// trace out 2nd and 3rd party
  		cx_mat B1 = TrX(B, {2, 3}, {2, 3, 2}); // explicitly written dimensions
  	      

• Partial transpose of a quantum state, where subsys is a uvec containing the indices of to be transposed subsystems.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		mat A = randn<mat>(8, 8); // 8x8 matrix
  		A *= A.t(); // make Hermitian
  		A /= trace(A); //normalize
  
  		cx_mat B = randRho(12); // random 12x12 density matrix, see randRho()
  		
  		// transpose 2nd and 3rd party
  		mat A1 = Tx(A, {2, 3}); // all dim are 2, total 3 parties
  		mat A2 = Tx(A, {2, 3}, 2); // same as before
  
  		// transpose 2nd and 3rd party
  		cx_mat B1 = Tx(B, {2, 3}, {2, 3, 2}); // explicitly written dimensions
  	      

• Permute subsystems of a quantum state, where perm is a uvec containing the permutation of subsystem indices.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in perm start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		mat A = randn<mat>(8, 8); // 8x8 matrix
  		A *= A.t(); // make hermitian
  		A /= trace(A); //normalise
  
  		cx_mat B = randRho(12); // random 12x12 density matrix, see randRho()
  		
  		// permute 2nd and 3rd party
  		mat A1 = sysperm(A, {1, 3, 2}); // all dim are 2, total 3 parties
  		mat A2 = sysperm(A, {1, 3, 2}, 2); // same as before
  		
  		// permute 2nd and 3rd party
  		cx_mat B1 = sysperm(B, {1, 3, 2}, {2, 3, 2}); // explicitly written dimensions
  		                                              // for B1, dim = {2, 2, 3} 
  	      

• Principal square root of real symmetric or Hermitian matrices.
• A has to be a square matrix (mat or cx_mat).
• Due to generality, always returns complex matrices (cx_mat), even for positive real matrices.
• Note: sqrtm_sym is faster than sqrtm_gen. For real symmetric or Hermitian matrices always use sqrtm_sym.
• Example:

  		cx_mat A = randn<cx_mat>(5,5);
  		A *= A.t(); // make hermitian
  		
  		mat B = randn<mat>(5,5);
  		B *= B.t(); // make hermitian
  		
  		cx_mat sqrtA = sqrtm_sym(A);
  
  		cx_mat sqrtB = sqrtm_sym(B);  
    	      

• Principal square root of general square matrices.
• A has to be a square matrix (mat or cx_mat).
• Due to generality, always returns complex matrices, even for positive real matrices.
• Note: sqrtm_sym is faster than sqrtm_gen. For real symmetric or Hermitian matrices always use sqrtm_sym.
• Note: Only works if A is diagonalizable matrix. To make sure, use is_diagonalizable function first.
• Example:

  		cx_mat A = randn<cx_mat>(5,5);
  		
  		mat B = randn<mat>(5,5);
  		
  		cx_mat sqrtA = sqrtm_gen(A);
  
  		cx_mat sqrtB = sqrtm_gen(B);  
    	      

• Matrix power for real symmetric or Hermitian matrices.
• A has to be a square matrix (mat or cx_mat). n can be any type of scaler, like double, complex<double>, int etc.
• Due to generality, always returns complex matrices (if n is not of integer type), even for positive real matrices. If n is of integer type matrix element type will be preserved.
• Note: powm_sym is faster than powm_gen. For real symmetric or Hermitian matrices always use powm_sym.
• Example:

  		mat A = randn<mat>(5, 5);
  		A *= A.t(); // make hermitian
  		
  		mat A_3 = powm_sym(A, 3);
  
  		cx_mat A_3.14 = powm_sym(A, 3.14);  
    	      

• Matrix power for general square matrices.
• A has to be a square matrix (mat or cx_mat). n can be any type of scaler, like double, complex<double>, int etc.
• Due to generality, always returns complex matrices (if n is not of integer type), even for positive real matrices. If n is of integer type matrix element type will be preserved.
• Note: powm_sym is faster than powm_gen. For real symmetric or Hermitian matrices always use powm_sym.
• Note: Only works if A is diagonalizable matrix. To make sure, use is_diagonalizable function first.
• Example:

  		mat A = randn<mat>(5, 5);
  		
  		mat A_3 = powm_gen(A, 3);
  
  		cx_mat A_3.14 = powm_gen(A, 3.14);  
    	      

• Matrix exponential for real symmetric or Hermitian matrices.expm_sym(A,n) calculates exp(n*A).
• A has to be a square real symmetric or Hermitian matrix (mat or cx_mat). n can be any type of scaler, like double, complex<double>, int etc.
• Note: expm_sym is often slower but more accurate than expm_gen or Armadillo's expmat. But for very large matrices it can be faster than the alternatives.
• Example:

  		mat A = randn<mat>(5,  5);
  		A *= A.t(); // make hermitian
  		
  		mat A_e = expm_sym(A);
  
  		complex<double> I (0.0, 1.0);
  		cx_mat AI_e = expm_sym(A, I); // exp(I*A)
  	      

• Matrix exponential for general square matrices.
• A has to be a square matrix (mat or cx_mat).
• Note: expm_gen is basically same as Armadillo's expmat, and only to be used with older versions of Armadillo.
• Example:

  		cx_mat A = randn<cx_mat>(5, 5);
  		
  		cx_mat A_e = expm_gen(A);
  
  		complex<double> I (0.0, 1.0);
  		cx_mat AI_e = expm_gen(I*A);
  	      

• Matrix function for real symmetric or Hermitian matrices.
• A has to be a square real symmetric or Hermitian matrix (mat or cx_mat). function pointer, functor or lambda function has to take one complex type and return one complex type.
• Due to generality, always returns complex matrices.
• If you want to use standard complex functions, use static member functions of Func<type> class.
• Note: funcm_sym is faster than funcm_gen. For real symmetric or Hermitian matrices always use funcm_sym.
• Note: If you want to transform the matrix element-wise, use .transform instead.
• Example:

  		cx_mat A = randn<cx_mat>(5, 5);
  		A *= A.t(); // make hermitian
  		
  		//define a lambda function
  		auto sinm = [](complex<double> a){return std::sin(a);};
  
  		cx_mat A_sin = funcm_sym(A, sinm); // calculate sin(A)
  
  		cx_mat A_sin2 = funcm_sym(A, func::sin); // same as A_sin, 
  		                                         //using Func<type> class
  	      

• Matrix function for general square matrices.
• A has to be a square matrix (mat or cx_mat). function pointer, functor or lambda function has to take one complex type and return one complex type.
• Due to generality, always returns complex matrices.
• If you want to use standard complex functions, use static member functions of Func<type> class.
• Note: funcm_sym is faster than funcm_gen. For real symmetric or Hermitian matrices always use funcm_sym.
• Note: Only works if A is diagonalizable matrix. To make sure, use is_diagonalizable function first.
• Note: If you want to transform the matrix element-wise, use .transform instead.
• Example:

  		cx_mat A = randn<cx_mat>(5, 5);
  		
  		//define a lambda function
  		auto sinm = [](complex<double> a){return std::sin(a);}; 
  
  		cx_mat A_sin = funcm_gen(A, sinm); // calculate sin(A)
  
  		cx_mat A_sin2 = funcm_gen(A, func::sin); // same as A_sin, 
  		                                         // using Func<type> class
  	      

• Tensor product of arbitrary no. of matrices A,B,.... Return type is deduced from the types of input matrices.
• Can also take std::vector and field of matrices.
• Example:

  		cx_mat A = randn<cx_mat>(5, 5); //complex matrix
  		cx_mat B = randn<cx_mat>(4, 4); //complex matrix
  		cx_mat C = randn<cx_mat>(3, 3); //complex matrix
  		mat D = randn<mat>(6, 6);       //real matrix
  		
  		cx_mat T1 = tensor(A, B, C, D); // tensor product of A,B,C,D 
  		                                // return type is deduced as cx_mat
  
  		std::vector<cx_mat> V = {A, B, C};
  		cx_mat T2 = tensor(V); // tensor product of A,B,C
  		
  		cx_mat T3 = tensor({A, B, C}); // same as above  
  	      

• Tensor product power of a matrix.
• Returns A^⊗n .

• Direct sum of arbitrary no. of matrices A,B,.... Return type is deduced from the types of input matrices.
• Can also take std::vector and field of matrices.
• Example:

  		cx_mat A = randn<cx_mat>(5, 5); // complex matrix 
  		cx_mat B = randn<cx_mat>(4, 4); // complex matrix
  		cx_mat C = randn<cx_mat>(3, 3); // complex matrix
  		mat D = randn<mat>(6, 6);       // real matrix
  		
  		cx_mat T1 = dsum(A, B, C, D); // direct sum of A, B, C, D
  		                              // return type is deduced to be cx_mat
  
  		std::vector<cx_mat> V = {A, B, C}; 
  		cx_mat T2 = dsum(V); // direct sum of A, B, C
  		
  		cx_mat T3 = dsum({A, B, C}); // same as above
  	      

• Direct sum power of a matrix.
• Returns A^⊕n .

• Matrix absolute value of a square matrix.
• A has to be a square matrix (mat or cx_mat).
• Returns sqrtm_sym(A.t()*A). Return type is same as that of A.

• Schatten p-norm of a matrix.
• p has to be greater than or equal to 0.

• Minimal purification of a density matrix. Returns vec or cx_vec depending on the type of A.
• tol is optional. By default, tol is equal to the floating point precision used in QIClib. See preprocessor macros in QIClib for details. Eigenvalues (of A) less than tol are treated zero.

• Modified Gram-Schmidt orthogonalization.
• A can be either a matrix (mat or cx_mat) or std::vector/field of vectors (vec or cx_vec).
• normalize is optional. By default, it is true. If it is false then orthogonal vectors will not be normalized.

• To convert density matrix of a pure state to a column vector (vec or cx_vec).
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• If A is a column vector, then it is returned back unchanged.
• If A is not a pure state, then eigenvector corresponding to highest eigenvalue is returned. To make sure, use is_pure function first.

• To convert 2-qubit density matrix from the standard basis to the Hilbert-Schmidt (or Pauli) basis.
• A has to be a square matrix (mat or cx_mat) of 4x4.
• Always returns a 4x4 fixed size real matrix (mat44).

• To convert 2-qubit density matrix from the Hilbert-Schmidt (or Pauli) basis to the standard basis.
• A has to be a real square matrix (mat) of 4x4.
• Always returns a 4x4 fixed size complex matrix (cx_mat44).

• To create multipartite qudit pure state in standard computational basis according to the uvec mask. (See following example).
• ELEM_TYPE is optional, by default, it is double. It specifies the element type of returned column vector. ELEM_TYPE can be real or complex types.
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Example:

  		vec A = mket({0, 0, 1}); // |001> state
  		                         // each subsystem is qubit
  
  		cx_vec B = mket<complex<double> >({0, 0, 1}); // same a above
  		                                              // but complex vector returned
  		
  		vec C = mket({0, 0, 1}, 2); // same a above
  		
  		vec D = mket({0, 2, 1}, {2, 3, 2}); //|021> state
  		                                    // 1st and 3rd party are qubit
  		                                    // 2nd one is qutrit
  
  		cx_vec E = mket<complex<double> >({0, 2, 1},{2, 3, 2}); //same as above
  		                                                        // but complex vector returned
  	      

• To create multipartite qudit projectors in standard computational basis according to the uvec mask (See following example).
• ELEM_TYPE is optional, by default, it is double. It specifies the element type of returned matrix. ELEM_TYPE can be real or complex types.
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Example:

  		mat A = mproj({0, 0, 1}); // |001><001| projector
  		                          // each subsystem is qubit
  
  		cx_mat B = mproj<complex<double> >({0, 0, 1}); // same a above
  		                                               // but complex matrix returned
  		
  		mat C = mproj({0, 0, 1}, 2); // same a above
  		
  		mat D = mproj({0, 2, 1}, {2, 3, 2}); //|021><021| projector
  		                                     // 1st and 3rd party are qubit
  		                                     // 2nd one is qutrit
  
  		cx_mat E = mproj<complex<double> >({0, 2, 1},{2, 3, 2}); //same as above
  		                                                         // but complex matrix returned
  	      

• To generate object with uniform random values in the range specified by the column vector range. By default, the range is in between 0 to 1.
• The first functions generate only one random element of type (real or complex) specified by ELEM_TYPE. By default, TYPE is double.
• The second functions generate a row or column vector of type (with real or complex elements) specified by VEC_TYPE, with n_elem number of elements. By default, VEC_TYPE is vec.
• The third function generates a matrix of type (with real or complex elements) specified by MAT_TYPE, with n_rows number of rows, n_cols number of columns. By default, MAT_TYPE is mat.
• To change the seed of the random number generator, use rdevs.set_seed(YOUR_SEED). To set random seed, use rdevs.set_seed_random(). By default, the random number generator is instantiated with random seed. See RandomDevices class for more details.
• Example:

  		double a = randU(); // random number between 0 and 1
  		
  		complex<double> b = randU<complex<double> >();  // complex random number 
  		                                                // between (0,0) and (1,1)
  		
  		double c = randU({-5, 5}); // random number between -5 and 5
  		
  		cx_vec V = randU<cx_vec>(10); // cx_vec with 10 elements 
  
  		cx_rowvec R = randU<cx_rowvec>(10, {-5, 5}); // cx_rowvec with 10 elements
  
  		mat M1 = randU(10, 10); // 10x10 mat
  
  		mat M2 = randU(10, 10, {-5, 5}); // 10x10 mat  
  
  		cx_mat M3 = randU<cx_mat>(10, 10, {-5, 5}); // 10x10 cx_mat  
  	      

• To generate object with random values from normal distribution with mean and standard deviation specified by the column vector mean_sd. By default, the mean and standard deviation are 0 and 1 respectively.
• The first function generate only one random element of type (real or complex) specified by ELEM_TYPE. By default TYPE is double.
• The second function generate a row or column vector of type (with real or complex elements) specified by VEC_TYPE, with n_elem number of elements. By default, VEC_TYPE is vec.
• The third function generate a matrix of type (with real or complex elements) specified by MAT_TYPE, with n_rows number of rows, n_cols number of columns. By default, MAT_TYPE is mat.
• To change the seed of the random number generator, use rdevs.set_seed(YOUR_SEED). To set random seed, use rdevs.set_seed_random(). By default, the random number generator is instantiated with random seed. See RandomDevices class for more details.
• Example:

  		double a = randN(); // random number from a normal distribution 
  		                    // with mean 0 and sd 1
  		
  		complex<double> b = randN<complex<double> >(); // complex random number 
  		                                               // from a normal distribution 
  		                                               // with mean 0 and sd 1
  		
  		double c = randN({-5, 5}); // random number from a normal distribution 
  		                           // with mean -5 and sd 5
  		
  		cx_vec V = randN<cx_vec>(10); // cx_vec with 10 elements 
  
  		cx_rowvec R = randN<cx_rowvec>(10, {-5, 5}); // cx_rowvec with 10 elements
  
  		mat M1 = randN(10, 10); // 10x10 mat
  
  		mat M2 = randN(10, 10, {-5, 5}); // 10x10 mat  
  
  		cx_mat M3 = randN<cx_mat>(10, 10, {-5, 5}); // 10x10 cx_mat  
  	      

• To generate object with uniform random integers in the range specified by the column vector range. By default, the range is in between 0 to 1000.
• The first function generates only one random integer of type (real or complex) specified by ELEM_TYPE. By default TYPE is sword.
• The second function generates a row or column vector of type (with real or complex elements) specified by VEC_TYPE, with n_elem number of elements. By default, VEC_TYPE is ivec.
• The third function generates a matrix of type (with real or complex elements) specified by MAT_TYPE, with n_rows number of rows, n_cols number of columns. By default, MAT_TYPE is imat.
• To change the seed of the random number generator, use rdevs.set_seed(YOUR_SEED). To set random seed, use rdevs.set_seed_random(). By default, the random number generator is instantiated with random seed. See RandomDevices class for more details.
• Example:

  		sword a = randI(); // random integer (arma::sword) between 0 and 10000
  		int a1 = randI<int>(); // random integer (int) between 0 and 10000
  
  		
  		complex<double> b = randI<complex<double> >(); // complex random integer 
  		                                               // (type-casted to complex<double>) 
  		                                               // between (0,0) and (1000,1000)
  		
  		sword c = randI({-100, 100}); // random integer (arma::sword) between -100 and 100
  		
  		ivec V = randI(10); // ivec with 10 elements 
  
  		irowvec R = randI<irowvec>(10, {-100, 100}); // irowvec with 10 elements
  
  		imat M1 = randI(10, 10); // 10x10 imat
  
  		mat M2 = randI<mat>(10, 10, {-100, 100}); // 10x10 mat  
  	      

• To generate random Hermitian (cx_mat) or real symmetric matrix (mat). dim is an unsigned integer (uword) specifying the dimension of the matrix.
• ELEM_TYPE specifies the element type of the matrix (complex (or real) floating point type). By default, it is complex<double>. If ELEM_TYPE is real number (double) the function generates random real symmetric matrix.
• To change the seed of the random number generator, use rdevs.set_seed(YOUR_SEED). To set random seed, use rdevs.set_seed_random(). By default, the random number generator is instantiated with random seed. See RandomDevices class for more details.
• Example:

  		cx_mat H1 = randHermitian(10); // 10x10 Hermitian cx_mat
  
  		mat H2 = randHermitian<double>(10); // 10x10 real symmetric matrix
  	      

• To generate random unitary (cx_mat) or real orthogonal matrix (mat), uniformly distributed according to the Haar measure. dim is an unsigned integer (uword) specifying the dimension of the matrix.
• ELEM_TYPE specifies the element type of the matrix (complex (or real) floating point type). By default, it is complex<double>. If ELEM_TYPE is real number (double) the function generates random real orthogonal matrix.
• To change the seed of the random number generator, use rdevs.set_seed(YOUR_SEED). To set random seed, use rdevs.set_seed_random(). By default, the random number generator is instantiated with random seed. See RandomDevices class for more details.


• Example:

  		cx_mat U1 = randUnitary(10); // 10x10 Unitary cx_mat
  
  		mat U2 = randUnitary<double>(10); // 10x10 real orthogonal matrix
  	      

• To generate random pure states (cx_vec or vec), uniformly distributed according to the Haar measure. dim is an unsigned integer (uword) specifying the dimension of the state.
• ELEM_TYPE specifies the element type of the matrix (complex (or real) floating point type). By default, it is complex<double>.
• To change the seed of the random number generator, use rdevs.set_seed(YOUR_SEED). To set random seed, use rdevs.set_seed_random(). By default, the random number generator is instantiated with random seed. See RandomDevices class for more details.
• Example:

  		cx_vec psi1 = randPsi(10); // pure state (cx_vec) of dimension 10
  
  		vec psi2 = randPsi<double>(10); // real pure state (vec) 
  	      

• To generate random permutation of unsigned (positive) integers of size n. Returns a uvec.
• start (uword) is optional. It specifies starting point of the permutation. By default, it is 0.
• To change the seed of the random number generator, use rdevs.set_seed(YOUR_SEED). To set random seed, use rdevs.set_seed_random(). By default, the random number generator is instantiated with random seed. See RandomDevices class for more details.
• Example:

  		uvec perm1 = randPerm(5); // random permutation of {0,1,2,3,4}
  		
  		uvec perm2 = randPerm(5, 2); // random permutation of {2,3,4,5,6}
  	      

• von Neumann entropy of a quantum state.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• If A is a column vector, return value is always zero.

• Shannon entropy of a probability distribution.
• V has to be a positive real column vector (vec or std::vector of positive real numbers).

• Renyi entropy of a quantum state.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec), and a has to to greater than or equal to zero.
• If A is a column vector, return value is always zero.

• Renyi entropy of a probability distribution.
• V has to be a positive real column vector (vec or std::vector of positive real numbers), and a has to to greater than or equal to zero.

• Tsallis entropy of a quantum state.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec), and a has to to greater than or equal to zero.
• If A is a column vector, return value is always zero.

• Tsallis entropy of a probability distribution.
• V has to be a positive real column vector (vec or std::vector of positive real numbers), and a has to to greater than or equal to zero.

• Quantum mutual information between 2 subsystems of a quantum state.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• In the first case, dim has to be a uvec of two elements, containing the dimensions of two subsystems.
• In the second case, dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem. subsys1 and subsys1 are uvec containing the indices parties that constitute of 1st partition and 2nd partition respectively.
• Indices in subsys1 and subsys2 start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		
  		cx_mat A = randRho(12); // 12x12 density matrix
  		
  		cx_mat B = randRho(16); // 16x16 density matrix
  
  		// mutual_info between two parties with dim 4 and 3
  		double A1 = mutual_info(A, {4, 3}); // 1st party has dim=4,
  		                                    // 2nd party has dim=3
  
  		// mutual_info between 1,3 and 2,4 parties
  		double B1 = mutual_info(B, {1, 3}, {2, 4});
  
  		double B2 = mutual_info(B, {1, 3}, {2, 4}, 2); //same as above
  		                                               // all party has dim = 2
  
  		double B1 = mutual_info(B, {1, 3}, {2, 4}, {2, 2, 2, 2}); // same as above
  		                                                          // explicitly written dimensions
  	      

• Relative entropy of between two quantum state.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec). Same for B.

• Relative entropy of between two probability distribution.
• V1 and V2 have to be positive real column vectors of same type (vec or std::vector of positive real numbers).

• Entanglement entropy of a pure quantum state.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec), and dim has to be a uvec of two elements, containing the dimensions of two subsystems..
• If A is a square matrix and not a density matrix of a pure state, then von Neumann entropy of traced out 2nd party density matrix is returned. To make sure use is_pure function first.

• Negativity of a quantum state, where subsys is a uvec containing the indices of transposed subsystems.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		cx_mat A = randRho(8); // 8x8 density matrix
  		
  		cx_mat B = randRho(12); // 12x12 density matrix
  		
  		// negativity by transposing 2nd and 3rd party
  		double A1 = neg(A, {2, 3}); // all dim are 2, total 3 parties
  
  		// negativity by transposing 2nd and 3rd party
  		double B1 = neg(B, {2, 3}, {2, 3, 2}); // explicitly written dimensions
  	      

• Logarithmic negativity of a quantum state, where subsys is a uvec containing the indices of transposed subsystems.
• A has to be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		cx_mat A = randRho(8); // 8x8 density matrix
  
  		cx_mat B = randRho(12); // 12x12 density matrix
  
  		// logarithmic negativity by transposing 2nd and 3rd party
  		double A1 = log_neg(A, {2, 3}); // all dim are 2, total 3 parties
  
  		// logarithmic negativity by transposing 2nd and 3rd party
  		double B1 = log_neg(B, {2, 3}, {2, 3, 2}); // explicitly written dimensions
  	      

• Concurrence of a two-qubit quantum state.
• A has to be a square matrix (mat or cx_mat) of dimension 4.

• Entanglement of formation of a two-qubit quantum state.
• A has to be a square matrix (mat or cx_mat) or column vector (vec or cx_vec) of dimension 4.

• To check if a bipartite state is entangled or not (based on covariant matrix criteria (Phys. Rev. A 78, 052319 (2008), Proposition IV.2), can be used to detect bound entangled states)
• A has to be a square matrix (mat or cx_mat).
• Returns true if A is entangled, else false.
• The first form assumes both the subsystems have same dimension. dim is an unsigned integer type (uword) stating the dimensions of both the subsystems.
• The second case applies if the subsystems have different dimensions. dim1 and dim2 are unsigned integer types (uword) stating the dimensions of 1st and 2nd party respectively.

• Schmidt decomposition of a bipartite pure state.
• A has to be a square matrix (mat or cx_mat) or column vector (vec or cx_vec).
• dim is a uvec of two elements, containing dimensions of each party.
• S is a real vector containing Schmidt coefficients. U and V are matrices that store Schmidt vectors column-wise, such that A = Σi S(i) * kron(U.col(i), V.col(i)).
• schmidt_full gives full Schmidt basis on both sides (i.e., U and V).
• If A is a square matrix, then it is translated into a column vector using conv_to_pure.
• If decomposition fails then vec S = schmidt(A,dim) throws std::runtime_error exception, and schmidt(A,dim,S,U,V) and schmidt_full(A,dim,S,U,V) returns false and resets S, U, and V.
• Example:

  		cx_vec A = randPsi(8); // 8 dimensional pure state
  
  		vec S1 = schmidt(A, {2, 4});
  
  		vec S;
  		cx_mat U,V;
  		schmidt(A, {2, 4}, S, U, V);
  	      

• Schmidt vectors of a bipartite pure state.
• A has to be a square matrix (mat or cx_mat) or column vector (vec or cx_vec).
• dim is a uvec of two elements, containing dimensions of each party.
• If A is a square matrix, then it is translated into a column vector using conv_to_pure.
• schmidtA(A,dim) returns truncated Schmidt basis on Alice's side in a column-wise manner, while schmidtA_full(A,dim) returns full Schmidt basis on Alice's side.
• schmidtB(A,dim) returns truncated Schmidt basis on Bob's side in a column-wise manner, while schmidtB_full(A,dim) returns full Schmidt basis on Bob's side.
• schmidtAB(A,dim) returns field of matrices containing truncated Schmidt basis on both sides in a column-wise manner, while schmidtAB_full(A,dim) returns field of matrices containing full Schmidt basis on both side.
• If decomposition fails then these functions throw std::runtime_error exception.
• Example:

  		cx_mat A = randPsi(8); // 8 dimensional pure state
  		
  		cx_mat V = schmidtB(A, {2, 4}); // V is 4x2 matrix
  		
  		cx_mat Vf = schmidtB_full(A, {2, 4}); // Vf is 4x4 matrix
  
  		field<cx_mat> AB = schmidtAB(A, {2, 4});
  		//AB(0) == Schmidt basis on Alice's side
  		//AB(1) == Schmidt basis on Bob's side
  	      

• l1-norm coherence of a quantum state.
• A can be either a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• U is a Unitary matrix (mat or cx_mat) to rotate the state for obtaining coherence with respect to the desired basis.

• Relative entropy of coherence of a quantum state.
• A can be either a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• U is a Unitary matrix (mat or cx_mat) to rotate the state for obtaining coherence with respect to the desired basis.

• Hilbert-Schmidt distance between two density matrices.
• A and B have to be square matrices (mat or cx_mat) of same dimension.

• Trace distance between two density matrices.
• A and B have to be square matrices (mat or cx_mat) of same dimension.

• Bures distance between two density matrices.
• A and B have to be square matrices (mat or cx_mat) of same dimension.

• Fidelity between two density matrices.
• A and B have to be square matrices (mat or cx_mat) of same dimension.

• To apply a quantum gate U or a set of Kraus operators Ksto a quantum state A.
• A has be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• U has to a square matrix (mat or cx_mat), and Ks can be a field or a std::vector of square matrices (mat or cx_mat) having same dimensions.
• In the first three cases, subsys is a uvec containing the indices of subsystems, where the gate or the Kraus operators will be applied.
• Indices in subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• In the first three cases, the dimension of the gate U or the Kraus operators Ks must match the dimension of subsystems specified in subsys.
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• In the last case, the set of Kraus operators Ks acts on the whole quantum state A. The dimension of the Kraus operators Ks must match the dimension of the quantum state A.
• Note: Return type is deduced from the element types of A and U/Ks. For simplicity use auto to deduce the return type.
• Example:

  		cx_mat A = randRho(8); // 8x8 density matrix
  		cx_mat B = randRho(2); // 2x2 density matrix
  
  		const mat& U1 = spm.S(1); // \sigma_x gate
  		cx_mat U2 = kron(spm.S(1),spm.S(1)); // kron(\sigma_x,\sigma_x) gate 
  
  		vector<cx_mat> BF = {(1 - 0.25) * spm.S(0), 0.25 * spm.S(1)};
  		// Kraus operators for Bit-Flip channel with p = 0.25 
  					
  		cx_mat A1 = apply(A, U1, {2}); // all dim are 2, total 3 parties
  		cx_mat A2 = apply(A, U2, {1, 3}, {2, 2, 2}); // explicitly written dimensions
  		
  		cx_mat A3 = apply(A, BF, {2}, {2, 2, 2}); // apply Bit-Flip on 2nd party
  		
  		cx_mat B1 = apply(B, BF); // apply Bit-Flip on B
  	      

• To apply a quantum gate U to a quantum state A as a controlled gate. subsys is a uvec containing the indices subsystems, where the gate will be applied, and ctrl is a uvec containing the indices of control subsystems.
• A has be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• U to has be a square matrix (mat or cx_mat) and the dimension of gate U must match the dimension of subsystems specified in subsys.
• All control subsystems must have same dimensions.
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in subsys and ctrl start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Note: Return type is deduced from the element types of A and U. For simplicity use auto to deduce the return type.
• Example:

  		cx_mat A = randRho(8); // 8x8 density matrix
  
  		const mat& U1 = gates.X; // \sigma_x gate
  		
  		mat U2 = kron(gates.X, gates.X); // kron(\sigma_x,\sigma_x) gate 
  
  		cx_mat A1 = apply_ctrl(A, U1, {1}, {2}); // all dim are 2, total 3 parties
  
  		cx_mat A2 = apply_ctrl(A, U2, {2}, {1, 3}, {2, 2, 2}); // explicitly written dimensions
  	      

• To convert a quantum gate U into a controlled quantum gate. subsys is a uvec containing the indices subsystems, where the gate will be applied, and ctrl is a uvec containing the indices of control subsystems.
• U to has be a square matrix (mat or cx_mat) and the dimension of gate U must match the dimension of subsystems specified in subsys.
• All control subsystems must have same dimensions.
• In the first case, n is a uword mentioning total number of subsystems, having same dimensions equal to dim (uword). By default, dim is 2.
• In the second case, dim is a uvec containing dimensions of every subsystem. These form is needed if different subsystems have different dimensions.
• Indices in subsys and ctrl start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		const mat& U = gates.X; // \sigma_x gate
  		
  		mat CU1 = make_ctrl(U, {1}, {2}, 2); // make two-party controlled X gate
  		                                     // where 1st party is control and
  		                                     // second party is target.
  		                                     // Here, n = 2 and 
  		                                     // dim = 2 (by default)
  
  		mat CU2 = make_ctrl(U, {1}, {2}, {2, 2}); //same as above 
  		                                          // explicitly written dimensions
  	      

• To measure a quantum state A using a set of Kraus or projection operators Ks. Returns std::tuple of: 1. result of the measurement (uword), 2. vector of outcome probabilities (vec), and 3. field of post-measurement normalized states (field<cx_mat> or field<mat>).
• A has be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• Ks can be a std::vector or a field of Kraus or projection operators (either cx_mat, mat, cx_vec or vec). It can also be a cx_mat or mat, where each column will be treated as projection operators.
• First two functions measure the whole quantum state A, and the dimension of the operators in Ks must match the dimension of A.
• Last two functions partially measure the subsystems of the state A. The indices of measured subsystems are specified in subsys (uvec). The dimension of the operators in Ks must match the dimension of subsystems specified in subsys.
• All Kraus or projection operators must have same dimensions.
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Note: MATRIX_TYPE is deduced from the element types of A and Ks. For simplicity use auto to deduce return type.
• Example:

  		cx_mat A1 = randRho(2); // 2x2 random density matrix
  
  		std::vector Ks = {spm.basis2(0,1), spm.basis2(1,1)};
  		  // std::vector of projection operators in \sigma_x direction
  
  		auto result1 = measure(A1, Ks); // measurement
  		std::get<0>(A1); // measurement result
  		std::get<1>(A1); // outcome probabilities
  		std::get<2>(A1); // post measurement normalised states
  			
  		cx_mat A2 = randRho(8); //8x8 random density matrix
  
  		auto result2 = measure(A2, Ks, {2}, {2,2,2}); measure the 2nd party
  	      

• To measure a quantum state A in the computational basis. Returns std::tuple of: 1. result of the measurement (uword) and 2. vector of outcome probabilities (vec).
• A has be a square matrix (mat or cx_mat) or a column vector (vec or cx_vec).
• The first function measures the whole quantum state, where as the last function partially measures the subsystems of the state A, where the indices of measured subsystems are specified in subsys.
• dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Indices in subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
• Example:

  		cx_mat A1 = randRho(2); // 2x2 random density matrix
  
  		auto result1 = measure_comp(A1); // measurement
  		std::get<0>(A1); // measurement result
  		std::get<1>(A1); //  outcome probabilities
  		    
  		cx_mat A2 = randRho(8); //8x8 random density matrix
  
  		auto result2 = measure_comp(A2, {2}, {2,2,2}); measure the 2nd party
  	      

Classes for Quantum Discord like features

• Note: These features depend on NLopt.
• Note: If you want to use old API's for discord like function, define QICLIB_USE_OLD_DISCORD macro before including QIClib header. See preprocessor macros in QIClib for details. Old API's for discord like functions are not maintained anymore.

• discord_space provides computational space for calculating quantum discord (actual as well as constrained-regular, as discussed in Phys. Rev. A 92, 062301 (2015), Sec. III.A, Case 4) for quantum states when measurement is done over a qubit or qutrit system, deficit_space does the same for one-way quantum work deficit. Both have same syntaxes, so for the demonstration purpose, we will talk about discord_space only.
• Works for real or complex density matrices (mat or cx_mat). MATRIX_TYPE denotes the the data type of the density matrix i.e., mat or cx_mat.
• Constructors:
  discord_space<MATRIX_TYPE>(A, subsys, dim)
discord_space<MATRIX_TYPE>(A, subsys)
  
  • A is a square matrix of type MATRIX_TYPE.
  • subsys is a unsigned integer type (uword), which denotes the index of measured subsystem. Indices for subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
  • dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Member functions:
  
  • .result(): Computes quantum discord, and returns the result. This function will keep the value cached in the memory. Further use of this function will return the cached value, rather than computing discord all over again.
  • .result_reg(): Computes constrained-regular quantum discord, and returns the result. This function will keep the value cached in the memory. Further use of this function will return the cached value, rather than computing constrained-regular quantum discord all over again.
  • .compute(): Computes quantum discord whenever called. This function will rewrite the cached result to the new one. One can use .result() to get the new value (see the following example).
  • .compute_reg(): Computes constrained-regular quantum discord whenever called. This function will rewrite the cached result to the new one. One can use .result_reg() to get the new value (see the following example).
  • .reset(): This will reset the state of discord_space to the default settings. All cached values will be forgotten.
  • .reset(subsys): This will reset the state of discord_space to the default settings. All cached values will be forgotten. subsys is a unsigned integer type (uword), which denotes the index of new measured subsystem. Indices for subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on.
  • .reset(A, subsys, dim): This will reset the state of discord_space to the default settings. All cached values will be forgotten. A is the new density matrix (square matrix) of type MATRIX_TYPE. subsys is a unsigned integer type (uword), which denotes the index of new measured subsystem. Indices for subsys start from 1, i.e., first party has the index 1, second party has the index 2 and so on. dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.dim is optional. It can be either unsigned integer type (uword) or uvec containing dimensions of every subsystem. When it is unsigned integer type, all subsystems are supposed to have same dimensions (equal to dim). By default, it is 2. If different subsystems have different dimensions, use dim as a uvec containing dimensions of every subsystem.
• Helper member functions: These functions are there to change the settings of the optimization processes done by NLopt, to fine tune the result. In general, users can overlook these functions, and default settings will be used.
  
  • .use_global_opt(bool a = true): Set whether global optimization will be performed (true) or not (false).
  • .global_algorithm(nlopt::algorithm a = nlopt::GN_DIRECT_L_RAND): Set global optimization algorithm.
  • .global_xtol(double a): Set global optimization xtol. By default, it is 4e-2 for qubit systems, and 0.25 for qutrit systems.
  • .global_ftol(double a = 0): Set global optimization ftol.
  • .local_algorithm(nlopt::algorithm a = nlopt::LN_COBYLA): Set local optimization algorithm.
  • .local_xtol(double a): Set local optimization xtol.
  • .local_ftol(double a = 0): Set local optimization ftol.
  • .angle_range(const vec& a): Set the ranges of angles (parameters of unitary matrices) as multiples of π. For qubit systems, default is {1, 2}, and for qutrit systems, default is {2, 2, 2, 2, 2}.
  • .initial_angle(const vec& a): Set the initial values of angles as multiples of π. For qubit systems, default is {0.1, 0.1}, and for qutrit systems, default is {0.1, 0.1, 0.1, 0.1, 0.1}.
• Example:

  		cx_mat A = randRho(8); // 8x8 density matrix (3-qubit)
  
  		// create discord_space, we measure the 2nd party to compute discord
  		discord_space DSpace(A, {2});		
  		
  		// calculate discord
  		double DResult = DSpace.result();
  
  		// change local optimization algorithm,
  		// then compute the discord,
  		// and get the result in one line --->
  		double DResult = DSpace.local_algorithm(nlopt::LN_SBPLX).compute().result().
  	
  		
  		cx_mat B = randRho(12); // new 2x3x2 density matrix
  		
  		// using same discord_space for new density matrix,
  		// measurement on 2nd party (qutrit)
  		DSpace.reset(B, {2}, {2, 3, 2});
  		
  		double new_DResult = Dspace.result(); 
  
  	      

Frequently asked questions

• What are the basic dependencies of QIClib?
  
  • Armadillo C++ linear algebra library (version 4.2 or later), with LAPACK and BLAS implementations.
  • A C++11 compliant compiler. gcc version 4.8 or later, or clang version 3.3 or later are recommended.
  • NLopt nonlinear optimization library for certain features. See getting started section.
• How do I use QIClib in my project?
  See getting started section.
• How do I report bugs?
  See this.
• How do I use XYZ function properly?
  See API documentation. If you still have concerns, see this.
• Can you implement certain features, which are useful for my work?
  Yes, see this.
• Instead of using other popular C++ linear algebra libraries, why did QIClib choose Armadillo?
  • Armadillo is fast, efficient and reliable.
  • It has nice API, which is deliberately similar to Matlab.
  • It supports modern C++11 features, like move semantics or initialiser list initialization, which indeed increase efficiency and usability.
  • It has rapid developement process.
  • It can make use of high performance multi-threaded LAPACK and BLAS replacements like OpenBLAS, Intel MKL, or AMD ACML.
• I know C++98, but do not know much about C++11. Can I use QIClib?
  Yes. Though QIClib is written in C++11 standard, you can still do C++98 style coding with QIClib functions. You just need a C++11 compliant compiler and make sure that you add necessary compiler flag to enable C++11 features, e.g., in gcc or clang add -std=c++11 flag during compilation.
• Can I use C++11 auto with QIClib functions?
  Though use of C++11 auto is not recommended with Armadillo expressions, due to the extensive use of template meta-programming, every QIClib functions return by value, so use of auto is safe there. In fact, use of auto is recommended with QIClib functions and objects, as it will always guarantee move operations. But again, care should be taken when auto is used with Armadillo expressions. See the following example:
  

  		cx_mat A = randn<cx_mat>(8,8);
  		A *= A.t();
  		A /= trace(A);
  
  		auto B = A; // OK, no temporary expression
  		auto C = TrX(A,{1},{2,2,2}); // OK, QIClib function
  		auto D = powm_sym(A,3); // OK, QIClib function
  
  		auto E = A * A; // NOT OK!! temporary expression
  		auto F = A + B + D ; // NOT OK!! temporary expression
  		auto G = A * A.t(); // NOT OK!! temporary expression
  
  		vec V = {0,1,2,3,4};
  
  		for(auto&& i : V) // OK, no temporary exression
  		{
  		// do something
  		}
  		
  		for(auto&& i : V+V) // NOT OK!! temporary expression
  		{
  		// do something
  		}
  
  		auto& S1 = spm.S(1); // OK and recommended
  		const cx_mat& S2 = spm.S(2); // also OK  
  		
  	      

• Does QIClib use C++14 features?
  No, only C++11 features have been used. C++14 features like auto return type deductions for functions have been avoided by using template meta-programming. But users are requested to use latest C++14/17 features, which increase efficiency and usability.
• Does QIClib use multi-threaded algorithms for computations?
  Yes and No. There are several ways to develop multi-threaded or parallel applications using QIClib.
  • Using high performance multi-threaded LAPACK and BLAS replacements like OpenBLAS, Intel MKL, or AMD ACML.
  • Loop parallelization using OpenMP or others similar libraries.
  • TrX, Tx, sysperm, apply, apply_ctrl, make_ctrl, measure, and measure_comp have optional multi-threaded algorithms if OpenMP is enabled in compiler argument, i.e., -fopenmp in gcc or clang. For this, define QICLIB_USE_OPENMP macro before before including QIClib header. One can individually enable multi-threading in those mentioned functions. For this, see preprocessor macros in QIClib for details.
  Note: These three procedures of parallelization do not mix well. Use only one of the procedures. For best performance, try to use only one of first two methods.

Helper preprocessor macros

API Additions

• Bug fix in tsallis_prob.
• Added gram_schmidt.
• Added dense_to_sparse, sparse_to_dense.
• Added discord_space, deficit_space.
• Added is_Normal.
• Added rel_entropy, rel_entropy_prob.
• Added coherence measures, l1_coh and rel_entropy_coh.
• Added schmidt_full.
• Added macro features to fine tune the library behavior.
• std::runtime_error exception for eig_sym/eig_gen decomposition.
• std::runtime_error exception for qr decomposition.
• Added randPerm.
• Added make_ctrl.
• Default template parameter of mket, mproj is now double.
• Optional multi-threaded algorithm for TrX, Tx, sysperm, apply, apply_ctrl, make_ctrl, measure, and measure_comp.

• Added purify, std_to_HS, and HS_to_std.
• Added measure and measure_comp.
• Added random number/matrix generators.
• Added is_diagonalizable.
• Faster is_equal, is_Hermitian, is_Unitary, is_pure, is_valid_state.
• Fixed a bug in discord like functions.

• dsum and dsum_pow added.
• Fixed a bug in tensor.

• absm added.
• schatten added.