[2806] | 1 | /*************************************************************************
|
---|
| 2 | AP library
|
---|
| 3 | Copyright (c) 2003-2009 Sergey Bochkanov (ALGLIB project).
|
---|
| 4 |
|
---|
| 5 | >>> LICENSE >>>
|
---|
| 6 | This program is free software; you can redistribute it and/or modify
|
---|
| 7 | it under the terms of the GNU General Public License as published by
|
---|
| 8 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 9 | License, or (at your option) any later version.
|
---|
| 10 |
|
---|
| 11 | This program is distributed in the hope that it will be useful,
|
---|
| 12 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 14 | GNU General Public License for more details.
|
---|
| 15 |
|
---|
| 16 | A copy of the GNU General Public License is available at
|
---|
| 17 | http://www.fsf.org/licensing/licenses
|
---|
| 18 |
|
---|
| 19 | >>> END OF LICENSE >>>
|
---|
| 20 | *************************************************************************/
|
---|
| 21 | namespace AP
|
---|
| 22 | {
|
---|
| 23 | /********************************************************************
|
---|
| 24 | Class defining a complex number with double precision.
|
---|
| 25 | ********************************************************************/
|
---|
| 26 | public struct Complex
|
---|
| 27 | {
|
---|
| 28 | public double x;
|
---|
| 29 | public double y;
|
---|
| 30 |
|
---|
| 31 | public Complex(double _x)
|
---|
| 32 | {
|
---|
| 33 | x = _x;
|
---|
| 34 | y = 0;
|
---|
| 35 | }
|
---|
| 36 | public Complex(double _x, double _y)
|
---|
| 37 | {
|
---|
| 38 | x = _x;
|
---|
| 39 | y = _y;
|
---|
| 40 | }
|
---|
| 41 | public static implicit operator Complex(double _x)
|
---|
| 42 | {
|
---|
| 43 | return new Complex(_x);
|
---|
| 44 | }
|
---|
| 45 | public static bool operator==(Complex lhs, Complex rhs)
|
---|
| 46 | {
|
---|
| 47 | return ((double)lhs.x==(double)rhs.x) & ((double)lhs.y==(double)rhs.y);
|
---|
| 48 | }
|
---|
| 49 | public static bool operator!=(Complex lhs, Complex rhs)
|
---|
| 50 | {
|
---|
| 51 | return ((double)lhs.x!=(double)rhs.x) | ((double)lhs.y!=(double)rhs.y);
|
---|
| 52 | }
|
---|
| 53 | public static Complex operator+(Complex lhs)
|
---|
| 54 | {
|
---|
| 55 | return lhs;
|
---|
| 56 | }
|
---|
| 57 | public static Complex operator-(Complex lhs)
|
---|
| 58 | {
|
---|
| 59 | return new Complex(-lhs.x,-lhs.y);
|
---|
| 60 | }
|
---|
| 61 | public static Complex operator+(Complex lhs, Complex rhs)
|
---|
| 62 | {
|
---|
| 63 | return new Complex(lhs.x+rhs.x,lhs.y+rhs.y);
|
---|
| 64 | }
|
---|
| 65 | public static Complex operator-(Complex lhs, Complex rhs)
|
---|
| 66 | {
|
---|
| 67 | return new Complex(lhs.x-rhs.x,lhs.y-rhs.y);
|
---|
| 68 | }
|
---|
| 69 | public static Complex operator*(Complex lhs, Complex rhs)
|
---|
| 70 | {
|
---|
| 71 | return new Complex(lhs.x*rhs.x-lhs.y*rhs.y, lhs.x*rhs.y+lhs.y*rhs.x);
|
---|
| 72 | }
|
---|
| 73 | public static Complex operator/(Complex lhs, Complex rhs)
|
---|
| 74 | {
|
---|
| 75 | Complex result;
|
---|
| 76 | double e;
|
---|
| 77 | double f;
|
---|
| 78 | if( System.Math.Abs(rhs.y)<System.Math.Abs(rhs.x) )
|
---|
| 79 | {
|
---|
| 80 | e = rhs.y/rhs.x;
|
---|
| 81 | f = rhs.x+rhs.y*e;
|
---|
| 82 | result.x = (lhs.x+lhs.y*e)/f;
|
---|
| 83 | result.y = (lhs.y-lhs.x*e)/f;
|
---|
| 84 | }
|
---|
| 85 | else
|
---|
| 86 | {
|
---|
| 87 | e = rhs.x/rhs.y;
|
---|
| 88 | f = rhs.y+rhs.x*e;
|
---|
| 89 | result.x = (lhs.y+lhs.x*e)/f;
|
---|
| 90 | result.y = (-lhs.x+lhs.y*e)/f;
|
---|
| 91 | }
|
---|
| 92 | return result;
|
---|
| 93 | }
|
---|
| 94 | public override int GetHashCode()
|
---|
| 95 | {
|
---|
| 96 | return x.GetHashCode() ^ y.GetHashCode();
|
---|
| 97 | }
|
---|
| 98 | public override bool Equals(object obj)
|
---|
| 99 | {
|
---|
| 100 | if( obj is byte)
|
---|
| 101 | return Equals(new Complex((byte)obj));
|
---|
| 102 | if( obj is sbyte)
|
---|
| 103 | return Equals(new Complex((sbyte)obj));
|
---|
| 104 | if( obj is short)
|
---|
| 105 | return Equals(new Complex((short)obj));
|
---|
| 106 | if( obj is ushort)
|
---|
| 107 | return Equals(new Complex((ushort)obj));
|
---|
| 108 | if( obj is int)
|
---|
| 109 | return Equals(new Complex((int)obj));
|
---|
| 110 | if( obj is uint)
|
---|
| 111 | return Equals(new Complex((uint)obj));
|
---|
| 112 | if( obj is long)
|
---|
| 113 | return Equals(new Complex((long)obj));
|
---|
| 114 | if( obj is ulong)
|
---|
| 115 | return Equals(new Complex((ulong)obj));
|
---|
| 116 | if( obj is float)
|
---|
| 117 | return Equals(new Complex((float)obj));
|
---|
| 118 | if( obj is double)
|
---|
| 119 | return Equals(new Complex((double)obj));
|
---|
| 120 | if( obj is decimal)
|
---|
| 121 | return Equals(new Complex((double)(decimal)obj));
|
---|
| 122 | return base.Equals(obj);
|
---|
| 123 | }
|
---|
| 124 | }
|
---|
| 125 |
|
---|
| 126 | /********************************************************************
|
---|
| 127 | AP math namespace
|
---|
| 128 | ********************************************************************/
|
---|
| 129 | public struct rcommstate
|
---|
| 130 | {
|
---|
| 131 | public int stage;
|
---|
| 132 | public int[] ia;
|
---|
| 133 | public bool[] ba;
|
---|
| 134 | public double[] ra;
|
---|
| 135 | public AP.Complex[] ca;
|
---|
| 136 | };
|
---|
| 137 |
|
---|
| 138 | /********************************************************************
|
---|
| 139 | AP math namespace
|
---|
| 140 | ********************************************************************/
|
---|
| 141 | public class Math
|
---|
| 142 | {
|
---|
| 143 | //public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond);
|
---|
| 144 | public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond + 1000*System.DateTime.Now.Second + 60*1000*System.DateTime.Now.Minute);
|
---|
| 145 |
|
---|
| 146 | public const double MachineEpsilon = 5E-16;
|
---|
| 147 | public const double MaxRealNumber = 1E300;
|
---|
| 148 | public const double MinRealNumber = 1E-300;
|
---|
| 149 |
|
---|
| 150 | public static double RandomReal()
|
---|
| 151 | {
|
---|
| 152 | double r = 0;
|
---|
| 153 | lock(RndObject){ r = RndObject.NextDouble(); }
|
---|
| 154 | return r;
|
---|
| 155 | }
|
---|
| 156 | public static int RandomInteger(int N)
|
---|
| 157 | {
|
---|
| 158 | int r = 0;
|
---|
| 159 | lock(RndObject){ r = RndObject.Next(N); }
|
---|
| 160 | return r;
|
---|
| 161 | }
|
---|
| 162 | public static double Sqr(double X)
|
---|
| 163 | {
|
---|
| 164 | return X*X;
|
---|
| 165 | }
|
---|
| 166 | public static double AbsComplex(Complex z)
|
---|
| 167 | {
|
---|
| 168 | double w;
|
---|
| 169 | double xabs;
|
---|
| 170 | double yabs;
|
---|
| 171 | double v;
|
---|
| 172 |
|
---|
| 173 | xabs = System.Math.Abs(z.x);
|
---|
| 174 | yabs = System.Math.Abs(z.y);
|
---|
| 175 | w = xabs>yabs ? xabs : yabs;
|
---|
| 176 | v = xabs<yabs ? xabs : yabs;
|
---|
| 177 | if( v==0 )
|
---|
| 178 | return w;
|
---|
| 179 | else
|
---|
| 180 | {
|
---|
| 181 | double t = v/w;
|
---|
| 182 | return w*System.Math.Sqrt(1+t*t);
|
---|
| 183 | }
|
---|
| 184 | }
|
---|
| 185 | public static Complex Conj(Complex z)
|
---|
| 186 | {
|
---|
| 187 | return new Complex(z.x, -z.y);
|
---|
| 188 | }
|
---|
| 189 | public static Complex CSqr(Complex z)
|
---|
| 190 | {
|
---|
| 191 | return new Complex(z.x*z.x-z.y*z.y, 2*z.x*z.y);
|
---|
| 192 | }
|
---|
| 193 |
|
---|
| 194 | }
|
---|
| 195 | }
|
---|