1 Text format
using System;
namespace Legalsoft.Truffer
{
public class Poissondev : Ran
{
private double lambda { get; set; }
private double sqlam { get; set; }
private double loglam { get; set; }
private double lamexp { get; set; }
private double lambold { get; set; }
private double[] logfact { get; set; }
public Poissondev(double llambda, ulong ix) : base(ix)
{
this.lambda = llambda;
this.logfact = new double[1024];
for (int i = 0; i < 1024; i + + )
{
logfact[i] = -1.0;
}
this.lambold = -1.0;
}
public intdev()
{
double v2 = 0.0;
int k;
if (lambda < 5.0)
{
if (lambda != lambold)
{
lamexp = Math.Exp(-lambda);
}
k = -1;
double t = 1.0;
do
{
+ + k;
t *= doub();
} while (t > lamexp);
}
else
{
if (lambda != lambold)
{
sqlam = Math.Sqrt(lambda);
loglam = Math.Log(lambda);
}
for (; ; )
{
double u = 0.64 * doub();
double v = -0.68 + 1.28 * doub();
if (lambda > 13.5)
{
v2 = Globals.SQR(v);
if (v >= 0.0)
{
if (v2 > 6.5 * u * (0.64 – u) * (u + 0.2))
{
continue;
}
}
else
{
if (v2 > 9.6 * u * (0.66 – u) * (u + 0.07))
{
continue;
}
}
}
k = (int)Math.Floor(sqlam * (v / u) + lambda + 0.5);
if (k < 0)
{
continue;
}
double u2 = Globals.SQR(u);
if (lambda > 13.5)
{
if (v >= 0.0)
{
if (v2 < 15.2 * u2 * (0.61 - u) * (0.8 - u))
{
break;
}
}
else
{
if (v2 < 6.76 * u2 * (0.62 - u) * (1.4 - u))
{
break;
}
}
}
double lfac;
if (k < 1024)
{
if (logfact[k] < 0.0)
{
logfact[k] = Globals.gammln(k + 1.0);
}
lfac = logfact[k];
}
else
{
lfac = Globals.gammln(k + 1.0);
}
double p = sqlam * Math.Exp(-lambda + k * loglam – lfac);
if (u2 < p)
{
break;
}
}
}
lambold = lambda;
return k;
}
public int dev(double llambda)
{
lambda = llambda;
return dev();
}
}
}
2 Code format
using System; namespace Legalsoft.Truffer { public class Poissondev : Ran { private double lambda { get; set; } private double sqlam { get; set; } private double loglam { get; set; } private double lamexp { get; set; } private double lambold { get; set; } private double[] logfact { get; set; } public Poissondev(double llambda, ulong ix) : base(ix) { this.lambda = llambda; this.logfact = new double[1024]; for (int i = 0; i < 1024; i + + ) { logfact[i] = -1.0; } this.lambold = -1.0; } public intdev() { double v2 = 0.0; int k; if (lambda < 5.0) { if (lambda != lambold) { lamexp = Math.Exp(-lambda); } k = -1; double t = 1.0; do { + + k; t *= doub(); } while (t > lamexp); } else { if (lambda != lambold) { sqlam = Math.Sqrt(lambda); loglam = Math.Log(lambda); } for (; ; ) { double u = 0.64 * doub(); double v = -0.68 + 1.28 * doub(); if (lambda > 13.5) { v2 = Globals.SQR(v); if (v >= 0.0) { if (v2 > 6.5 * u * (0.64 - u) * (u + 0.2)) { continue; } } else { if (v2 > 9.6 * u * (0.66 - u) * (u + 0.07)) { continue; } } } k = (int)Math.Floor(sqlam * (v / u) + lambda + 0.5); if (k < 0) { continue; } double u2 = Globals.SQR(u); if (lambda > 13.5) { if (v >= 0.0) { if (v2 < 15.2 * u2 * (0.61 - u) * (0.8 - u)) { break; } } else { if (v2 < 6.76 * u2 * (0.62 - u) * (1.4 - u)) { break; } } } double lfac; if (k < 1024) { if (logfact[k] < 0.0) { logfact[k] = Globals.gammln(k + 1.0); } lfac = logfact[k]; } else { lfac = Globals.gammln(k + 1.0); } double p = sqlam * Math.Exp(-lambda + k * loglam - lfac); if (u2 < p) { break; } } } lambold = lambda; return k; } public int dev(double llambda) { lambda = llambda; return dev(); } } }
The knowledge points of the article match the official knowledge files, and you can further learn related knowledge. Algorithm skill tree Home page Overview 54745 people are learning the system