Conjugate transpose: Difference between revisions

m
(Added 11l)
m (→‎{{header|Wren}}: Minor tidy)
 
(10 intermediate revisions by 7 users not shown)
Line 38:
{{trans|Nim}}
 
<langsyntaxhighlight lang="11l">-V eps = 1e-10
 
F to_str(m)
Line 131:
test(M3)
print("\n")
test(M4)</langsyntaxhighlight>
 
{{out}}
Line 185:
 
=={{header|Ada}}==
<langsyntaxhighlight Adalang="ada">with Ada.Text_IO; use Ada.Text_IO;
with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;
Line 231:
Put_Line("nmat:"); Examine(nmat); New_Line;
Put_Line("umat:"); Examine(umat);
end ConTrans;</langsyntaxhighlight>
{{out}}
<pre>hmat:
Line 269:
=={{header|ALGOL 68}}==
Uses the same test cases as the Ada sample.
<langsyntaxhighlight lang="algol68">BEGIN # find and classify the complex conjugate transpose of a complex matrix #
# returns the conjugate transpose of m #
OP CONJUGATETRANSPOSE = ( [,]COMPL m )[,]COMPL:
Line 367:
)
)
END</langsyntaxhighlight>
{{out}}
<pre>
Line 399:
 
=={{header|C}}==
<langsyntaxhighlight lang="c">/* Uses C99 specified complex.h, complex datatype has to be defined and operation provided if used on non-C99 compilers */
 
#include<stdlib.h>
Line 579:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>
Line 600:
 
=={{header|C++}}==
<langsyntaxhighlight lang="cpp">#include <cassert>
#include <cmath>
#include <complex>
Line 790:
test(matrix3);
return 0;
}</langsyntaxhighlight>
 
{{out}}
Line 832:
 
=={{header|Common Lisp}}==
<syntaxhighlight lang="lisp">
<lang Lisp>
(defun matrix-multiply (m1 m2)
(mapcar
Line 863:
(defun unitary-p (m)
(identity-p (matrix-multiply m (conjugate-transpose m))) )
</syntaxhighlight>
</lang>
 
{{out}}
Line 886:
=={{header|D}}==
{{trans|Python}} A well typed and mostly imperative version:
<langsyntaxhighlight lang="d">import std.stdio, std.complex, std.math, std.range, std.algorithm,
std.numeric;
 
Line 992:
writefln("Unitary? %s.\n", isUnitary(mat, ct));
}
}</langsyntaxhighlight>
{{out}}
<pre>Matrix:
Line 1,031:
===Alternative Version===
A more functional version that contains some typing problems (same output).
<langsyntaxhighlight lang="d">import std.stdio, std.complex, std.math, std.range, std.algorithm,
std.numeric, std.exception, std.traits;
 
Line 1,120:
writefln("Unitary? %s.\n", isUnitary(mat, ct));
}
}</langsyntaxhighlight>
 
