C# complex类代码示例

 public static void Main()
 {
     var a = new complex(1,2);
     var b = new complex(2,3);
     var c = a+b;
     if(c == new complex(3,5)) c.print("test complex addition passed");
     else c.print("test complex addition failed");
 }

        public static void fftc1d(ref complex[] a) {
            int n;


            n = ap.len(a);
            fft.fftc1d(ref a, n);

            return;
        }

 public void AddTest()
 {
     IComplex compl = new complex();
       complex a = new complex(2, 3);
       complex b = new complex(3, 4);
       complex actual = compl.Add(a, b);
       complex expected = a + b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

 public void SubtractTest()
 {
     IComplex compl = new complex();
       complex a = new complex(20, 30);
       complex b = new complex(30, 40);
       complex actual = compl.Subtract(a, b);
       complex expected = a - b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

 public void MultiplyTest()
 {
     IComplex compl = new complex();
       complex a = new complex(20, 30);
       complex b = new complex(30, 40);
       complex actual = compl.Multiply(a, b);
       complex expected = a * b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

 public void DivideTest()
 {
     IComplex compl = new complex();
       complex a = new complex(20, 30);
       complex b = new complex(30, 40);
       complex actual = compl.Divide(a, b);
       complex expected = a / b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

        private void Test_conj() {
            // test for common matrix 4x3
            ILArray<complex> A = new complex[,] {
                {new complex(1,2), new complex(3,4), new complex(5,6)},
                {new complex(7,8), new complex(9,10), new complex(11,12)},
                {new complex(13,14), new complex(15,16), new complex(17,18)},
                {new complex(19,20), new complex(21,22), new complex(23,24)}};
            ILArray<fcomplex> Af = ILMath.tofcomplex(A); 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af));
            // test on reference matrix
            Test_isConj(A,ILMath.conj(A.R)); 
            Test_isfConj(Af,ILMath.conj(Af.R));
            // test on 2x2x3 matrix (reference)
            A = A.Reshape(new ILDimension(2,2,3)).R;
            if (!A.IsReference) 
                throw new Exception("Unable to test conj() for reference 2x2x3!"); 
            Af = Af.Reshape(new ILDimension(2,2,3)).R; 
            if (!Af.IsReference) 
                throw new Exception("Unable to test conj() for reference 2x2x3!"); 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af));
            // test those on dense array 
            Test_isConj(A,ILMath.conj(A.C)); 
            Test_isfConj(Af,ILMath.conj(Af.C));
            // test for empty
            A = ILArray<complex>.empty();
            Af = ILArray<fcomplex>.empty();
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af)); 
            // test on vector
            A = new complex[] { new complex(1,2), new complex(3,4), new complex(5,6)}; 
            Af = ILMath.tofcomplex(A); 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af)); 
            // transpose 
            Test_isConj(A.T,ILMath.conj(A.T)); 
            Test_isfConj(Af.T,ILMath.conj(Af.T)); 
            // test on scalar 
            A = A[2]; Af = Af[2]; 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af)); 

        }

    /*************************************************************************
    Dense solver. Same as RMatrixMixedSolve(), but for complex matrices.

    Algorithm features:
    * automatic detection of degenerate cases
    * condition number estimation
    * iterative refinement
    * O(N^2) complexity

    INPUT PARAMETERS
        A       -   array[0..N-1,0..N-1], system matrix
        LUA     -   array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
        P       -   array[0..N-1], pivots array, CMatrixLU result
        N       -   size of A
        B       -   array[0..N-1], right part

    OUTPUT PARAMETERS
        Info    -   same as in RMatrixSolveM
        Rep     -   same as in RMatrixSolveM
        X       -   same as in RMatrixSolveM

      -- ALGLIB --
         Copyright 27.01.2010 by Bochkanov Sergey
    *************************************************************************/
    public static void cmatrixmixedsolve(complex[,] a, complex[,] lua, int[] p, int n, complex[] b, out int info, out densesolverreport rep, out complex[] x)
    {
        info = 0;
        rep = new densesolverreport();
        x = new complex[0];
        densesolver.cmatrixmixedsolve(a, lua, p, n, b, ref info, rep.innerobj, ref x);
        return;
    }

 public static void smp_cmatrixsolve(complex[,] a, int n, complex[] b, out int info, out densesolverreport rep, out complex[] x)
 {
     info = 0;
     rep = new densesolverreport();
     x = new complex[0];
     densesolver._pexec_cmatrixsolve(a, n, b, ref info, rep.innerobj, ref x);
     return;
 }

        /*************************************************************************
        Basic LU solver for ScaleA*PLU*x = y.

        This subroutine assumes that:
        * L is well-scaled, and it is U which needs scaling by ScaleA.
        * A=PLU is well-conditioned, so no zero divisions or overflow may occur

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        private static void cbasiclusolve(complex[,] lua,
            int[] p,
            double scalea,
            int n,
            ref complex[] xb,
            ref complex[] tmp)
        {
            int i = 0;
            complex v = 0;
            int i_ = 0;

            for(i=0; i<=n-1; i++)
            {
                if( p[i]!=i )
                {
                    v = xb[i];
                    xb[i] = xb[p[i]];
                    xb[p[i]] = v;
                }
            }
            for(i=1; i<=n-1; i++)
            {
                v = 0.0;
                for(i_=0; i_<=i-1;i_++)
                {
                    v += lua[i,i_]*xb[i_];
                }
                xb[i] = xb[i]-v;
            }
            xb[n-1] = xb[n-1]/(scalea*lua[n-1,n-1]);
            for(i=n-2; i>=0; i--)
            {
                for(i_=i+1; i_<=n-1;i_++)
                {
                    tmp[i_] = scalea*lua[i,i_];
                }
                v = 0.0;
                for(i_=i+1; i_<=n-1;i_++)
                {
                    v += tmp[i_]*xb[i_];
                }
                xb[i] = (xb[i]-v)/(scalea*lua[i,i]);
            }
        }

        /*************************************************************************
        Internal LU solver

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        private static void cmatrixlusolveinternal(complex[,] lua,
            int[] p,
            double scalea,
            int n,
            complex[,] a,
            bool havea,
            complex[,] b,
            int m,
            ref int info,
            densesolverreport rep,
            ref complex[,] x)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            int rfs = 0;
            int nrfs = 0;
            complex[] xc = new complex[0];
            complex[] y = new complex[0];
            complex[] bc = new complex[0];
            complex[] xa = new complex[0];
            complex[] xb = new complex[0];
            complex[] tx = new complex[0];
            double[] tmpbuf = new double[0];
            complex v = 0;
            double verr = 0;
            double mxb = 0;
            double scaleright = 0;
            bool smallerr = new bool();
            bool terminatenexttime = new bool();
            int i_ = 0;

            info = 0;
            x = new complex[0,0];

            alglib.ap.assert((double)(scalea)>(double)(0));
            
            //
            // prepare: check inputs, allocate space...
            //
            if( n<=0 || m<=0 )
            {
                info = -1;
                return;
            }
            for(i=0; i<=n-1; i++)
            {
                if( p[i]>n-1 || p[i]<i )
                {
                    info = -1;
                    return;
                }
            }
            x = new complex[n, m];
            y = new complex[n];
            xc = new complex[n];
            bc = new complex[n];
            tx = new complex[n];
            xa = new complex[n+1];
            xb = new complex[n+1];
            tmpbuf = new double[2*n+2];
            
            //
            // estimate condition number, test for near singularity
            //
            rep.r1 = rcond.cmatrixlurcond1(lua, n);
            rep.rinf = rcond.cmatrixlurcondinf(lua, n);
            if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) || (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
            {
                for(i=0; i<=n-1; i++)
                {
                    for(j=0; j<=m-1; j++)
                    {
                        x[i,j] = 0;
                    }
                }
                rep.r1 = 0;
                rep.rinf = 0;
                info = -3;
                return;
            }
            info = 1;
            
            //
            // solve
            //
            for(k=0; k<=m-1; k++)
            {
                
                //
                // copy B to contiguous storage
                //
                for(i_=0; i_<=n-1;i_++)
                {
//.........这里部分代码省略.........

        /*************************************************************************
        Dense solver. Same as RMatrixLUSolve(), but for  HPD matrices  represented
        by their Cholesky decomposition.

        Algorithm features:
        * automatic detection of degenerate cases
        * O(N^2) complexity
        * condition number estimation
        * matrix is represented by its upper or lower triangle

        No iterative refinement is provided because such partial representation of
        matrix does not allow efficient calculation of extra-precise  matrix-vector
        products for large matrices. Use RMatrixSolve or RMatrixMixedSolve  if  you
        need iterative refinement.

        INPUT PARAMETERS
            CHA     -   array[0..N-1,0..N-1], Cholesky decomposition,
                        SPDMatrixCholesky result
            N       -   size of A
            IsUpper -   what half of CHA is provided
            B       -   array[0..N-1], right part

        OUTPUT PARAMETERS
            Info    -   same as in RMatrixSolve
            Rep     -   same as in RMatrixSolve
            X       -   same as in RMatrixSolve

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void hpdmatrixcholeskysolve(complex[,] cha,
            int n,
            bool isupper,
            complex[] b,
            ref int info,
            densesolverreport rep,
            ref complex[] x)
        {
            complex[,] bm = new complex[0,0];
            complex[,] xm = new complex[0,0];
            int i_ = 0;

            info = 0;
            x = new complex[0];

            if( n<=0 )
            {
                info = -1;
                return;
            }
            bm = new complex[n, 1];
            for(i_=0; i_<=n-1;i_++)
            {
                bm[i_,0] = b[i_];
            }
            hpdmatrixcholeskysolvem(cha, n, isupper, bm, 1, ref info, rep, ref xm);
            x = new complex[n];
            for(i_=0; i_<=n-1;i_++)
            {
                x[i_] = xm[i_,0];
            }
        }

 /*************************************************************************
 Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
 *************************************************************************/
 public static void _pexec_hpdmatrixsolve(complex[,] a,
     int n,
     bool isupper,
     complex[] b,
     ref int info,
     densesolverreport rep,
     ref complex[] x)
 {
     hpdmatrixsolve(a,n,isupper,b,ref info,rep,ref x);
 }

        /*************************************************************************
        1-dimensional real inverse FFT.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            F   -   array[0..floor(N/2)] - frequencies from forward real FFT
            N   -   problem size

        OUTPUT PARAMETERS
            A   -   inverse DFT of a input array, array[0..N-1]

        NOTE:
            F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
        half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
        is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
        F[floor(N/2)] has no special properties.

        Relying on properties noted above, FFTR1DInv subroutine uses only elements
        from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
        N is even it ignores imaginary part of F[floor(N/2)] too.

        When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
        - you can pass either either frequencies array with N elements or  reduced
        array with roughly N/2 elements - subroutine will  successfully  transform
        both.

        If you call this function using reduced arguments list -  "FFTR1DInv(F,A)"
        - you must pass FULL array with N elements (although higher  N/2 are still
        not used) because array size is used to automatically determine FFT length


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftr1dinv(complex[] f,
            int n,
            ref double[] a)
        {
            int i = 0;
            double[] h = new double[0];
            complex[] fh = new complex[0];

            a = new double[0];

            ap.assert(n>0, "FFTR1DInv: incorrect N!");
            ap.assert(ap.len(f)>=(int)Math.Floor((double)n/(double)2)+1, "FFTR1DInv: Length(F)<Floor(N/2)+1!");
            ap.assert(math.isfinite(f[0].x), "FFTR1DInv: F contains infinite or NAN values!");
            for(i=1; i<=(int)Math.Floor((double)n/(double)2)-1; i++)
            {
                ap.assert(math.isfinite(f[i].x) & math.isfinite(f[i].y), "FFTR1DInv: F contains infinite or NAN values!");
            }
            ap.assert(math.isfinite(f[(int)Math.Floor((double)n/(double)2)].x), "FFTR1DInv: F contains infinite or NAN values!");
            if( n%2!=0 )
            {
                ap.assert(math.isfinite(f[(int)Math.Floor((double)n/(double)2)].y), "FFTR1DInv: F contains infinite or NAN values!");
            }
            
            //
            // Special case: N=1, FFT is just identity transform.
            // After this block we assume that N is strictly greater than 1.
            //
            if( n==1 )
            {
                a = new double[1];
                a[0] = f[0].x;
                return;
            }
            
            //
            // inverse real FFT is reduced to the inverse real FHT,
            // which is reduced to the forward real FHT,
            // which is reduced to the forward real FFT.
            //
            // Don't worry, it is really compact and efficient reduction :)
            //
            h = new double[n];
            a = new double[n];
            h[0] = f[0].x;
            for(i=1; i<=(int)Math.Floor((double)n/(double)2)-1; i++)
            {
                h[i] = f[i].x-f[i].y;
                h[n-i] = f[i].x+f[i].y;
            }
            if( n%2==0 )
            {
                h[(int)Math.Floor((double)n/(double)2)] = f[(int)Math.Floor((double)n/(double)2)].x;
            }
            else
            {
                h[(int)Math.Floor((double)n/(double)2)] = f[(int)Math.Floor((double)n/(double)2)].x-f[(int)Math.Floor((double)n/(double)2)].y;
                h[(int)Math.Floor((double)n/(double)2)+1] = f[(int)Math.Floor((double)n/(double)2)].x+f[(int)Math.Floor((double)n/(double)2)].y;
            }
            fftr1d(h, n, ref fh);
            for(i=0; i<=n-1; i++)
            {
                a[i] = (fh[i].x-fh[i].y)/n;
            }
        }

    /*************************************************************************
    1-dimensional complex inverse FFT.

    Array size N may be arbitrary number (composite or prime).  Algorithm  has
    O(N*logN) complexity for any N (composite or prime).

    See FFTC1D() description for more information about algorithm performance.

    INPUT PARAMETERS
        A   -   array[0..N-1] - complex array to be transformed
        N   -   problem size

    OUTPUT PARAMETERS
        A   -   inverse DFT of a input array, array[0..N-1]
                A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)


      -- ALGLIB --
         Copyright 29.05.2009 by Bochkanov Sergey
    *************************************************************************/
    public static void fftc1dinv(ref complex[] a, int n)
    {

        fft.fftc1dinv(ref a, n);
        return;
    }