Imported Upstream version 4.0.0~alpha1

Former-commit-id: 806294f5ded97629b74c85c09952f2a74fe182d9

external/referencesource/System.Numerics/System/Numerics/Complex.cs

@@ -0,0 +1,413 @@
+// ==++==
+// 
+//   Copyright (c) Microsoft Corporation.  All rights reserved.
+// 
+// ==--==
+/*=========================================================================
+**
+** Class: Complex
+**
+**
+** Purpose: 
+** This feature is intended to create Complex Number as a type 
+** that can be a part of the .NET framework (base class libraries).  
+** A complex number z is a number of the form z = x + yi, where x and y 
+** are real numbers, and i is the imaginary unit, with the property i2= -1.
+**
+**
+===========================================================================*/
+
+using System;
+using System.Globalization;
+using System.Diagnostics.CodeAnalysis;
+using System.Diagnostics.Contracts;
+
+namespace System.Numerics {
+
+#if !SILVERLIGHT
+    [Serializable]
+#endif // !SILVERLIGHT
+    public struct Complex : IEquatable<Complex>, IFormattable {
+
+        // --------------SECTION: Private Data members ----------- //
+
+        private Double m_real;
+        private Double m_imaginary;
+
+        // ---------------SECTION: Necessary Constants ----------- //
+
+        private const Double LOG_10_INV = 0.43429448190325;
+
+
+        // --------------SECTION: Public Properties -------------- //
+
+        public Double Real {
+            get {
+                return m_real;
+            }
+        }
+
+        public Double Imaginary {
+            get {
+                return m_imaginary;
+            }
+        }
+
+        public Double Magnitude {
+            get {
+                return Complex.Abs(this);
+            }
+        }
+
+        public Double Phase {
+            get {
+                return Math.Atan2(m_imaginary, m_real);
+            }
+        }
+
+        // --------------SECTION: Attributes -------------- //
+
+        public static readonly Complex Zero = new Complex(0.0, 0.0);
+        public static readonly Complex One = new Complex(1.0, 0.0);
+        public static readonly Complex ImaginaryOne = new Complex(0.0, 1.0);
+
+        // --------------SECTION: Constructors and factory methods -------------- //
+
+        public Complex(Double real, Double imaginary)  /* Constructor to create a complex number with rectangular co-ordinates  */
+        {
+            this.m_real = real;
+            this.m_imaginary = imaginary;
+        }
+
+        public static Complex FromPolarCoordinates(Double magnitude, Double phase) /* Factory method to take polar inputs and create a Complex object */
+        {
+            return new Complex((magnitude * Math.Cos(phase)), (magnitude * Math.Sin(phase)));
+        }
+
+        public static Complex Negate(Complex value) {
+            return -value;
+        }
+
+        public static Complex Add(Complex left, Complex right) {
+            return left + right;
+        }
+
+        public static Complex Subtract(Complex left, Complex right) {
+            return left - right;
+        }
+
+        public static Complex Multiply(Complex left, Complex right) {
+            return left * right;
+        }
+
+        public static Complex Divide(Complex dividend, Complex divisor) {
+            return dividend / divisor;
+        }
+
+        // --------------SECTION: Arithmetic Operator(unary) Overloading -------------- //
+        public static Complex operator -(Complex value)  /* Unary negation of a complex number */
+        {
+
+            return (new Complex((-value.m_real), (-value.m_imaginary)));
+        }
+
+        // --------------SECTION: Arithmetic Operator(binary) Overloading -------------- //       
+        public static Complex operator +(Complex left, Complex right) {
+            return (new Complex((left.m_real + right.m_real), (left.m_imaginary + right.m_imaginary)));
+
+        }
+
+        public static Complex operator -(Complex left, Complex right) {
+            return (new Complex((left.m_real - right.m_real), (left.m_imaginary - right.m_imaginary)));
+        }
+
+        public static Complex operator *(Complex left, Complex right) {
+            // Multiplication:  (a + bi)(c + di) = (ac -bd) + (bc + ad)i
+            Double result_Realpart = (left.m_real * right.m_real) - (left.m_imaginary * right.m_imaginary);
+            Double result_Imaginarypart = (left.m_imaginary * right.m_real) + (left.m_real * right.m_imaginary);
+            return (new Complex(result_Realpart, result_Imaginarypart));
+        }
+
+        public static Complex operator /(Complex left, Complex right) {
+            // Division : Smith's formula.
+            double a = left.m_real;
+            double b = left.m_imaginary;
+            double c = right.m_real;
+            double d = right.m_imaginary;
+
+            if (Math.Abs(d) < Math.Abs(c)) {
+                double doc = d / c;
+                return new Complex((a + b * doc) / (c + d * doc), (b - a * doc) / (c + d * doc));
+            } else {
+                double cod = c / d;
+                return new Complex((b + a * cod) / (d + c * cod), (-a + b * cod) / (d + c * cod));
+            }
+        }
+
+
+        // --------------SECTION: Other arithmetic operations  -------------- //
+
+        public static Double Abs(Complex value) {
+
+            if(Double.IsInfinity(value.m_real) || Double.IsInfinity(value.m_imaginary)) {
+                return double.PositiveInfinity;
+            }
+
+            // |value| == sqrt(a^2 + b^2)
+            // sqrt(a^2 + b^2) == a/a * sqrt(a^2 + b^2) = a * sqrt(a^2/a^2 + b^2/a^2)
+            // Using the above we can factor out the square of the larger component to dodge overflow.
+
+
+            double c = Math.Abs(value.m_real);
+            double d = Math.Abs(value.m_imaginary);
+
+            if (c > d) {
+                double r = d / c;
+                return c * Math.Sqrt(1.0 + r * r);
+            } else if (d == 0.0) {
+                return c;  // c is either 0.0 or NaN
+            } else {
+                double r = c / d;
+                return d * Math.Sqrt(1.0 + r * r);
+            }
+        }
+        public static Complex Conjugate(Complex value) {
+            // Conjugate of a Complex number: the conjugate of x+i*y is x-i*y 
+
+            return (new Complex(value.m_real, (-value.m_imaginary)));
+
+        }
+        public static Complex Reciprocal(Complex value) {
+            // Reciprocal of a Complex number : the reciprocal of x+i*y is 1/(x+i*y)
+            if ((value.m_real == 0) && (value.m_imaginary == 0)) {
+                return Complex.Zero;
+            }
+
+            return Complex.One / value;
+        }
+
+        // --------------SECTION: Comparison Operator(binary) Overloading -------------- //
+
+        public static bool operator ==(Complex left, Complex right) {
+            return ((left.m_real == right.m_real) && (left.m_imaginary == right.m_imaginary));
+
+
+        }
+        public static bool operator !=(Complex left, Complex right) {
+            return ((left.m_real != right.m_real) || (left.m_imaginary != right.m_imaginary));
+
+        }
+
+        // --------------SECTION: Comparison operations (methods implementing IEquatable<ComplexNumber>,IComparable<ComplexNumber>) -------------- //
+
+        public override bool Equals(object obj) {
+            if (!(obj is Complex)) return false;
+            return this == ((Complex)obj);
+        }
+        public bool Equals(Complex value) {
+            return ((this.m_real.Equals(value.m_real)) && (this.m_imaginary.Equals(value.m_imaginary)));
+
+        }
+
+        // --------------SECTION: Type-casting basic numeric data-types to ComplexNumber  -------------- //
+
+        public static implicit operator Complex(Int16 value) {
+            return (new Complex(value, 0.0));
+        }
+        public static implicit operator Complex(Int32 value) {
+            return (new Complex(value, 0.0));
+        }
+        public static implicit operator Complex(Int64 value) {
+            return (new Complex(value, 0.0));
+        }
+        [CLSCompliant(false)]
+        public static implicit operator Complex(UInt16 value) {
+            return (new Complex(value, 0.0));
+        }
+        [CLSCompliant(false)]
+        public static implicit operator Complex(UInt32 value) {
+            return (new Complex(value, 0.0));
+        }
+        [CLSCompliant(false)]
+        public static implicit operator Complex(UInt64 value) {
+            return (new Complex(value, 0.0));
+        }
+        [CLSCompliant(false)]
+        public static implicit operator Complex(SByte value) {
+            return (new Complex(value, 0.0));
+        }
+        public static implicit operator Complex(Byte value) {
+            return (new Complex(value, 0.0));
+        }
+        public static implicit operator Complex(Single value) {
+            return (new Complex(value, 0.0));
+        }
+        public static implicit operator Complex(Double value) {
+            return (new Complex(value, 0.0));
+        }
+        public static explicit operator Complex(BigInteger value) {
+            return (new Complex((Double)value, 0.0));
+        }
+        public static explicit operator Complex(Decimal value) {
+            return (new Complex((Double)value, 0.0));
+        }
+
+
+        // --------------SECTION: Formattig/Parsing options  -------------- //
+
+        public override String ToString() {
+            return (String.Format(CultureInfo.CurrentCulture, "({0}, {1})", this.m_real, this.m_imaginary));
+        }
+
+        public String ToString(String format) {
+            return (String.Format(CultureInfo.CurrentCulture, "({0}, {1})", this.m_real.ToString(format, CultureInfo.CurrentCulture), this.m_imaginary.ToString(format, CultureInfo.CurrentCulture)));
+        }
+
+        public String ToString(IFormatProvider provider) {
+            return (String.Format(provider, "({0}, {1})", this.m_real, this.m_imaginary));
+        }
+
+        public String ToString(String format, IFormatProvider provider) {
+            return (String.Format(provider, "({0}, {1})", this.m_real.ToString(format, provider), this.m_imaginary.ToString(format, provider)));
+        }
+
+
+        public override Int32 GetHashCode() {
+            Int32 n1 = 99999997;
+            Int32 hash_real = this.m_real.GetHashCode() % n1;
+            Int32 hash_imaginary = this.m_imaginary.GetHashCode();
+            Int32 final_hashcode = hash_real ^ hash_imaginary;
+            return (final_hashcode);
+        }
+
+
+
+        // --------------SECTION: Trigonometric operations (methods implementing ITrigonometric)  -------------- //
+
+        public static Complex Sin(Complex value) {
+            double a = value.m_real;
+            double b = value.m_imaginary;
+            return new Complex(Math.Sin(a) * Math.Cosh(b), Math.Cos(a) * Math.Sinh(b));
+        }
+
+        [SuppressMessage("Microsoft.Naming", "CA1704:IdentifiersShouldBeSpelledCorrectly", MessageId = "Sinh", Justification = "[....]: Existing Name")]
+        public static Complex Sinh(Complex value) /* Hyperbolic sin */
+        {
+            double a = value.m_real;
+            double b = value.m_imaginary;
+            return new Complex(Math.Sinh(a) * Math.Cos(b), Math.Cosh(a) * Math.Sin(b));
+
+        }
+        public static Complex Asin(Complex value) /* Arcsin */
+        {
+            return (-ImaginaryOne) * Log(ImaginaryOne * value + Sqrt(One - value * value));
+        }
+
+        public static Complex Cos(Complex value) {
+            double a = value.m_real;
+            double b = value.m_imaginary;
+            return new Complex(Math.Cos(a) * Math.Cosh(b), - (Math.Sin(a) * Math.Sinh(b)));
+        }
+
+        [SuppressMessage("Microsoft.Naming", "CA1704:IdentifiersShouldBeSpelledCorrectly", MessageId = "Cosh", Justification = "[....]: Existing Name")]
+        public static Complex Cosh(Complex value) /* Hyperbolic cos */
+        {
+            double a = value.m_real;
+            double b = value.m_imaginary;
+            return new Complex(Math.Cosh(a) * Math.Cos(b), Math.Sinh(a) * Math.Sin(b));
+        }
+        public static Complex Acos(Complex value) /* Arccos */
+        {
+            return (-ImaginaryOne) * Log(value + ImaginaryOne*Sqrt(One - (value * value)));
+
+        }
+        public static Complex Tan(Complex value) {
+            return (Sin(value) / Cos(value));
+        }
+
+        [SuppressMessage("Microsoft.Naming", "CA1704:IdentifiersShouldBeSpelledCorrectly", MessageId = "Tanh", Justification = "[....]: Existing Name")]
+        public static Complex Tanh(Complex value) /* Hyperbolic tan */
+        {
+            return (Sinh(value) / Cosh(value));
+        }
+        public static Complex Atan(Complex value) /* Arctan */
+        {
+            Complex Two = new Complex(2.0, 0.0);
+            return (ImaginaryOne / Two) * (Log(One - ImaginaryOne * value) - Log(One + ImaginaryOne * value));
+        }
+
+        // --------------SECTION: Other numerical functions  -------------- //        
+
+        public static Complex Log(Complex value) /* Log of the complex number value to the base of 'e' */
+        {
+            return (new Complex((Math.Log(Abs(value))), (Math.Atan2(value.m_imaginary, value.m_real))));
+
+        }
+        public static Complex Log(Complex value, Double baseValue) /* Log of the complex number to a the base of a double */
+        {
+            return (Log(value) / Log(baseValue));
+        }
+        public static Complex Log10(Complex value) /* Log to the base of 10 of the complex number */
+        {
+
+            Complex temp_log = Log(value);
+            return (Scale(temp_log, (Double)LOG_10_INV));
+
+        }
+        public static Complex Exp(Complex value) /* The complex number raised to e */
+        {
+            Double temp_factor = Math.Exp(value.m_real);
+            Double result_re = temp_factor * Math.Cos(value.m_imaginary);
+            Double result_im = temp_factor * Math.Sin(value.m_imaginary);
+            return (new Complex(result_re, result_im));
+        }
+
+        [SuppressMessage("Microsoft.Naming", "CA1704:IdentifiersShouldBeSpelledCorrectly", MessageId = "Sqrt", Justification = "[....]: Existing Name")]
+        public static Complex Sqrt(Complex value) /* Square root ot the complex number */
+        {
+            return Complex.FromPolarCoordinates(Math.Sqrt(value.Magnitude), value.Phase / 2.0);
+        }
+
+        public static Complex Pow(Complex value, Complex power) /* A complex number raised to another complex number */
+        {
+
+            if (power == Complex.Zero) {
+                return Complex.One;
+            }
+
+            if (value == Complex.Zero) {
+                return Complex.Zero;
+            }
+
+            double a = value.m_real;
+            double b = value.m_imaginary;
+            double c = power.m_real;
+            double d = power.m_imaginary;
+
+            double rho = Complex.Abs(value);
+            double theta = Math.Atan2(b, a);
+            double newRho = c * theta + d * Math.Log(rho);
+
+            double t = Math.Pow(rho, c) * Math.Pow(Math.E, -d * theta);
+
+            return new Complex(t * Math.Cos(newRho), t * Math.Sin(newRho));
+        }
+
+        public static Complex Pow(Complex value, Double power) // A complex number raised to a real number 
+        {
+            return Pow(value, new Complex(power, 0));
+        }
+
+
+
+        //--------------- SECTION: Private member functions for internal use -----------------------------------//
+
+        private static Complex Scale(Complex value, Double factor) {
+
+            Double result_re = factor * value.m_real;
+            Double result_im = factor * value.m_imaginary;
+            return (new Complex(result_re, result_im));
+        }
+
+    }
+}