=={{header|F_Sharp|F#}}==
<syntaxhighlight lang="fsharp">
// Conjugate transpose. Nigel Galloway: January 10th., 2022
let fN g=let g=g|>List.map(List.map(fun(n,g)->System.Numerics.Complex(n,g)))|>MathNet.Numerics.LinearAlgebra.MatrixExtensions.matrix in (g,g.ConjugateTranspose())
let fG n g=(MathNet.Numerics.LinearAlgebra.Matrix.inverse n-g)|>MathNet.Numerics.LinearAlgebra.Matrix.forall(fun(n:System.Numerics.Complex)->abs n.Real<1e-14&&abs n.Imaginary<1e-14)
let test=[fN [[(3.0,0.0);(2.0,1.0)];[(2.0,-1.0);(1.0,0.0)]];fN [[(1.0,0.0);(1.0,0.0);(0.0,0.0)];[(0.0,0.0);(1.0,0.0);(1.0,0.0)];[(1.0,0.0);(0.0,0.0);(1.0,0.0)]];fN [[(1.0/sqrt 2.0,0.0);(1.0/sqrt 2.0,0.0);(0.0,0.0)];[(0.0,1.0/sqrt 2.0);(0.0,-1.0/sqrt 2.0);(0.0,0.0)];[(0.0,0.0);(0.0,0.0);(0.0,1.0)]]]
test|>List.iter(fun(n,g)->printfn $"Matrix\n------\n%A{n}\nConjugate transposed\n--------------------\n%A{g}\nIs hermitian: %A{n.IsHermitian()}\nIs normal: %A{n*g=g*n}\nIs unitary: %A{fG n g}\n")
</syntaxhighlight>
{{out}}
<pre>
Matrix
------
DenseMatrix 2x2-Complex
(3, 0) (2, 1)
(2, -1) (1, 0)
 
Conjugate transposed
--------------------
DenseMatrix 2x2-Complex
(3, -0) (2, 1)
(2, -1) (1, -0)
 
Is hermitian: true
Is normal: true
Is unitary: false
 
Matrix
------
DenseMatrix 3x3-Complex
(1, 0) (1, 0) (0, 0)
(0, 0) (1, 0) (1, 0)
(1, 0) (0, 0) (1, 0)
 
Conjugate transposed
--------------------
DenseMatrix 3x3-Complex
(1, -0) (0, -0) (1, -0)
(1, -0) (1, -0) (0, -0)
(0, -0) (1, -0) (1, -0)
 
Is hermitian: false
Is normal: true
Is unitary: false
 
Matrix
------
DenseMatrix 3x3-Complex
(0.707107, 0) (0.707107, 0) (0, 0)
(0, 0.707107) (0, -0.707107) (0, 0)
(0, 0) (0, 0) (0, 1)
 
Conjugate transposed
--------------------
DenseMatrix 3x3-Complex
(0.707107, -0) (0, -0.707107) (0, -0)
(0.707107, -0) (0, 0.707107) (0, -0)
(0, -0) (0, -0) (0, -1)
 
Is hermitian: false
Is normal: true
Is unitary: true
</pre>
=={{header|Factor}}==
Before the fix to [https://github.com/slavapestov/factor/issues/484 Factor bug #484], <code>m.</code> gave the wrong answer and this code failed. Factor 0.94 is too old to work.
 
{{works with|Factor|development (future 0.95)}}
<langsyntaxhighlight lang="factor">USING: kernel math.functions math.matrices sequences ;
IN: rosetta.hermitian
 
Line 1,139 ⟶ 1,201:
 
: unitary-matrix? ( matrix -- ? )
[ dup conj-t m. ] [ length identity-matrix ] bi = ;</langsyntaxhighlight>
 
Usage:
Line 1,157 ⟶ 1,219:
The examples and algorithms are taken from the j solution, except for UnitaryQ. The j solution uses the matrix inverse verb. Compilation on linux, assuming the program is file f.f08 :<pre>
gfortran -std=f2008 -Wall -fopenmp -ffree-form -fall-intrinsics -fimplicit-none f.f08 -o f</pre>
<syntaxhighlight lang="fortran">
<lang FORTRAN>
program conjugate_transpose
 
Line 1,222 ⟶ 1,284:
 
end program conjugate_transpose
</syntaxhighlight>
</lang>
<pre>
-*- mode: compilation; default-directory: "/tmp/" -*-
Line 1,269 ⟶ 1,331:
 
=={{header|FreeBASIC}}==
<langsyntaxhighlight lang="freebasic">'complex type and operators for it
type complex
real as double
Line 1,398 ⟶ 1,460:
print is_hermitian(A), is_normal(A), is_unitary(A)
print is_hermitian(B), is_normal(B), is_unitary(B)
print is_hermitian(C), is_normal(C), is_unitary(C)</langsyntaxhighlight>
{{out}}
<pre>true true true
Line 1,405 ⟶ 1,467:
 
=={{header|Go}}==
<langsyntaxhighlight lang="go">package main
 
import (
Line 1,516 ⟶ 1,578:
}
return m3
}</langsyntaxhighlight>
Output:
<pre>
Line 1,559 ⟶ 1,621:
=={{header|Haskell}}==
Slow implementation using lists.
<langsyntaxhighlight lang="haskell">import Data.Complex (Complex(..), conjugate)
import Data.List (transpose)
 
Line 1,633 ⟶ 1,695:
:: Num a
=> Matrix (Complex a) -> Matrix (Complex a)
conjTranspose = map (map conjugate) . transpose</langsyntaxhighlight>
Output:
<pre>
Line 1,673 ⟶ 1,735:
=={{header|J}}==
 
'''SolutionConjugate transpose:''': <langsyntaxhighlight lang="j"> ct =: +@|: NB. Conjugate transpose (ct A is A_ct)</langsyntaxhighlight>
 
'''Examples''': <lang j> X =: +/ . * NB. Matrix Multiply (x)
<code>+</code> when used without a left argument is conjugate, <code>|:</code> is transpose, and <code>@</code> composes functions.
 
'''Examples''': <syntaxhighlight lang="j"> X =: +/ . * NB. Matrix Multiply (x)
 
HERMITIAN =: 3 2j1 ,: 2j_1 1
Line 1,686 ⟶ 1,751:
UNITARY =: (-:%:2) * 1 1 0 , 0j_1 0j1 0 ,: 0 0 0j1 * %:2
(ct -: %.) UNITARY NB. A_ct = A^-1
1</langsyntaxhighlight>
 
'''Reference''' (with example matrices for other langs to use):<syntaxhighlight lang ="j"> HERMITIAN;NORMAL;UNITARY
+--------+-----+--------------------------+
| 3 2j1|1 1 0| 0.707107 0.707107 0|
Line 1,703 ⟶ 1,768:
+-----+-----+-----+
|1 1 0|0 1 0|0 1 1|
+-----+-----+-----+</langsyntaxhighlight>
 
=={{header|Java}}==
<syntaxhighlight lang ="java">
 
import java.util.Arrays;
import java.util.List;
 
public final class ConjugateTranspose {
 
public static void main(String[] aArgs) {
ComplexMatrix one = new ComplexMatrix( new Complex[][] { { new Complex(0, 4), new Complex(-1, 1) },
{ new Complex(1, -1), new Complex(0, 4) } } );
 
ComplexMatrix two = new ComplexMatrix(
new Complex[][] { { new Complex(1, 0), new Complex(1, 1), new Complex(0, 2) },
{ new Complex(1, -1), new Complex(5, 0), new Complex(-3, 0) },
{ new Complex(0, -2), new Complex(-3, 0), new Complex(0, 0) } } );
 
final double term = 1.0 / Math.sqrt(2.0);
ComplexMatrix three = new ComplexMatrix( new Complex[][] { { new Complex(term, 0), new Complex(term, 0) },
{ new Complex(0, term), new Complex(0, -term) } } );
List<ComplexMatrix> matricies = List.of( one, two, three );
for ( ComplexMatrix matrix : matricies ) {
System.out.println("Matrix:");
matrix.display();
System.out.println("Conjugate transpose:");
matrix.conjugateTranspose().display();
System.out.println("Hermitian: " + matrix.isHermitian());
System.out.println("Normal: " + matrix.isNormal());
System.out.println("Unitary: " + matrix.isUnitary() + System.lineSeparator());
}
}
 
}
 
final class ComplexMatrix {
public ComplexMatrix(Complex[][] aData) {
rowCount = aData.length;
colCount = aData[0].length;
data = Arrays.stream(aData).map( row -> Arrays.copyOf(row, row.length) ).toArray(Complex[][]::new);
}
public ComplexMatrix multiply(ComplexMatrix aOther) {
if ( colCount != aOther.rowCount ) {
throw new RuntimeException("Incompatible matrix dimensions.");
}
Complex[][] newData = new Complex[rowCount][aOther.colCount];
Arrays.stream(newData).forEach( row -> Arrays.fill(row, new Complex(0, 0)) );
for ( int row = 0; row < rowCount; row++ ) {
for ( int col = 0; col < aOther.colCount; col++ ) {
for ( int k = 0; k < colCount; k++ ) {
newData[row][col] = newData[row][col].add(data[row][k].multiply(aOther.data[k][col]));
}
}
}
return new ComplexMatrix(newData);
}
 
public ComplexMatrix conjugateTranspose() {
if ( rowCount != colCount ) {
throw new IllegalArgumentException("Only applicable to a square matrix");
}
Complex[][] newData = new Complex[colCount][rowCount];
for ( int row = 0; row < rowCount; row++ ) {
for ( int col = 0; col < colCount; col++ ) {
newData[col][row] = data[row][col].conjugate();
}
}
return new ComplexMatrix(newData);
}
public static ComplexMatrix identity(int aSize) {
Complex[][] data = new Complex[aSize][aSize];
for ( int row = 0; row < aSize; row++ ) {
for ( int col = 0; col < aSize; col++ ) {
data[row][col] = ( row == col ) ? new Complex(1, 0) : new Complex(0, 0);
}
}
return new ComplexMatrix(data);
}
public boolean equals(ComplexMatrix aOther) {
if ( aOther.rowCount != rowCount || aOther.colCount != colCount ) {
return false;
}
for ( int row = 0; row < rowCount; row++ ) {
for ( int col = 0; col < colCount; col++ ) {
if ( data[row][col].subtract(aOther.data[row][col]).modulus() > EPSILON ) {
return false;
}
}
}
return true;
}
public void display() {
for ( int row = 0; row < rowCount; row++ ) {
System.out.print("[");
for ( int col = 0; col < colCount - 1; col++ ) {
System.out.print(data[row][col] + ", ");
}
System.out.println(data[row][colCount - 1] + " ]");
}
}
public boolean isHermitian() {
return equals(conjugateTranspose());
}
public boolean isNormal() {
ComplexMatrix conjugateTranspose = conjugateTranspose();
return multiply(conjugateTranspose).equals(conjugateTranspose.multiply(this));
}
public boolean isUnitary() {
ComplexMatrix conjugateTranspose = conjugateTranspose();
return multiply(conjugateTranspose).equals(identity(rowCount)) &&
conjugateTranspose.multiply(this).equals(identity(rowCount));
}
private final int rowCount;
private final int colCount;
private final Complex[][] data;
private static final double EPSILON = 0.000_000_000_001;
}
 
final class Complex {
public Complex(double aReal, double aImag) {
real = aReal;
imag = aImag;
}
public Complex add(Complex aOther) {
return new Complex(real + aOther.real, imag + aOther.imag);
}
public Complex multiply(Complex aOther) {
return new Complex(real * aOther.real - imag * aOther.imag, real * aOther.imag + imag * aOther.real);
}
public Complex negate() {
return new Complex(-real, -imag);
}
public Complex subtract(Complex aOther) {
return this.add(aOther.negate());
}
public Complex conjugate() {
return new Complex(real, -imag);
}
public double modulus() {
return Math.hypot(real, imag);
}
public boolean equals(Complex aOther) {
return real == aOther.real && imag == aOther.imag;
}
@Override
public String toString() {
String prefix = ( real < 0.0 ) ? "" : " ";
String realPart = prefix + String.format("%.3f", real);
String sign = ( imag < 0.0 ) ? " - " : " + ";
return realPart + sign + String.format("%.3f", Math.abs(imag)) + "i";
}
private final double real;
private final double imag;
}
</syntaxhighlight>
{{ out }}
<pre>
Matrix:
[ 0.000 + 4.000i, -1.000 + 1.000i ]
[ 1.000 - 1.000i, 0.000 + 4.000i ]
Conjugate transpose:
[ 0.000 - 4.000i, 1.000 + 1.000i ]
[-1.000 - 1.000i, 0.000 - 4.000i ]
Hermitian: false
Normal: true
Unitary: false
 
Matrix:
[ 1.000 + 0.000i, 1.000 + 1.000i, 0.000 + 2.000i ]
[ 1.000 - 1.000i, 5.000 + 0.000i, -3.000 + 0.000i ]
[ 0.000 - 2.000i, -3.000 + 0.000i, 0.000 + 0.000i ]
Conjugate transpose:
[ 1.000 + 0.000i, 1.000 + 1.000i, 0.000 + 2.000i ]
[ 1.000 - 1.000i, 5.000 + 0.000i, -3.000 + 0.000i ]
[ 0.000 - 2.000i, -3.000 + 0.000i, 0.000 + 0.000i ]
Hermitian: true
Normal: true
Unitary: false
 
Matrix:
[ 0.707 + 0.000i, 0.707 + 0.000i ]
[ 0.000 + 0.707i, 0.000 - 0.707i ]
Conjugate transpose:
[ 0.707 + 0.000i, 0.000 - 0.707i ]
[ 0.707 + 0.000i, 0.000 + 0.707i ]
Hermitian: false
Normal: true
Unitary: true
</pre>
 
=={{header|jq}}==
Line 1,715 ⟶ 1,992:
 
If your jq does not have "transpose" then the following may be used:
<langsyntaxhighlight lang="jq"># transpose/0 expects its input to be a rectangular matrix
# (an array of equal-length arrays):
def transpose:
if (.[0] | length) == 0 then []
else [map(.[0])] + (map(.[1:]) | transpose)
end ;</langsyntaxhighlight>
'''(2) Operations on real/complex numbers'''
<langsyntaxhighlight lang="jq"># x must be real or complex, and ditto for y;
# always return complex
def plus(x; y):
Line 1,748 ⟶ 2,025:
if type == "number" then [.,0]
else [.[0], -(.[1]) ]
end;</langsyntaxhighlight>
'''(3) Array operations'''
<langsyntaxhighlight lang="jq"># a and b are arrays of real/complex numbers
def dot_product(a; b):
a as $a | b as $b
| reduce range(0;$a|length) as $i
(0; . as $s | plus($s; multiply($a[$i]; $b[$i]) ));</langsyntaxhighlight>
'''(4) Matrix operations'''
<langsyntaxhighlight lang="jq"># convert a matrix of mixed real/complex entries to all complex entries
def to_complex:
def toc: if type == "number" then [.,0] else . end;
Line 1,793 ⟶ 2,070:
reduce range(0;M|length) as $i
(0; reduce range(0; M[0]|length) as $j
(.; 0 + sqdiff( M[$i][$j]; N[$i][$j] ) ) ) <= epsilon;</langsyntaxhighlight>
====Conjugate transposition====
<langsyntaxhighlight lang="jq"># (entries may be real and/or complex)
def conjugate_transpose:
map( map(conjugate) ) | transpose;
Line 1,817 ⟶ 2,094:
| complex_identity(length) as $I
| approximately_equal( $I; matrix_multiply($H;$M); 1e-10)
and approximately_equal( $I ; matrix_multiply($M;$H); 1e-10) ; </langsyntaxhighlight>
 
====Examples====
<langsyntaxhighlight lang="jq">def hermitian_example:
[ [ 3, [2,1]],
[[2,-1], 1 ] ];
Line 1,848 ⟶ 2,125:
;
 
demo</langsyntaxhighlight>
{{out}}
<langsyntaxhighlight lang="sh">$ jq -r -c -n -f Conjugate_transpose.jq
Hermitian example:
 
Line 1,865 ⟶ 2,142:
Normal example: true
 
Unitary example: true</langsyntaxhighlight>
 
=={{header|Julia}}==
Julia has a built-in matrix type, and the conjugate-transpose of a complex matrix <code>A</code> is simply:
<syntaxhighlight lang ="julia">A'</langsyntaxhighlight>
(similar to Matlab). You can check whether <code>A</code> is Hermitian via the built-in function
<syntaxhighlight lang ="julia">ishermitian(A)</langsyntaxhighlight>
Ignoring the possibility of roundoff errors for floating-point matrices (like most of the examples in the other languages), you can check whether a matrix is normal or unitary by the following functions
<langsyntaxhighlight lang="julia">eye(A) = A^0
isnormal(A) = size(A,1) == size(A,2) && A'*A == A*A'
isunitary(A) = size(A,1) == size(A,2) && A'*A == eye(A)</langsyntaxhighlight>
 
=={{header|Kotlin}}==
As Kotlin doesn't have built in classes for complex numbers or matrices, some basic functionality needs to be coded in order to tackle this task:
<langsyntaxhighlight lang="scala">// version 1.1.3
 
typealias C = Complex
Line 2,003 ⟶ 2,280:
println("Unitary? ${mct.isUnitary()}\n")
}
}</langsyntaxhighlight>
 
