C# Rational.Subtract方法代码示例

        public static BigDecimal Zeta(int n, MathContext mc)
        {
            if (n <= 0)
                throw new NotSupportedException("Zeta at negative argument " + n + " not supported");
            if (n == 1)
                throw new ArithmeticException("Pole at zeta(1) ");

            if (n%2 == 0) {
                /* Even indices. Abramowitz-Stegun 23.2.16. Start with 2^(n-1)*B(n)/n!
                        */
                Rational b = Bernoulli.Default[n].Abs();
                b = b.Divide(Factorial.Default[n]);
                b = b.Multiply(BigInteger.One.ShiftLeft(n - 1));

                /* to be multiplied by pi^n. Absolute error in the result of pi^n is n times
                        * error in pi times pi^(n-1). Relative error is n*error(pi)/pi, requested by mc.
                        * Need one more digit in pi if n=10, two digits if n=100 etc, and add one extra digit.
                        */
                var mcpi = new MathContext(mc.Precision + (int) (System.Math.Log10(10.0*n)));
                BigDecimal piton = PiRound(mcpi).Pow(n, mc);
                return MultiplyRound(piton, b);
            }
            if (n == 3) {
                /* Broadhurst BBP <a href="http://arxiv.org/abs/math/9803067">arXiv:math/9803067</a>
                        * Error propagation: S31 is roughly 0.087, S33 roughly 0.131
                        */
                int[] a31 = {1, -7, -1, 10, -1, -7, 1, 0};
                int[] a33 = {1, 1, -1, -2, -1, 1, 1, 0};
                BigDecimal S31 = BroadhurstBbp(3, 1, a31, mc);
                BigDecimal S33 = BroadhurstBbp(3, 3, a33, mc);
                S31 = S31.Multiply(new BigDecimal(48));
                S33 = S33.Multiply(new BigDecimal(32));
                return S31.Add(S33).Divide(new BigDecimal(7), mc);
            }
            if (n == 5) {
                /* Broadhurst BBP <a href=http://arxiv.org/abs/math/9803067">arXiv:math/9803067</a>
                        * Error propagation: S51 is roughly -11.15, S53 roughly 22.165, S55 is roughly 0.031
                        * 9*2048*S51/6265 = -3.28. 7*2038*S53/61651= 5.07. 738*2048*S55/61651= 0.747.
                        * The result is of the order 1.03, so we add 2 digits to S51 and S52 and one digit to S55.
                        */
                int[] a51 = {31, -1614, -31, -6212, -31, -1614, 31, 74552};
                int[] a53 = {173, 284, -173, -457, -173, 284, 173, -111};
                int[] a55 = {1, 0, -1, -1, -1, 0, 1, 1};
                BigDecimal S51 = BroadhurstBbp(5, 1, a51, new MathContext(2 + mc.Precision));
                BigDecimal S53 = BroadhurstBbp(5, 3, a53, new MathContext(2 + mc.Precision));
                BigDecimal S55 = BroadhurstBbp(5, 5, a55, new MathContext(1 + mc.Precision));
                S51 = S51.Multiply(new BigDecimal(18432));
                S53 = S53.Multiply(new BigDecimal(14336));
                S55 = S55.Multiply(new BigDecimal(1511424));
                return S51.Add(S53).Subtract(S55).Divide(new BigDecimal(62651), mc);
            }
            /* Cohen et al Exp Math 1 (1) (1992) 25
                        */
            var betsum = new Rational();
            var bern = new Bernoulli();
            var fact = new Factorial();
            for (int npr = 0; npr <= (n + 1)/2; npr++) {
                Rational b = bern[2*npr].Multiply(bern[n + 1 - 2*npr]);
                b = b.Divide(fact[2*npr]).Divide(fact[n + 1 - 2*npr]);
                b = b.Multiply(1 - 2*npr);
                if (npr%2 == 0)
                    betsum = betsum.Add(b);
                else
                    betsum = betsum.Subtract(b);
            }
            betsum = betsum.Divide(n - 1);
            /* The first term, including the facor (2pi)^n, is essentially most
                        * of the result, near one. The second term below is roughly in the range 0.003 to 0.009.
                        * So the precision here is matching the precisionn requested by mc, and the precision
                        * requested for 2*pi is in absolute terms adjusted.
                        */
            var mcloc = new MathContext(2 + mc.Precision + (int) (System.Math.Log10(n)));
            BigDecimal ftrm = PiRound(mcloc).Multiply(new BigDecimal(2));
            ftrm = ftrm.Pow(n);
            ftrm = MultiplyRound(ftrm, betsum.ToBigDecimal(mcloc));
            var exps = new BigDecimal(0);

            /* the basic accuracy of the accumulated terms before multiplication with 2
                        */
            double eps = System.Math.Pow(10d, -mc.Precision);

            if (n%4 == 3) {
                /* since the argument n is at least 7 here, the drop
                                * of the terms is at rather constant pace at least 10^-3, for example
                                * 0.0018, 0.2e-7, 0.29e-11, 0.74e-15 etc for npr=1,2,3.... We want 2 times these terms
                                * fall below eps/10.
                                */
                int kmax = mc.Precision/3;
                eps /= kmax;
                /* need an error of eps for 2/(exp(2pi)-1) = 0.0037
                                * The absolute error is 4*exp(2pi)*err(pi)/(exp(2pi)-1)^2=0.0075*err(pi)
                                */
                BigDecimal exp2p = PiRound(new MathContext(3 + ErrorToPrecision(3.14, eps/0.0075)));
                exp2p = Exp(exp2p.Multiply(new BigDecimal(2)));
                BigDecimal c = exp2p.Subtract(BigDecimal.One);
                exps = DivideRound(1, c);
                for (int npr = 2; npr <= kmax; npr++) {
                    /* the error estimate above for npr=1 is the worst case of
                                        * the absolute error created by an error in 2pi. So we can
                                        * safely re-use the exp2p value computed above without
//.........这里部分代码省略.........