external/referencesource/System.Numerics/System/Numerics/NumericsHelpers.cs

@@ -0,0 +1,360 @@
+// ==++==
+// 
+//   Copyright (c) Microsoft Corporation.  All rights reserved.
+// 
+// ==--==
+
+using System;
+using System.Diagnostics;
+using System.Diagnostics.Contracts;
+using System.Collections.Generic;
+using System.Runtime.InteropServices;
+using System.Text;
+
+namespace System.Numerics {
+
+    [StructLayout(LayoutKind.Explicit)]
+    internal struct DoubleUlong {
+      [FieldOffset(0)]
+      public double dbl;
+      [FieldOffset(0)]
+      public ulong uu;
+    }
+
+
+    internal static class NumericsHelpers {
+    private const int kcbitUint = 32;
+
+    public static void GetDoubleParts(double dbl, out int sign, out int exp, out ulong man, out bool fFinite) {
+      Contract.Ensures(Contract.ValueAtReturn(out sign) == +1 || Contract.ValueAtReturn(out sign) == -1);
+
+      DoubleUlong du;
+      du.uu = 0;
+      du.dbl = dbl;
+
+      sign = 1 - ((int)(du.uu >> 62) & 2);
+      man = du.uu & 0x000FFFFFFFFFFFFF;
+      exp = (int)(du.uu >> 52) & 0x7FF;
+      if (exp == 0) {
+        // Denormalized number.
+        fFinite = true;
+        if (man != 0)
+          exp = -1074;
+      }
+      else if (exp == 0x7FF) {
+        // NaN or Inifite.
+        fFinite = false;
+        exp = int.MaxValue;
+      }
+      else {
+        fFinite = true;
+        man |= 0x0010000000000000;
+        exp -= 1075;
+      }
+    }
+
+    public static double GetDoubleFromParts(int sign, int exp, ulong man) {
+      DoubleUlong du;
+      du.dbl = 0;
+
+      if (man == 0)
+        du.uu = 0;
+      else {
+        // Normalize so that 0x0010 0000 0000 0000 is the highest bit set.
+        int cbitShift = CbitHighZero(man) - 11;
+        if (cbitShift < 0)
+          man >>= -cbitShift;
+        else
+          man <<= cbitShift;
+        exp -= cbitShift;
+        Contract.Assert((man & 0xFFF0000000000000) == 0x0010000000000000);
+
+        // Move the point to just behind the leading 1: 0x001.0 0000 0000 0000
+        // (52 bits) and skew the exponent (by 0x3FF == 1023).
+        exp += 1075;
+
+        if (exp >= 0x7FF) {
+          // Infinity.
+          du.uu = 0x7FF0000000000000;
+        }
+        else if (exp <= 0) {
+          // Denormalized.
+          exp--;
+          if (exp < -52) {
+            // Underflow to zero.
+            du.uu = 0;
+          }
+          else {
+            du.uu = man >> -exp;
+            Contract.Assert(du.uu != 0);
+          }
+        }
+        else {
+          // Mask off the implicit high bit.
+          du.uu = (man & 0x000FFFFFFFFFFFFF) | ((ulong)exp << 52);
+        }
+      }
+
+      if (sign < 0)
+        du.uu |= 0x8000000000000000;
+
+      return du.dbl;
+    }
+
+
+
+    // Do an in-place twos complement of d and also return the result.
+    // "Dangerous" because it causes a mutation and needs to be used
+    // with care for immutable types
+    public static uint[] DangerousMakeTwosComplement(uint[] d) {
+        // first do complement and +1 as long as carry is needed
+        int i = 0;
+        uint v = 0;
+        for (; i < d.Length; i++) {
+            v = ~d[i] + 1;
+            d[i] = v;
+            if (v != 0) { i++; break; }
+        }
+        if (v != 0) {
+            // now ones complement is sufficient
+            for (; i < d.Length; i++) {
+                d[i] = ~d[i];
+            }
+        } else {
+            //??? this is weird
+            d = resize(d, d.Length + 1);
+            d[d.Length - 1] = 1;
+        }
+        return d;
+    }
+
+    public static uint[] resize(uint[] v, int len) {
+        if (v.Length == len) return v;
+        uint[] ret = new uint[len];
+        int n = System.Math.Min(v.Length, len);
+        for (int i = 0; i < n; i++) {
+            ret[i] = v[i];
+        }
+        return ret;
+    }
+
+    public static void Swap<T>(ref T a, ref T b) {
+      T tmp = a;
+      a = b;
+      b = tmp;
+    }
+
+    public static uint GCD(uint u1, uint u2) {
+      const int cvMax = 32;
+      if (u1 < u2)
+        goto LOther;
+    LTop:
+      Contract.Assert(u2 <= u1);
+      if (u2 == 0)
+        return u1;
+      for (int cv = cvMax; ; ) {
+        u1 -= u2;
+        if (u1 < u2)
+          break;
+        if (--cv == 0) {
+          u1 %= u2;
+          break;
+        }
+      }
+    LOther:
+      Contract.Assert(u1 < u2);
+      if (u1 == 0)
+        return u2;
+      for (int cv = cvMax; ; ) {
+        u2 -= u1;
+        if (u2 < u1)
+          break;
+        if (--cv == 0) {
+          u2 %= u1;
+          break;
+        }
+      }
+      goto LTop;
+    }
+
+    public static ulong GCD(ulong uu1, ulong uu2) {
+      const int cvMax = 32;
+      if (uu1 < uu2)
+        goto LOther;
+    LTop:
+      Contract.Assert(uu2 <= uu1);
+      if (uu1 <= uint.MaxValue)
+        goto LSmall;
+      if (uu2 == 0)
+        return uu1;
+      for (int cv = cvMax; ; ) {
+        uu1 -= uu2;
+        if (uu1 < uu2)
+          break;
+        if (--cv == 0) {
+          uu1 %= uu2;
+          break;
+        }
+      }
+    LOther:
+      Contract.Assert(uu1 < uu2);
+      if (uu2 <= uint.MaxValue)
+        goto LSmall;
+      if (uu1 == 0)
+        return uu2;
+      for (int cv = cvMax; ; ) {
+        uu2 -= uu1;
+        if (uu2 < uu1)
+          break;
+        if (--cv == 0) {
+          uu2 %= uu1;
+          break;
+        }
+      }
+      goto LTop;
+
+    LSmall:
+      uint u1 = (uint)uu1;
+      uint u2 = (uint)uu2;
+      if (u1 < u2)
+        goto LOtherSmall;
+    LTopSmall:
+      Contract.Assert(u2 <= u1);
+      if (u2 == 0)
+        return u1;
+      for (int cv = cvMax; ; ) {
+        u1 -= u2;
+        if (u1 < u2)
+          break;
+        if (--cv == 0) {
+          u1 %= u2;
+          break;
+        }
+      }
+    LOtherSmall:
+      Contract.Assert(u1 < u2);
+      if (u1 == 0)
+        return u2;
+      for (int cv = cvMax; ; ) {
+        u2 -= u1;
+        if (u2 < u1)
+          break;
+        if (--cv == 0) {
+          u2 %= u1;
+          break;
+        }
+      }
+      goto LTopSmall;
+    }
+
+    public static ulong MakeUlong(uint uHi, uint uLo) {
+      return ((ulong)uHi << kcbitUint) | uLo;
+    }
+
+    public static uint GetLo(ulong uu) {
+      return (uint)uu;
+    }
+
+    public static uint GetHi(ulong uu) {
+      return (uint)(uu >> kcbitUint);
+    }
+
+    public static uint Abs(int a) {
+      uint mask = (uint)(a >> 31);
+      return ((uint)a ^ mask) - mask;
+    }
+
+//    public static ulong Abs(long a) {
+//      ulong mask = (ulong)(a >> 63);
+//      return ((ulong)a ^ mask) - mask;
+//    }
+
+    public static uint CombineHash(uint u1, uint u2) {
+      return ((u1 << 7) | (u1 >> 25)) ^ u2;
+    }
+
+    public static int CombineHash(int n1, int n2) {
+      return (int)CombineHash((uint)n1, (uint)n2);
+    }
+    public static int CbitHighZero(uint u) {
+      if (u == 0)
+        return 32;
+
+      int cbit = 0;
+      if ((u & 0xFFFF0000) == 0) {
+        cbit += 16;
+        u <<= 16;
+      }
+      if ((u & 0xFF000000) == 0) {
+        cbit += 8;
+        u <<= 8;
+      }
+      if ((u & 0xF0000000) == 0) {
+        cbit += 4;
+        u <<= 4;
+      }
+      if ((u & 0xC0000000) == 0) {
+        cbit += 2;
+        u <<= 2;
+      }
+      if ((u & 0x80000000) == 0)
+        cbit += 1;
+      return cbit;
+    }
+
+    public static int CbitLowZero(uint u) {
+      if (u == 0)
+        return 32;
+
+      int cbit = 0;
+      if ((u & 0x0000FFFF) == 0) {
+        cbit += 16;
+        u >>= 16;
+      }
+      if ((u & 0x000000FF) == 0) {
+        cbit += 8;
+        u >>= 8;
+      }
+      if ((u & 0x0000000F) == 0) {
+        cbit += 4;
+        u >>= 4;
+      }
+      if ((u & 0x00000003) == 0) {
+        cbit += 2;
+        u >>= 2;
+      }
+      if ((u & 0x00000001) == 0)
+        cbit += 1;
+      return cbit;
+    }
+
+    public static int CbitHighZero(ulong uu) {
+      if ((uu & 0xFFFFFFFF00000000) == 0)
+        return 32 + CbitHighZero((uint)uu);
+      return CbitHighZero((uint)(uu >> 32));
+    }
+
+//    public static int CbitLowZero(ulong uu) {
+//      if ((uint)uu == 0)
+//        return 32 + CbitLowZero((uint)(uu >> 32));
+//      return CbitLowZero((uint)uu);
+//    }
+//
+//    public static int Cbit(uint u) {
+//      u = (u & 0x55555555) + ((u >> 1) & 0x55555555);
+//      u = (u & 0x33333333) + ((u >> 2) & 0x33333333);
+//      u = (u & 0x0F0F0F0F) + ((u >> 4) & 0x0F0F0F0F);
+//      u = (u & 0x00FF00FF) + ((u >> 8) & 0x00FF00FF);
+//      return (int)((ushort)u + (ushort)(u >> 16));
+//    }
+//
+//    static int Cbit(ulong uu) {
+//      uu = (uu & 0x5555555555555555) + ((uu >> 1) & 0x5555555555555555);
+//      uu = (uu & 0x3333333333333333) + ((uu >> 2) & 0x3333333333333333);
+//      uu = (uu & 0x0F0F0F0F0F0F0F0F) + ((uu >> 4) & 0x0F0F0F0F0F0F0F0F);
+//      uu = (uu & 0x00FF00FF00FF00FF) + ((uu >> 8) & 0x00FF00FF00FF00FF);
+//      uu = (uu & 0x0000FFFF0000FFFF) + ((uu >> 16) & 0x0000FFFF0000FFFF);
+//      return (int)((uint)uu + (uint)(uu >> 32));
+//    }
+  }
+}