{{out}}
Line 2,050 ⟶ 2,327:
=={{header|Maple}}==
The commands <code>HermitianTranspose</code> and <code>IsUnitary</code> are provided by the <code>LinearAlgebra</code> package.
<langsyntaxhighlight Maplelang="maple">M:=<<3|2+I>,<2-I|1>>:
 
with(LinearAlgebra):
Line 2,059 ⟶ 2,336:
type(M,'Matrix'(hermitian));
IsNormal(M);
IsUnitary(M);</langsyntaxhighlight>
Output:
<pre> [ 3 2 + I]
Line 2,076 ⟶ 2,353:
 
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<langsyntaxhighlight Mathematicalang="mathematica">NormalMatrixQ[a_List?MatrixQ] := Module[{b = Conjugate@Transpose@a},a.b === b.a]
UnitaryQ[m_List?MatrixQ] := (Conjugate@Transpose@m.m == IdentityMatrix@Length@m)
 
Line 2,092 ⟶ 2,369:
 
{HermitianMatrixQ@#, NormalMatrixQ@#, UnitaryQ@#}&@m
-> {False, False, False}</langsyntaxhighlight>
 
=={{header|Nim}}==
Line 2,098 ⟶ 2,375:
The complex type is defined as generic regarding the type of real an imaginary part. We have chosen to use Complex[float] and make only our Matrix type generic regarding the dimensions. Thus, a Matrix has a two dimensions M and N which are static, i.e. known at compile time. We have enforced the condition M = N for square matrices (also at compile time).
 
<langsyntaxhighlight Nimlang="nim">import complex, strformat
 
type Matrix[M, N: static Positive] = array[M, array[N, Complex[float]]]
Line 2,250 ⟶ 2,527:
test(M2)
test(M3)
test(M4)</langsyntaxhighlight>
 
{{out}}
Line 2,314 ⟶ 2,591:
 
=={{header|PARI/GP}}==
<syntaxhighlight lang="text">conjtranspose(M)=conj(M~)
isHermitian(M)=M==conj(M~)
isnormal(M)=my(H=conj(M~));H*M==M*H
isunitary(M)=M*conj(M~)==1</langsyntaxhighlight>
 
=={{header|Perl}}==
In general, using two or more modules which overload operators can be problematic. For this task, using both Math::Complex and Math::MatrixReal gives us the behavior we want for everything except matrix I/O, i.e. parsing and stringification.
<langsyntaxhighlight lang="perl">use strict;
use warnings;
use English;
use Math::Complex;
Line 2,364 ⟶ 2,642:
sub identity {
my $N = shift;
my $m = new Math::MatrixReal->new($N, $N);
$m->one();
return $m;
Line 2,370 ⟶ 2,648:
 
sub example1 {
my $m = new Math::MatrixReal->new(2, 2);
$m->assign(1, 1, cplx(3, 0));
$m->assign(1, 2, cplx(2, 1));
Line 2,379 ⟶ 2,657:
 
sub example2 {
my $m = new Math::MatrixReal->new(3, 3);
$m->assign(1, 1, cplx(1, 0));
$m->assign(1, 2, cplx(1, 0));
Line 2,393 ⟶ 2,671:
 
sub example3 {
my $m = new Math::MatrixReal->new(3, 3);
$m->assign(1, 1, cplx(0.70710677, 0));
$m->assign(1, 2, cplx(0.70710677, 0));
Line 2,404 ⟶ 2,682:
$m->assign(3, 3, cplx(0, 1));
return $m;
}</langsyntaxhighlight>
{{out}}
<pre>
Line 2,443 ⟶ 2,721:
 
=={{header|Phix}}==
Phix has no support for complex numbers, so roll our own, ditto matrix maths. Note this code has no testing for non-square matrices.
<!--<syntaxhighlight lang="phix">(phixonline)-->
<lang Phix>enum REAL, IMAG
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">include</span> <span style="color: #004080;">complex</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
type complex(sequence s)
return length(s)=2 and atom(s[REAL]) and atom(s[IMAG])
<span style="color: #008080;">procedure</span> <span style="color: #000000;">m_print</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
end type
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
function c_add(complex a, complex b)
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
return sq_add(a,b)
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
end function
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">complex_sprint</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
function c_mul(complex a, complex b)
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"["</span><span style="color: #0000FF;">&</span><span style="color: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #008000;">","</span><span style="color: #0000FF;">)&</span><span style="color: #008000;">"]"</span>
return {a[REAL] * b[REAL] - a[IMAG] * b[IMAG],
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
a[REAL] * b[IMAG] + a[IMAG] * b[REAL]}
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)&</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 
function c_conj(complex a)
<span style="color: #008080;">function</span> <span style="color: #000000;">conjugate_transpose</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
return {a[REAL],-a[IMAG]}
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
 
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
function c_print(complex a)
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
if a[IMAG]=0 then return sprintf("%g",a[REAL]) end if
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">complex_conjugate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">][</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
return sprintf("%g%+gi",a)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
procedure m_print(sequence a)
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
integer l = length(a)
for i=1 to l do
<span style="color: #008080;">function</span> <span style="color: #000000;">m_unitary</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">act</span><span style="color: #0000FF;">)</span>
for j=1 to l do
<span style="color: #000080;font-style:italic;">-- note: a was normal and act = a*ct already</span>
a[i][j] = c_print(a[i][j])
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">act</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
a[i] = "["&join(a[i],",")&"]"
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
end for
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">re</span><span style="color: #0000FF;">,</span><span style="color: #000000;">im</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">act</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
puts(1,join(a,"\n")&"\n")
<span style="color: #000080;font-style:italic;">-- round to nearest billionth
end procedure
-- (powers of 2 help the FPU out)</span>
 
<span style="color: #000000;">re</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">round</span><span style="color: #0000FF;">(</span><span style="color: #000000;">re</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1024</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1024</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1024</span><span style="color: #0000FF;">)</span>
 
<span style="color: #000000;">im</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">round</span><span style="color: #0000FF;">(</span><span style="color: #000000;">im</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1024</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1024</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1024</span><span style="color: #0000FF;">)</span>
function conjugate_transpose(sequence a)
<span style="color: #008080;">if</span> <span style="color: #000000;">im</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">0</span>
sequence res = a
<span style="color: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">j</span> <span style="color: #008080;">and</span> <span style="color: #000000;">re</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
integer l = length(a)
<span style="color: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">j</span> <span style="color: #008080;">and</span> <span style="color: #000000;">re</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
for i=1 to l do
<span style="color: #008080;">return</span> <span style="color: #000000;">0</span>
for j=1 to l do
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
res[i][j] = c_conj(a[j][i])
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end for
<span style="color: #008080;">return</span> <span style="color: #000000;">1</span>
return res
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<span style="color: #008080;">function</span> <span style="color: #000000;">m_mul</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>
function m_unitary(sequence act)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
-- note: a was normal and act = a*ct already
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
integer l = length(act)
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
for i=1 to l do
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
for j=1 to l do
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
atom {re,im} = act[i,j]
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">complex_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span><span style="color: #7060A8;">complex_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]))</span>
-- round to nearest billionth
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
-- (powers of 2 help the FPU out)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
re = round(re,1024*1024*1024)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
im = round(im,1024*1024*1024)
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
if im!=0
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
or (i=j and re!=1)
or (i!=j and re!=0) then
<span style="color: #008080;">procedure</span> <span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
return 0
<span style="color: #004080;">sequence</span> <span style="color: #000000;">ct</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">conjugate_transpose</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
end if
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Original matrix:\n"</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #000000;">m_print</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Conjugate transpose:\n"</span><span style="color: #0000FF;">)</span>
return 1
<span style="color: #000000;">m_print</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ct</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #000080;font-style:italic;">-- note: rounding similar to that in m_unitary may be rqd (in a similar
 
-- loop in a new m_equal function) on these two equality tests,
function m_mul(sequence a, sequence b)
-- but as it is, all tests pass with the builtin = operator.</span>
sequence res = sq_mul(a,0)
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Hermitian?: %t\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">=</span><span style="color: #000000;">ct</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (this one)</span>
integer l = length(a)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">act</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ct</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">cta</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ct</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
for i=1 to l do
<span style="color: #004080;">bool</span> <span style="color: #000000;">normal</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">act</span><span style="color: #0000FF;">=</span><span style="color: #000000;">cta</span> <span style="color: #000080;font-style:italic;">-- (&this one)</span>
for j=1 to l do
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Normal?: %t\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">normal</span><span style="color: #0000FF;">)</span>
for k=1 to l do
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Unitary?: %t\n\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">normal</span> <span style="color: #008080;">and</span> <span style="color: #000000;">m_unitary</span><span style="color: #0000FF;">(</span><span style="color: #000000;">act</span><span style="color: #0000FF;">))</span>
res[i][j] = c_add(res[i][j],c_mul(a[i][k],b[k][j]))
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
end for
end for
<span style="color: #008080;">constant</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span>
end for
<span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}},</span>
return res
<span style="color: #0000FF;">{{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}}},</span>
end function
 
