private long SampleUsingBTPE()

in Standard/src/Emulation/Math.cs [115:230]

• 98 lines of code
• 20 McCabe index (conditional complexity)

        private long SampleUsingBTPE()
        {
            // The following code is ported from NumPy's binomial
            // sampling implementation.
            double r,q,fm,p1,xm,xl,xr,c,laml,lamr,p2,p3,p4;
            double a,u,v,s,F,rho,t,A,nrq,x1,x2,f1,f2,z,z2,w,w2,x;
            long m,y,k,i;

            /* initialize */
            r = Min(SuccessProbability, FailureProbability);
            q = 1.0 - r;
            fm = NSamples * r + r;
            m = (long)Floor(fm);
            p1 = Floor(2.195 * Sqrt(NSamples * r * q) - 4.6 * q) + 0.5;
            xm = m + 0.5;
            xl = xm - p1;
            xr = xm + p1;
            c = 0.134 + 20.5/(15.3 + m);
            a = (fm - xl)/(fm-xl*r);
            laml = a*(1.0 + a/2.0);
            a = (xr - fm)/(xr*q);
            lamr = a*(1.0 + a/2.0);
            p2 = p1*(1.0 + 2.0*c);
            p3 = p2 + c/laml;
            p4 = p3 + c/lamr;


        /* sigh ... */
        Step10:
            nrq = NSamples * r * q;
            u = RandomState.NextDouble() * p4;
            v = RandomState.NextDouble();
            if (u > p1) goto Step20;
            y = (long)Floor(xm - p1*v + u);
            goto Step60;

        Step20:
            if (u > p2) goto Step30;
            x = xl + (u - p1)/c;
            v = v*c + 1.0 - Abs(m - x + 0.5)/p1;
            if (v > 1.0) goto Step10;
            y = (long)Floor(x);
            goto Step50;

        Step30:
            if (u > p3) goto Step40;
            y = (long)Floor(xl + Log(v)/laml);
            /* Reject if v == 0.0 since cast of inf not well defined */
            if ((y < 0) || (v == 0.0)) goto Step10;
            v = v*(u-p2)*laml;
            goto Step50;

        Step40:
            y = (long)Floor(xr - Log(v)/lamr);
            /* Reject if v == 0.0 since cast of inf not well defined */
            if ((y > NSamples) || (v == 0.0)) goto Step10;
            v = v*(u-p3)*lamr;

        Step50:
            k = Abs(y - m);
            if ((k > 20) && (k < ((nrq)/2.0 - 1))) goto Step52;

            s = r/q;
            a = s*(NSamples+1);
            F = 1.0;
            if (m < y)
            {
                for (i=m+1; i<=y; i++)
                {
                    F *= (a/i - s);
                }
            }
            else if (m > y)
            {
                for (i=y+1; i<=m; i++)
                {
                    F /= (a/i - s);
                }
            }
            if (v > F) goto Step10;
            goto Step60;

        Step52:
            rho = (k/(nrq))*((k*(k/3.0 + 0.625) + 0.16666666666666666)/nrq + 0.5);
            t = -k*k/(2*nrq); // lgtm [cs/loss-of-precision]
            A = Log(v);
            if (A < (t - rho)) goto Step60;
            if (A > (t + rho)) goto Step10;

            x1 = y+1;
            f1 = m+1;
            z = NSamples + 1-m;
            w = NSamples - y+1;
            x2 = x1*x1;
            f2 = f1*f1;
            z2 = z*z;
            w2 = w*w;
            if (A > (xm*Log(f1/x1)
                + (NSamples - m+0.5)*Log(z/w)
                + (y-m)*Log(w*r/(x1*q))
                + (13680.0-(462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0
                + (13680.0-(462.0-(132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0
                + (13680.0-(462.0-(132.0-(99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0
                + (13680.0-(462.0-(132.0-(99.0-140.0/w2)/w2)/w2)/w2)/w/166320.0))
            {
                goto Step10;
            }

        Step60:
            if (SuccessProbability > 0.5)
            {
                y = NSamples - y;
            }

            return y;
        }