complex_algorithm.cs

Numerical algorithm for Complex numbers

        public struct Complex
        {
            public double real;
            public double imaginary;

            public Complex(double real, double imaginary)  //constructor
            {
                this.real = real;
                this.imaginary = imaginary;
            }

            // Declare which operator to overload (+),
            // the types that can be added (two Complex objects),
            // and the return type (Complex):
            public static Complex operator +(Complex c1, Complex c2)
            {
                return new Complex(c1.real + c2.real, c1.imaginary + c2.imaginary);
            }
            public static Complex operator +(double a , Complex c)
            {
                return new Complex(a + c.real, c.imaginary);
            }

            public static Complex operator -(Complex c1, Complex c2)
            {
                return new Complex(c1.real - c2.real, c1.imaginary - c2.imaginary);
            }
            public static Complex operator -(double a , Complex c)
            {
                return new Complex(a - c.real, c.imaginary );
            }

            public static Complex operator *(Complex c1, Complex c2)
            {
                double ac = c1.real * c2.real;
                double bd = c1.imaginary * c2.imaginary;
                return new Complex(ac - bd, (c1.real + c1.imaginary) * (c2.real + c2.imaginary) - ac - bd);
            }

            public static Complex operator *(double a, Complex c)
            {
                return new Complex(a * c.real, a * c.imaginary);
            }

            public static Complex operator /(Complex c1, Complex c2)
            {
                double dc = 0.0, cdc = 1.0;
                if (Math.Abs(c2.real) >= Math.Abs(c2.imaginary))
                {
                    dc = c2.imaginary / c2.real;
                    cdc = c2.real + c2.imaginary * dc;
                    return new Complex((c1.real + c1.imaginary * dc) / cdc, (c1.imaginary - c1.real * dc) / cdc);
                }
                else
                {
                    dc = c2.real / c2.imaginary;
                    cdc = c2.real * dc + c2.imaginary;
                    return new Complex((c1.real * dc + c1.imaginary) / cdc, (c1.imaginary * dc - c1.real) / cdc);
                }
            }

            public Complex conjg(Complex c)
            {
                return new Complex(c.real, -c.imaginary);
            }

            public double cabs(Complex c)
            {
                double lr = Math.Max(Math.Abs(c.real), Math.Abs(c.imaginary));
                double sr = Math.Min(Math.Abs(c.real), Math.Abs(c.imaginary));
                double ca = lr * Math.Sqrt(1.0 + (sr * sr) / (lr * lr));
                return ca;
            }

            public Complex csqrt(Complex c)
            {
                double a = c.real;
                double b = c.imaginary;
                double ab = 0.0, C = 0.0, D = 0.0;
                if (a == 0.0)
                {
                    if (b == 0)
                    {
                        return new Complex(0.0, 0.0);
                    }
                    else
                    {
                        C = Math.Sqrt(Math.Abs(b) / 2.0);

                        if (b > 0)
                        {
                            return new Complex(C, C);
                        }
                        else
                        {
                            return new Complex(C, -C);
                        }
                    }
                }
                else
                {
                    if (a > 0)
                    {
                        if (Math.Abs(a) >= Math.Abs(b))
                        {
                            ab = b / a;
                            C = Math.Abs(a) * Math.Sqrt((1.0 + Math.Sqrt(1.0 + ab * ab)) / 2.0);
                        }
                        else
                        {
                            ab = a / b;
                            C = Math.Abs(b) * Math.Sqrt((Math.Abs(ab) + Math.Sqrt(1.0 + ab * ab)) / 2.0);
                        }

                        D = b / (2.0 * C);
                        return new Complex(C, D);
                    }
                    else
                    {
                        if (Math.Abs(a) >= Math.Abs(b))
                        {
                            ab = b / a;
                            D = Math.Sign(b) * Math.Abs(a) * Math.Sqrt((1.0 + Math.Sqrt(1.0 + ab * ab)) / 2.0);
                        }
                        else
                        {
                            ab = a / b;
                            D = Math.Sign(b) * Math.Abs(b) * Math.Sqrt((Math.Abs(ab) + Math.Sqrt(1.0 + ab * ab)) / 2.0);
                        }

                        C = Math.Abs(b) / (2.0 * D);
                        return new Complex(C, D);
                    }

                }

            }

            public Complex cexp(Complex c)
            {
                double r = 0.0, p = 0.0;
                r = Math.Exp(c.real);
                p = c.imaginary;
                return new Complex(r * Math.Cos(p), r * Math.Sin(p));
            }

            public static Complex ctanh(Complex c)
            {
                Complex C = new Complex(0.0, 0.0);
               // Complex cm = new Complex(-1.0, 0.0);

                C = (cexp(c) - cexp(-c)) / (cexp(c) + cexp(-c));
                return C;
            }
        }