<span style="color: #0000FF;">{{{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">}},</span>
procedure test(sequence a)
<span style="color: #0000FF;">{{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span>
sequence ct = conjugate_transpose(a)
<span style="color: #0000FF;">{{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}}},</span>
printf(1,"Original matrix:\n")
m_print(a)
<span style="color: #0000FF;">{{{</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,+</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">}},</span>
printf(1,"Conjugate transpose:\n")
<span style="color: #0000FF;">{{</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,+</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">}}},</span>
m_print(ct)
-- note: rounding similar to that in m_unitary may be rqd (in a similar
<span style="color: #0000FF;">{{{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span>
-- loop in a new m_equal function) on these two equality tests,
<span style="color: #0000FF;">{{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span>
-- but as it is, all tests pass with the builtin = operator.
<span style="color: #0000FF;">{{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}}},</span>
printf(1,"Hermitian?: %s\n",{iff(a=ct?"TRUE":"FALSE")}) -- (this one)
sequence act = m_mul(a,ct), cta = m_mul(ct,a)
<span style="color: #0000FF;">{{{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span>
bool normal = act=cta -- (&this one)
<span style="color: #0000FF;">{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">x</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span>
printf(1,"Normal?: %s\n",{iff(normal?"TRUE":"FALSE")})
<span style="color: #0000FF;">{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">}}},</span>
printf(1,"Unitary?: %s\n\n",{iff(normal and m_unitary(act)?"TRUE":"FALSE")})
end procedure
<span style="color: #0000FF;">{{{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">7</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">}},</span>
 
<span style="color: #0000FF;">{{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">}}}}</span>
constant x = sqrt(2)/2
 
<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">,</span><span style="color: #000000;">test</span><span style="color: #0000FF;">)</span>
constant tests = {{{{3, 0},{2,1}},
<!--</syntaxhighlight>-->
{{2,-1},{1,0}}},
 
{{{ 1, 0},{ 1, 1},{ 0, 2}},
{{ 1,-1},{ 5, 0},{-3, 0}},
{{ 0,-2},{-3, 0},{ 0, 0}}},
 
{{{0.5,+0.5},{0.5,-0.5}},
{{0.5,-0.5},{0.5,+0.5}}},
 
{{{ 1, 0},{ 1, 0},{ 0, 0}},
{{ 0, 0},{ 1, 0},{ 1, 0}},
{{ 1, 0},{ 0, 0},{ 1, 0}}},
 
{{{x, 0},{x, 0},{0, 0}},
{{0,-x},{0, x},{0, 0}},
{{0, 0},{0, 0},{0, 1}}},
 
{{{2,7},{9,-5}},
{{3,4},{8,-6}}}}
 
for i=1 to length(tests) do test(tests[i]) end for</lang>
{{out}}
<pre>
Original matrix:
[3,2+1ii]
[2-1ii,1]
Conjugate transpose:
[3,2+1ii]
[2-1ii,1]
Hermitian?: TRUEtrue
Normal?: TRUEtrue
Unitary?: FALSEfalse
 
Original matrix:
[1,1+1ii,0+2i]
[1-1ii,5,-3]
[0-2i,-3,0]
Conjugate transpose:
[1,1+1ii,0+2i]
[1-1ii,5,-3]
[0-2i,-3,0]
Hermitian?: TRUEtrue
Normal?: TRUEtrue
Unitary?: FALSEfalse
 
Original matrix:
Line 2,594 ⟶ 2,852:
[0.5-0.5i,0.5+0.5i]
[0.5+0.5i,0.5-0.5i]
Hermitian?: FALSEfalse
Normal?: TRUEtrue
Unitary?: TRUEtrue
 
Original matrix:
Line 2,606 ⟶ 2,864:
[1,1,0]
[0,1,1]
Hermitian?: FALSEfalse
Normal?: TRUEtrue
Unitary?: FALSEfalse
 
Original matrix:
[0.707107,0.707107,0]
[0-0.707107i,0+0.707107i,0]
[0,0,0+1ii]
Conjugate transpose:
[0.707107,0+0.707107i,0]
[0.707107,0-0.707107i,0]
[0,0,0-1ii]
Hermitian?: FALSEfalse
Normal?: TRUEtrue
Unitary?: TRUEtrue
 
Original matrix:
Line 2,628 ⟶ 2,886:
[2-7i,3-4i]
[9+5i,8+6i]
Hermitian?: FALSEfalse
Normal?: FALSEfalse
Unitary?: FALSEfalse
</pre>
 
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">
<lang PL/I>
test: procedure options (main); /* 1 October 2012 */
declare n fixed binary;
Line 2,684 ⟶ 2,942:
end MMULT;
end test;
</syntaxhighlight>
</lang>
Outputs from separate runs:
<pre>
Line 2,728 ⟶ 2,986:
 
=={{header|PowerShell}}==
<syntaxhighlight lang="powershell">
<lang PowerShell>
function conjugate-transpose($a) {
$arr = @()
Line 2,805 ⟶ 3,063:
"Normal? `$m = $(are-eq $mhm $hmm)"
"Unitary? `$m = $((are-eq $id2 $hmm) -and (are-eq $id2 $mhm))"
</syntaxhighlight>
</lang>
<b>Output:</b>
<pre>
Line 2,831 ⟶ 3,089:
=={{header|Python}}==
Internally, matrices must be represented as rectangular tuples of tuples of complex numbers.
<langsyntaxhighlight lang="python">def conjugate_transpose(m):
return tuple(tuple(n.conjugate() for n in row) for row in zip(*m))
 
Line 2,902 ⟶ 3,160:
print('Hermitian? %s.' % ishermitian(matrix, ct))
print('Normal? %s.' % isnormal(matrix, ct))
print('Unitary? %s.' % isunitary(matrix, ct))</langsyntaxhighlight>
 
{{out}}
Line 2,940 ⟶ 3,198:
 
=={{header|Racket}}==
<langsyntaxhighlight lang="racket">
#lang racket
(require math)
Line 2,957 ⟶ 3,215:
(define (hermitian? M)
(equal? (H M) M))
</syntaxhighlight>
</lang>
Test:
<langsyntaxhighlight lang="racket">
(define M (matrix [[3.000+0.000i +2.000+1.000i]
[2.000-1.000i +1.000+0.000i]]))
Line 2,966 ⟶ 3,224:
(unitary? M)
(hermitian? M)
</syntaxhighlight>
</lang>
Output:
<langsyntaxhighlight lang="racket">
(array #[#[3.0-0.0i 2.0+1.0i] #[2.0-1.0i 1.0-0.0i]])
#t
#f
#f
</syntaxhighlight>
</lang>
 
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2015-12-13}}
<syntaxhighlight lang="raku" perl6line>for [ # Test Matrices
[ 1, 1+i, 2i],
[ 1-i, 5, -3],
Line 3,029 ⟶ 3,287:
}
 
sub say-it (@array) { $_».fmt("%9s").say for @array }</langsyntaxhighlight>
{{out}}
<pre>Matrix:
Line 3,072 ⟶ 3,330:
 
=={{header|REXX}}==
<langsyntaxhighlight lang="rexx">/*REXX program performs a conjugate transpose on a complex square matrix. */
parse arg N elements; if N==''|N=="," then N=3 /*Not specified? Then use the default.*/
k= 0; do r=1 for N
Line 3,150 ⟶ 3,408:
numeric digits; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g *.5'e'_ % 2
m.=9; do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/; return g</langsyntaxhighlight>
{{out|output|text=&nbsp; when using the default input:}}
<pre>
Line 3,189 ⟶ 3,447:
MMH is Unary: yes
MHM is Unary: yes
</pre>
 
=={{header|RPL}}==
Although basic, RPL's matrix handling capabilities help to keep the code compact but still readable.
{{works with|Halcyon Calc|4.2.7}}
≪ - ABS 1E-10 < ≫
´SAME?´ STO
≪ DUP TRN → m mh
≪ m mh SAME? "Hermitian. " "" IFTE
m mh * mh m * SAME? "Normal. " "" IFTE +
m INV mh SAME? "Unitary. " "" IFTE +
´CNJTRN’ STO
[[(3,0) (2,1)][(2,-1) (1,0)]]
'Hm' STO
[[1 1 0][0 1 1][1 0 1]]
'Nm' STO
≪ [[1 1 0][(0,1) (0,-1) 0][0 0 0]] 2 √ 2 / * {3 3} (0,1) PUT ≫
'Um' STO
Hm CNJTRN
Nm CNJTRN
Um CNJTRN
{{out}}
<pre>
3: "Hermitian. Normal. "
2: "Normal. "
1: "Normal. Unitary. "
</pre>
 
=={{header|Ruby}}==
{{works with|Ruby|2.0}}
<langsyntaxhighlight lang="ruby">require 'matrix'
 
# Start with some matrix.
Line 3,217 ⟶ 3,505:
print ' normal? false'
print ' unitary? false'
end</langsyntaxhighlight>
Note: Ruby 1.9 had a bug in the Matrix#hermitian? method. It's fixed in 2.0.
 
=={{header|Rust}}==
Uses external crate 'num', version 0.1.34
<langsyntaxhighlight lang="rust">
extern crate num; // crate for complex numbers
 
Line 3,302 ⟶ 3,590:
println!("Unitary?: FALSE");
}
}</langsyntaxhighlight>
Output:
<pre>
Line 3,326 ⟶ 3,614:
 
=={{header|Scala}}==
<langsyntaxhighlight Scalalang="scala">object ConjugateTranspose {
case class Complex(re: Double, im: Double) {
Line 3,402 ⟶ 3,690:
}
}</langsyntaxhighlight>
{{out}}
<pre>
Line 3,422 ⟶ 3,710:
=={{header|Sidef}}==
{{trans|Raku}}
<langsyntaxhighlight lang="ruby">func is_Hermitian (Array m, Array t) -> Bool { m == t }
 
func mat_mult (Array a, Array b, Number ε = -3) {
Line 3,487 ⟶ 3,775:
say "Is Normal?\t#{is_Normal(m, t)}"
say "Is Unitary?\t#{is_Unitary(m, t)}"
}</langsyntaxhighlight>
{{out}}
<pre>
Line 3,533 ⟶ 3,821:
Sparkling has support for basic complex algebraic operations, but complex matrix operations are not in the standard library.
 
<langsyntaxhighlight lang="sparkling"># Computes conjugate transpose of M
let conjTransp = function conjTransp(M) {
return map(range(sizeof M[0]), function(row) {
Line 3,633 ⟶ 3,921:
print("U x U* = ");
printCplxMat(cplxMatMul(U, conjTransp(U)));
print();</langsyntaxhighlight>
 
=={{header|Stata}}==
In Mata, the ' operator is always the conjugate transpose. To get only the matrix transpose without complex conjugate, use the [ transposeonly] function.
 
<langsyntaxhighlight lang="stata">
: a=1,2i\3i,4
 
Line 3,670 ⟶ 3,958:
: a'*a==I(rows(a))
0
</syntaxhighlight>
</lang>
 
=={{header|Tcl}}==
Line 3,676 ⟶ 3,964:
{{tcllib|math::complexnumbers}}
{{tcllib|struct::matrix}}
<langsyntaxhighlight lang="tcl">package require struct::matrix
package require math::complexnumbers
 
Line 3,732 ⟶ 4,020:
}
return $mat
}</langsyntaxhighlight>
Using these tools to test for the properties described in the task:
<langsyntaxhighlight lang="tcl">proc isHermitian {matrix {epsilon 1e-14}} {
if {[$matrix rows] != [$matrix columns]} {
# Must be square!
Line 3,785 ⟶ 4,073:
$mmh destroy
return $result
}</langsyntaxhighlight>
<!-- Wot, no demonstration? -->
 
Line 3,796 ⟶ 4,084:
 
However, if we use the ''almostEquals'' method with the default tolerance of 1.0e-14, then we do get a ''true'' result.
<langsyntaxhighlight ecmascriptlang="wren">import "./complex" for Complex, CMatrix
import "./fmt" for Fmt
 
var cm1 = CMatrix.new(
Line 3,829 ⟶ 4,117:
var cm4 = cm3 * cm3.conjTranspose
var id = CMatrix.identity(3)
System.print("Unitary : %(cm4.almostEquals(id))")</langsyntaxhighlight>
 
{{out}}
9,479

edits