1 文本格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// 二项式偏差
    /// Binomial Deviates
    /// </summary>
    public class Binomialdev : Ran
    {
        private double pp { get; set; }
        private double p { get; set; }
        private double pb { get; set; }
        private double expnp { get; set; }
        private double np { get; set; }
        private double glnp { get; set; }
        private double plog { get; set; }
        private double pclog { get; set; }
        private double sq { get; set; }
        private int n { get; set; }
        private int swch { get; set; }
        private ulong uz { get; set; }
        private ulong uo { get; set; }
        private ulong unfin { get; set; }
        private ulong diff { get; set; }
        private ulong rltp { get; set; }
        private int[] pbits { get; set; } = new int[5];
        private double[] cdf { get; set; } = new double[64];
        private double[] logfact { get; set; } = new double[1024];

        public Binomialdev(int nn, double ppp, ulong i) : base(i)
        {
            this.pp = ppp;
            this.n = nn;
            pb = p = (pp <= 0.5 ? pp : 1.0 - pp);
            if (n <= 64)
            {
                uz = 0;
                uo = 0xffffffffffffffffL;
                rltp = 0;
                for (int j = 0; j < 5; j++)
                {
                    pbits[j] = 1 & ((int)(pb *= 2.0));
                }
                pb -= Math.Floor(pb);
                swch = 0;
            }
            else if (n * p < 30.0)
            {
                cdf[0] = Math.Exp(n * Math.Log(1 - p));
                for (int j = 1; j < 64; j++)
                {
                    cdf[j] = cdf[j - 1] + Math.Exp(Globals.gammln(n + 1.0) - Globals.gammln(j + 1.0) - Globals.gammln(n - j + 1.0) + j * Math.Log(p) + (n - j) * Math.Log(1.0 - p));
                }
                swch = 1;
            }
            else
            {
                np = n * p;
                glnp = Globals.gammln(n + 1.0);
                plog = Math.Log(p);
                pclog = Math.Log(1.0 - p);
                sq = Math.Sqrt(np * (1.0 - p));
                if (n < 1024)
                {
                    for (int j = 0; j <= n; j++)
                    {
                        logfact[j] = Globals.gammln(j + 1.0);
                    }
                }
                swch = 2;
            }
        }

        public int dev()
        {
            int k;
            if (swch == 0)
            {
                unfin = uo;
                for (int j = 0; j < 5; j++)
                {
                    diff = unfin & (int64() ^ (pbits[j] > 0 ? uo : uz));
                    if (pbits[j] > 0)
                    {
                        rltp |= diff;
                    }
                    else
                    {
                        rltp = rltp & ~diff;
                    }
                    unfin = unfin & ~diff;
                }
                k = 0;
                for (int j = 0; j < n; j++)
                {
                    if ((unfin & 1) > 0)
                    {
                        if (doub() < pb)
                        {
                            ++k;
                        }
                    }
                    else
                    {
                        if ((rltp & 1) > 0)
                        {
                            ++k;
                        }
                    }
                    unfin >>= 1;
                    rltp >>= 1;
                }
            }
            else if (swch == 1)
            {
                double y = doub();
                int kl = -1;
                k = 64;
                while (k - kl > 1)
                {
                    int km = (kl + k) / 2;
                    if (y < cdf[km])
                    {
                        k = km;
                    }
                    else
                    {
                        kl = km;
                    }
                }
            }
            else
            {
                for (; ; )
                {
                    double u = 0.645 * doub();
                    double v = -0.63 + 1.25 * doub();
                    double v2 = Globals.SQR(v);
                    if (v >= 0.0)
                    {
                        if (v2 > 6.5 * u * (0.645 - u) * (u + 0.2))
                        {
                            continue;
                        }
                    }
                    else
                    {
                        if (v2 > 8.4 * u * (0.645 - u) * (u + 0.1))
                        {
                            continue;
                        }
                    }
                    k = (int)(Math.Floor(sq * (v / u) + np + 0.5));
                    if (k < 0 || k > n)
                    {
                        continue;
                    }
                    double u2 = Globals.SQR(u);
                    if (v >= 0.0)
                    {
                        if (v2 < 12.25 * u2 * (0.615 - u) * (0.92 - u))
                        {
                            break;
                        }
                    }
                    else
                    {
                        if (v2 < 7.84 * u2 * (0.615 - u) * (1.2 - u))
                        {
                            break;
                        }
                    }
                    double b = sq * Math.Exp(glnp + k * plog + (n - k) * pclog - (n < 1024 ? logfact[k] + logfact[n - k] : Globals.gammln(k + 1.0) + Globals.gammln(n - k + 1.0)));
                    if (u2 < b)
                    {
                        break;
                    }
                }
            }
            if (p != pp)
            {
                k = n - k;
            }
            return k;
        }
    }
}
 

