📄 random.c
字号:
}
else if( g + v + a * x < 0.0 )
{
y = g + v + a * x + Rangauss(Rand);
accept = (y > 0.0) ? -2.0 * v * y : 0.0;
}
else
{
y = Rangauss(Rand);
accept = (g + a * x) * y - v * fabs(y);
}
} while( accept < 0.0 && accept < log(Randouble(Rand)) );
// Exit
*x0 = x;
*y0 = y;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Posneg2A
//
// Purpose: Set bi-exponential sampler for Posneg2 in case a > 0, g >= v
//
// | exp - x^2 / 2 + f x - u x + (g - v + a x)^2 / 2
// x = 0 |
// | exp - x^2 / 2 + f x + u x + (g - v + a x)^2 / 2
// x2 = -(g-v)/a |
// | exp - x^2 / 2 + f x + u x
// x1 = -(g+v)/a |
// | exp - x^2 / 2 + f x + u x + (g + v + a x)^2 / 2
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Posneg2A(
double fpu, // I f + u > 0
double fmu, // I f - u <= gmv
double gpv, // I g + v > 0
double gmv, // I g - v >= 0
double a, // I > 0
double height, // I height of sampler cusp above target maximum
double* xc, // O target maximum
double* p0, // O sampler exp( yc + p0 * (x-xc) ), x < xc (p0 > 0)
double* q0) // O sampler exp( yc + q0 * (x-xc) ), x > xc (q0 < 0)
{
double x1, x2, y0, y1, y2, m0, m1, m2, m9, xcap, ycap;
double d, p, q, r, s, t, u;
d = 1.0 - a * a;
if( d < DBL_EPSILON )
d = DBL_EPSILON;
// Set boundaries ...(x1,y1)...(x2,y2)...(0,y0)...
// and slopes m1 m2 m9,m0
x1 = -gpv / a;
y1 = x1 * (fpu - 0.5 * x1);
m1 = fpu - x1;
x2 = -gmv / a;
y2 = x2 * (fpu - 0.5 * x2);
m2 = fpu - x2;
y0 = 0.5 * gmv * gmv;
m9 = fpu + a * gmv; // +ve
m0 = fmu + a * gmv;
// Find apex from slope=0
if( m0 <= 0.0 ) // apex at 0
{
xcap = 0.0;
ycap = height + y0;
if(ycap <= y2 - m2*x2) // Lslope(xcap,ycap,-d/2,fpu+a*gmv,gmv*gmv/2)
p = m9 + sqrt(2.0 * d * height);
else
{
r = (ycap - y1) + m1 * x1;
if( r <= 0.0 ) // Lslope(xcap, ycap, -1/2., fpu, 0)
p = fpu + sqrt(2.0 * ycap);
else // Lslope(xcap,ycap,-d/2,fpu+a*gpv,gpv*gpv/2)
p = (fpu + a * gpv) + sqrt(2.0 * d * (d * ycap + a * a * r));
} // Rslope(xcap,ycap,-d/2,fmu+a*gmv,gmv*gmv/2)
q = m0 - sqrt(2.0 * d * height);
}
else // apex in (0,inf)
{
xcap = m0 / d;
ycap = height + 0.5 * (d * xcap * xcap + gmv * gmv);
if(ycap - y0 <= m0*xcap) // Lslope(xcap,ycap,-d/2,fmu+a*gmv,gmv*gmv/2)
p = sqrt(2.0 * d * height);
else
{
r = (ycap - y0) - m9 * xcap;
if( r <= 0.0 )
p = (ycap - y0) / xcap;
else
{
s = (ycap - y2) - m2 * (xcap - x2);
if( s <= 0.0 ) // Lslope(xcap,ycap,-d/2,fpu+a*gmv,gmv*gmv/2)
p = (fpu - fmu) + sqrt(d * (2.0 * r + d * xcap * xcap));
else
{
t = (ycap - y1) - m1 * (xcap - x1);
if( t <= 0.0 )
{ // Lslope(xcap,ycap,-1/2.,fpu,0)
u = 2.0 * s - 2.0 * xcap * x2 + x2 * x2;
p = fpu + u / (xcap + sqrt(u + xcap * xcap));
}
else
{ // Lslope(xcap,ycap,-d/2,fpu+a*gpv,gpv*gpv/2)
u = fmu + a * gmv + d * gpv / a;
p = ((fpu - fmu) + a * (gpv - gmv))
+ sqrt(2.0 * d * t + u * u);
}
}
}
} // Rslope(xcap,ycap,-d/2,fmu+a*gmv,gmv*gmv/2)
q = -sqrt(2.0 * d * height);
}
*xc = xcap;
*p0 = p;
*q0 = q;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Posneg2B
//
// Purpose: Set bi-exponential sampler for Posneg2 in case a > 0, g <= v
//
// | exp - x^2 / 2 + f x - u x + (g - v + a x)^2 / 2
// x2 = (v-g)/a |
// | exp - x^2 / 2 + f x - u x
// x = 0 |
// | exp - x^2 / 2 + f x + u x
// x1 = -(v+g)/a |
// | exp - x^2 / 2 + f x + u x + (g + v + a x)^2 / 2
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Posneg2B(
double fpu, // I f + u > 0
double fmu, // I f - u <= gmv
double gpv, // I g + v > 0
double gmv, // I g - v <= 0
double a, // I > 0
double height, // I height of sampler cusp above target maximum
double* xc, // O target maximum
double* p0, // O sampler exp( yc + p0 * (x-xc) ), x < xc (p0 > 0)
double* q0) // O sampler exp( yc + q0 * (x-xc) ), x > xc (q0 < 0)
{
double x1, x2, y1, y2, m1, m2;
double d, p, q, t;
d = 1.0 - a * a;
if( d < DBL_EPSILON )
d = DBL_EPSILON;
// Set boundaries ...(x1,y1)...(0,y0)...(x2,y2)...
// and slopes m1 +ve,-ve m2
x1 = -gpv / a;
y1 = x1 * (fpu - 0.5 * x1);
m1 = fpu - x1;
x2 = -gmv / a;
y2 = x2 * (fmu - 0.5 * x2);
m2 = fmu - x2;
// Apex necessarily at 0
t = m1 * x1 - (y1 - height);
if( t <= 0.0 ) // Lslope(0,height,-1/2.,fpu,0)
p = fpu + sqrt(2.0 * height);
else // Lslope(0,height,-d/2,fpu+a*gpv,gpv*gpv/2)
p = (fpu + a * gpv) + sqrt(d * (2.0 * t + d * x1 * x1));
t = m2 * x2 - (y2 - height);
if( t <= 0.0 ) // Rslope(0,height,-1/2.,fmu,0)
q = fmu - sqrt(2.0 * height);
else // Rslope(0,height,-d/2,fmu+a*gmv,gmv*gmv/2)
q = (fmu + a * gmv) - sqrt(d * (2.0 * t + d * x2 * x2));
*xc = 0.0;
*p0 = p;
*q0 = q;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Posneg2C
//
// Purpose: Set bi-exponential sampler for Posneg2 in case a < 0, g >= v
//
// | exp - x^2 / 2 + f x - u x + (g + v + a x)^2 / 2
// x2 = (g+v)/-a |
// | exp - x^2 / 2 + f x - u x
// x1 = (g-v)/-a |
// | exp - x^2 / 2 + f x - u x + (g - v + a x)^2 / 2
// x = 0 |
// | exp - x^2 / 2 + f x + u x + (g - v + a x)^2 / 2
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Posneg2C(
double fpu, // I f + u > 0
double fmu, // I f - u <= gmv
double gpv, // I g + v > 0
double gmv, // I g - v >= 0
double a, // I < 0
double height, // I height of sampler cusp above target maximum
double* xc, // O target maximum
double* p0, // O sampler exp( yc + p0 * (x-xc) ), x < xc (p0 > 0)
double* q0) // O sampler exp( yc + q0 * (x-xc) ), x > xc (q0 < 0)
{
double x1, x2, y0, y1, y2, m0, m1, m2, m9, xcap, ycap;
double d, p, q, r, s, t, u, v, w;
d = 1.0 - a * a;
if( d < DBL_EPSILON )
d = DBL_EPSILON;
// Set boundaries ...(0,y0)...(x1,y1)...(x2,y2)...
// and slopes m9,m0 m1 m2
x1 = -gmv / a;
y1 = x1 * (fmu - 0.5 * x1);
m1 = fmu - x1; // -ve
x2 = -gpv / a;
y2 = x2 * (fmu - 0.5 * x2);
m2 = fmu - x2;
y0 = 0.5 * gmv * gmv;
m9 = fpu + a * gmv;
m0 = fmu + a * gmv;
// Find apex from slope=0
if( m9 < 0.0 ) // apex in (-inf,0)
{
xcap = m9 / d;
ycap = height + 0.5 * (d * xcap * xcap + gmv * gmv);
p = sqrt(2. * d * height);// Lslope(xcap,ycap,-d/2,fpu+a*gmv,gmv*gmv/2)
r = (ycap - y0) - m9 * xcap;
if( r <= 0.0 ) // Rslope(xcap,ycap,-d/2,fpu+a*gmv,gmv*gmv/2)
q = -sqrt(2.0 * d * height);
else
{
s = ycap - y0 - m0 * xcap;
if( s <= 0.0 )
q = (ycap - y0) / xcap;
else
{
t = (ycap - y1) - m1 * (xcap - x1);
if( t <= 0.0 ) // Rslope(xcap,ycap,-d/2,fmu+a*gmv,gmv*gmv/2)
q = (fmu - fpu) - sqrt(2.0 * d * s + m9 * m9);
else
{
u = (ycap - y2) - m2 * (xcap - x2);
if( u <= 0.0 )
{ // Rslope(xcap,ycap,-1/2.,fmu,0)
v = 2.0 * t + (x1 - xcap) * (x1 - xcap);
w = (x1 * (1.0 - a) - 2.0 * xcap) * x1 * (1.0 + a);
q = (fmu - gmv)
- (2.0 * t + w) / (gmv - xcap + sqrt(v));
}
else
{ // Rslope(xcap,ycap,-d/2,fmu+a*gpv,gpv*gpv/2)
v = d * (2.0 * u + d * (x2 - xcap) * (x2 - xcap));
w = (fmu - fpu) - a * (gmv - gpv);
q = w - sqrt(v);
}
}
}
}
}
else if( m0 <= 0.0 ) // apex at 0
{
xcap = 0.0;
ycap = height + 0.5 * gmv * gmv;
// Lslope(xcap,ycap,-d/2,fpu+a*gmv,gmv*gmv/2)
p = m9 + sqrt(2.0 * d * height);
r = ycap - y1 + m1 * x1;
if( r <= 0.0 ) // Rslope(xcap,ycap,-d/2,fmu+a*gmv,gmv*gmv/2)
q = m0 - sqrt(2.0 * d * height);
else
{
s = ycap - y2 + m2 * x2;
if( s <= 0.0 ) // Rslope(xcap,ycap,-1/2.,fmu,0)
q = (fmu - gmv) - 2.0 * height / (gmv + sqrt(2.0 * ycap));
else // Rslope(xcap,ycap,-d/2,fmu+a*gpv,gpv*gpv/2)
q = (a * (gpv - gmv) + m0) - sqrt(d * (2.0 * s + d * x2 * x2));
}
}
else // apex in (0,x1)
{
xcap = m0 / d;
ycap = height + 0.5 * (d * xcap * xcap + gmv * gmv);
r = ycap - y0 - m0 * xcap;
if( r <= 0.0 ) // Lslope(xcap,ycap,-d/2,fmu+a*gmv,gmv*gmv/2)
p = sqrt(2.0 * d * height);
else
{
s = ycap - y0 - m9 * xcap;
if( s <= 0.0 )
p = (ycap - y0) / xcap;
else // Lslope(xcap,ycap,-d/2,fpu+a*gmv,gmv*gmv/2)
p = (fpu - fmu) + sqrt(d * (2.0 * s + d * xcap * xcap));
}
r = ycap - y1 - m1 * (xcap - x1);
if( r <= 0.0 ) // Rslope(xcap,ycap,-d/2,fmu+a*gmv,gmv*gmv/2)
q = -sqrt(2.0 * d * height);
else
{
s = ycap - y2 - m2 * (xcap - x2);
if( s <= 0.0 )
{ // Rslope(xcap,ycap,-1/2.,fmu,0)
v = 2.0 * r + (fmu - x1) * (fmu - x1) * (1.0 + a * a) / d;
w = (a * (gmv - fmu) - (1.0 + a) * gmv) * a / d;
q = -v / (w + sqrt(2.0 * r + (xcap - x1) * (xcap - x1)));
}
else // Rslope(xcap,ycap,-d/2,fmu+a*gpv,gpv*gpv/2)
q = a * (gpv - gmv)
- sqrt(d * (2.0 * s + d * (xcap - x2) * (xcap - x2)));
}
}
*xc = xcap;
*p0 = p;
*q0 = q;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Posneg2D
//
// Purpose: Set bi-exponential sampler for Posneg2 in case a < 0, g <= v
//
// | exp - x^2 / 2 + f x - u x + (g + v + a x)^2 / 2
// x2 = (v+g)/-a |
// | exp - x^2 / 2 + f x - u x
// x = 0 |
// | exp - x^2 / 2 + f x + u x
// x1 = -(v-g)/-a |
// | exp - x^2 / 2 + f x + u x + (g - v + a x)^2 / 2
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Posneg2D(
double fpu, // I f + u > 0
double fmu, // I f - u <= gmv
double gpv, // I g + v > 0
double gmv, // I g - v <= 0
double a, // I < 0
double height, // I height of sampler cusp above target maximum
double* xc, // O target maximum
double* p0, // O sampler exp( yc + p0 * (x-xc) ), x < xc (p0 > 0)
double* q0) // O sampler exp( yc + q0 * (x-xc) ), x > xc (q0 < 0)
{
double x1, x2, y1, y2, m1, m2;
double d, p, q, t;
d = 1.0 - a * a;
if( d < DBL_EPSILON )
d = DBL_EPSILON;
// Set boundaries ...(x1,y1)...(0,y0)...(x2,y2)...
// and slopes m1 +ve,-ve m2
x1 = -gmv / a;
y1 = x1 * (fpu - 0.5 * x1);
m1 = fpu -⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -