C# BigInteger.Sqrt方法代码示例

        public static BigInteger Run(BigInteger N, int B, int M)
        {
            // Find basis for sieving
            var primes = nt.Primes(B).Where(p => nt.IsQuadResidue(p, N)).ToList();
            var sqN = N.Sqrt();
            Func<int, BigInteger> Q = x => (x + sqN) * (x + sqN) - N;

            // use quadsieve on V
            BigInteger[] V = Enumerable.Range(0, M).Select(x => Q(x)).ToArray();
            int[,] coeff; // coeff[i,j] : Q(i) = Prod p_j ^e_ij
            Sieve(N, primes, ref V, out coeff);

            int K = primes.Count;
            int L = V.Count(w => w.IsOne);
            int[] indices = new int[L]; // indices[i] = x_i, Q(x_i) : B-smooth
            int[,] vectors = new int[K, L];
            for(int i=0, cnt=0; i<V.Length && cnt < L; i++)
            {
                if (!V[i].IsOne) continue;

                for(int j=0; j<K; j++)
                {
                    indices[cnt] = i;
                    vectors[j, cnt] = coeff[i, j] % 2; // transposed
                }
                cnt++;
            }

            var exponent = GaussianElimination(vectors);
            BigInteger ans1, ans2;
            ans1 = ans2 = BigInteger.One;
            for (int i = 0; i < exponent.Length; i++)
            {
                if (exponent[i] == 0) continue;
                ans1 *= (indices[i] + sqN);
                ans2 *= Q(indices[i]);
            }
            ans2 = ans2.Sqrt();
            var a = BigInteger.GreatestCommonDivisor(ans1 - ans2, N);
            if (a.IsOne)
                a = BigInteger.GreatestCommonDivisor(ans1 + ans2, N);
            return a;
        }

        /// <summary>  v    cQ1
        ///  Initializes a SievingRequest based on given parameters, and finds the quadratic residues for primes specified
        /// </summary>
        /// <param name="N">The number to factor</param>
        /// <param name="B">The limit for smooth numbers</param>
        /// <param name="f">The polynomial to sieve</param>
        /// <param name="sievereq">The SieveRequest that will be initialized</param>
        public static void InitSievingRequest(BigInteger N, int B, PolynomialFunction f, SieveRequest sievereq)
        {
            sievereq.AStart = (int)N.Sqrt() + 1;
            sievereq.StartIdx = 0;
            sievereq.polyFunction = f;
            sievereq.L = primeSupply[B];

            SievingData sievedat = new SievingData();

            EvaluatePoly(sievereq, sievedat);

            List<List<int>> tmpPrimeStarts = new List<List<int>>();
            List<int> tmpPrimeIntervals = new List<int>();

            for (int pI = 0; pI < B; pI++)
            {
                int p = primeSupply[pI];
                List<int> tmp = new List<int>();
                for (int a = 0; a < p; a++)
                {
                    if (sievedat.V[a] % p == 0)
                    {
                        tmp.Add(a);
                    }
                }
                if (tmp.Count > 0)
                {
                    tmpPrimeIntervals.Add(p);
                    tmpPrimeStarts.Add(tmp);
                }
            }

            sievereq.PrimeIntervals = tmpPrimeIntervals.ToArray();
            sievereq.PrimeStarts = tmpPrimeStarts.ToArray();
            sievereq.B = sievereq.PrimeIntervals.Length;
        }

 private static void Sieve(BigInteger N, List<int> primes, ref BigInteger[] V, out int[,] coeff)
 {
     var sqN = N.Sqrt();
     int M = V.Length;
     int K = primes.Count;
     // exponent vectors
     coeff = new int[M, K];
     for (int k = 0; k < K; k++)
     {
         int p = primes[k];
         int a;
         if (p == 2)
         {
             a = (sqN.IsEven ^ N.IsEven) ? 1 : 0;
             for (int i = a; i < M; i += 2)
             {
                 if (!V[i].IsEven) throw new ArithmeticException();
                 V[i] /= 2;
                 coeff[i, k]++;
             }
             continue;
         }
         int sr = nt.ModSqrt(N, p);
         a = (p+sr - sqN.Mod(p)) % p;
         int b = (2*p - sr - sqN.Mod(p)) % p;
         for (int i = a; i < M; i += p)
         {
             if (V[i] % p != 0) throw new ArithmeticException();
             V[i] /= p;
             coeff[i, k]++;
         }
         for (int i = b; i < M; i += p)
         {
             if (V[i] % p != 0) throw new ArithmeticException();
             V[i] /= p;
             coeff[i, k]++;
         }
     }
 }

        /// <summary>
        /// Checks whether a number is a prime number or not.
        /// </summary>
        /// <param name="n">The number to test.</param>
        /// <param name="checkCache">Whether to check cache for primality.</param>
        /// <returns>Whether the number is prime.</returns>
        public static bool IsPrime(BigInteger n, bool checkCache = true)
        {
            if (n <= 1) return false;
            if (checkCache && n <= BagMath.LargestSmallPrime) return Array.BinarySearch<int>(BagMath.SmallPrimes, (int)n) >= 0;
            if (n < 4) return true;
            if (n % 2 == 0) return false;
            if (n < 9) return true;
            if (n % 3 == 0) return false;

            BigInteger r = n.Sqrt(),
                       f = 5;

            while (f <= r)
            {
                if (n % f == 0) return false;
                if (n % (f + 2) == 0) return false;
                f += 6;
            }

            return true;
        }

        public static void Test_Sqrt_Neg1throwsException()
        {
            var neg1 = new BigInteger(-1);

            var sqrt = neg1.Sqrt();
        }