2 代码格式

using System;namespace Legalsoft.Truffer
{/// <summary>/// 二项式偏差/// Binomial Deviates/// </summary>public class Binomialdev : Ran{private double pp { get; set; }private double p { get; set; }private double pb { get; set; }private double expnp { get; set; }private double np { get; set; }private double glnp { get; set; }private double plog { get; set; }private double pclog { get; set; }private double sq { get; set; }private int n { get; set; }private int swch { get; set; }private ulong uz { get; set; }private ulong uo { get; set; }private ulong unfin { get; set; }private ulong diff { get; set; }private ulong rltp { get; set; }private int[] pbits { get; set; } = new int[5];private double[] cdf { get; set; } = new double[64];private double[] logfact { get; set; } = new double[1024];public Binomialdev(int nn, double ppp, ulong i) : base(i){this.pp = ppp;this.n = nn;pb = p = (pp <= 0.5 ? pp : 1.0 - pp);if (n <= 64){uz = 0;uo = 0xffffffffffffffffL;rltp = 0;for (int j = 0; j < 5; j++){pbits[j] = 1 & ((int)(pb *= 2.0));}pb -= Math.Floor(pb);swch = 0;}else if (n * p < 30.0){cdf[0] = Math.Exp(n * Math.Log(1 - p));for (int j = 1; j < 64; j++){cdf[j] = cdf[j - 1] + Math.Exp(Globals.gammln(n + 1.0) - Globals.gammln(j + 1.0) - Globals.gammln(n - j + 1.0) + j * Math.Log(p) + (n - j) * Math.Log(1.0 - p));}swch = 1;}else{np = n * p;glnp = Globals.gammln(n + 1.0);plog = Math.Log(p);pclog = Math.Log(1.0 - p);sq = Math.Sqrt(np * (1.0 - p));if (n < 1024){for (int j = 0; j <= n; j++){logfact[j] = Globals.gammln(j + 1.0);}}swch = 2;}}public int dev(){int k;if (swch == 0){unfin = uo;for (int j = 0; j < 5; j++){diff = unfin & (int64() ^ (pbits[j] > 0 ? uo : uz));if (pbits[j] > 0){rltp |= diff;}else{rltp = rltp & ~diff;}unfin = unfin & ~diff;}k = 0;for (int j = 0; j < n; j++){if ((unfin & 1) > 0){if (doub() < pb){++k;}}else{if ((rltp & 1) > 0){++k;}}unfin >>= 1;rltp >>= 1;}}else if (swch == 1){double y = doub();int kl = -1;k = 64;while (k - kl > 1){int km = (kl + k) / 2;if (y < cdf[km]){k = km;}else{kl = km;}}}else{for (; ; ){double u = 0.645 * doub();double v = -0.63 + 1.25 * doub();double v2 = Globals.SQR(v);if (v >= 0.0){if (v2 > 6.5 * u * (0.645 - u) * (u + 0.2)){continue;}}else{if (v2 > 8.4 * u * (0.645 - u) * (u + 0.1)){continue;}}k = (int)(Math.Floor(sq * (v / u) + np + 0.5));if (k < 0 || k > n){continue;}double u2 = Globals.SQR(u);if (v >= 0.0){if (v2 < 12.25 * u2 * (0.615 - u) * (0.92 - u)){break;}}else{if (v2 < 7.84 * u2 * (0.615 - u) * (1.2 - u)){break;}}double b = sq * Math.Exp(glnp + k * plog + (n - k) * pclog - (n < 1024 ? logfact[k] + logfact[n - k] : Globals.gammln(k + 1.0) + Globals.gammln(n - k + 1.0)));if (u2 < b){break;}}}if (p != pp){k = n - k;}return k;}}
}