yummy
2023-01-11 10:40:31
昨天晚上半夜发现自己写的东西假掉了，今天早上赶紧起来更新。
然后我发现昨天那玩意可以扩展到全体多项式，虽然求解
最近翻了翻自己非常著名的 Kotori 的数据生成器，发现当时实际上已经利用了一个不均匀随机数生成器，但是当时对这个随机数生成器并没有仔细研究。正好现在是放假，于是我想花一点时间来仔细思考下这个问题，但是后来发现事情并不是这么简单。
因为 yummy 没有系统学习过概率论，也没系统学习过多项式，还没系统学习过线性代数，因此很多理解可能是很 trivial 的，一些记号也可能是错的。
本文将提供一种方法，
前置知识：概率密度函数。 建议大家看 这篇日报 的“概念 - 分布”这一小块。这篇日报介绍的分布都比较高级，但很多使用了逆分布法和拒绝采样法。
普通高次方程不可解，因此逆分布法失效了；同时拒绝采样法没有高超的技巧的话期望调用次数很高，最坏调用次数为正无穷。
所以我们要采用别的方式。本文分为下面几个部分：
例 1：构造随机数生成器，使
P(x)=0.5+x,x\in (0,1) 。
线性分布中，如果作出
如果我们在这个直角梯形内均匀选取一个点（此处“均匀”的意思是，任取两块面积相同的区域，选中点在这两块区域内概率相等），那么感性理解一下，这个点的横坐标就是服从
问题转化成怎么在一个直角梯形内选取一个点。
用小学知识我们知道，两个直角梯形可以拼成一个矩形。
以例 1 为例，用两个直角梯形我们拼出了
可以证明，这种选取方式是均匀的，因此返回值服从
显然地，这种方法可以推广到任意合法的概率分布
小练习 1：作出
P(x)=\dfrac{2}{\pi}\sqrt{1-x^2},x\in (-1,1) 的图像，并生成服从该分布的随机变量。
例 2：
a,b,c 是三个(0,1) 内均匀独立的随机变量，分别求它们最小值、最大值、中位数的概率密度函数。
我们可以求累积分布函数
先看最大值。
例 6：构造随机数生成器，使
P(x)=\dfrac{3}{19}(2x-1)^2(x^2+x+1) 。
接下来将解决重根和虚根的问题。
首先有一个经典结论：实系数多项式的虚根总是两两共轭。
不难发现，
换言之，一个多项式总能被分解成若干实系数一次式和若干无实根实系数二次式的积，只不过高次的咱不会。
回顾例 4，我们可以以对称轴为界，将
在例 6 中，
选完落在哪一段后，对于
对于
例如例 6 中，
我大概写了一下，长这个样子：
#include<bits/stdc++.h>
using namespace std;
const double Eps=1e-6;
mt19937 rnd(time(0));
uniform_real_distribution<> Get(0,1);
double genpwr(double l,double u,int x_l,int r_x)
{
vector<double> res;
for(int i=0;i<x_l+r_x+1;i++)
res.push_back(Get(rnd));
nth_element(res.begin(),res.begin()+x_l,res.end());
return res[x_l]*(u-l)+l;
}
bool cmp_real(complex<double> x,complex<double> y){return real(x)<real(y);}
bool cmp_sec(pair<double,double> x,pair<double,double> y){return x.second<y.second;}
double sqrabs(complex<double> x){return real(x)*real(x)+imag(x)*imag(x);}
vector<double> polyint(vector<double> a)
{
vector<double> res={0};
for(int i=0;i<a.size();i++)
res.emplace_back(a[i]/(i+1));
return res;
}
double polyval(vector<double> &a,double x)//有条件的可以用多项式多点求值
{
double tot=0,curx=1;
for(int i=0;i<a.size();i++)
{
tot+=a[i]*curx;
curx*=x;
}
return tot;
}
vector<double> polydivlinear(vector<double> a,double b)// makes sense when (b+x) | a
{
if(-Eps<b and b<Eps)
{
assert(-Eps<a[0] and a[0]<Eps);
a.erase(a.begin());
return a;
}
vector<double> res;
for(int i=0;i<a.size()-1;i++)
{
double topush=a[i]/b;
a[i+1]-=topush;
res.emplace_back(topush);
}
return res;
}
vector<double> polydivquad(vector<double> a,double b,double k)// makes sense when (b+kx+x^2) | a
{
if(-Eps<b and b<Eps)
{
assert(-Eps<a[0] and a[0]<Eps);
a.erase(a.begin());
return a;
}
vector<double> res;
for(int i=0;i<a.size()-2;i++)
{
double topush=a[i]/b;
a[i+1]-=topush*k;
a[i+2]-=topush;
res.emplace_back(topush);
}
return res;
}
vector<double> polymult(vector<double> a,vector<double> b)
//这点规模懒得搬我 Karatsuba 板子 or FFT 板子了
{
vector<double> res(a.size()+b.size()-1,0);
for(int i=0;i<a.size();i++)
for(int j=0;j<b.size();j++)
res[i+j]+=a[i]*b[j];
return res;
}
struct polygenerator
{
double ov_l,ov_u,Int;
vector<double> re;
vector<complex<double> > im;
vector<pair<double,double> > x_F;
vector<vector<double> > rints;
vector<vector<pair<double,double> > > iints;//前原积分，后乘上 (x-L) 积分
polygenerator(double ll,double uu,vector<double> re_rts,vector<complex<double> > im_pos_rts)
// [l,u]、定义域外实根、共轭根
{
ov_l=ll;ov_u=uu;Int=0;
re=re_rts;
im=im_pos_rts;
vector<double> f={1},F;
for(double i:re)f=polymult(f,{-i,1});
for(auto i:im)f=polymult(f,{sqrabs(i),-2*real(i),1});
F=polyint(f);
if(polyval(F,ov_u)<0)
for(int i=0;i<F.size();i++)
F[i]=-F[i];
F[0]-=polyval(F,ov_l);
x_F.emplace_back(ov_l,0);
for(int i=0;i<im.size();i++)//对于共轭根将对称轴拿出来
{
double x=im[i].real();
if(ov_l<=x and x<ov_u)
x_F.emplace_back(x,polyval(F,x));
}
sort(x_F.begin(),x_F.end());
x_F.emplace_back(ov_u,polyval(F,ov_u));
Int=x_F.back().second;
//接下来每个函数每个段求一个定积分
rints.resize(x_F.size()-1);
iints.resize(x_F.size()-1);
for(int i=0;i<re.size();i++)
{
auto tmp=polyint(polydivlinear(f,-re[i]));
for(int j=0;j+1<x_F.size();j++)
rints[j].emplace_back(abs(polyval(tmp,x_F[j+1].first)-polyval(tmp,x_F[j].first)));
}
for(int i=0;i<im.size();i++)
{
double mod=sqrabs(im[i]),P=real(im[i]);
auto fdv=polydivquad(f,mod,-2*real(im[i])),tmp=polyint(fdv);
for(int j=0;j+1<x_F.size();j++)
{
double Lrg=x_F[j].first,Rrg=x_F[j+1].first;
double Cint=abs(polyval(tmp,Rrg)-polyval(tmp,Lrg));
vector<double> tp1;
if(P<Lrg+Eps)
{
tp1=polyint(polymult(fdv,{P*P,-2*P,1}));
iints[j].emplace_back(Cint,2*(Lrg-P)*abs(polyval(tp1,Rrg)-polyval(tp1,Lrg)));
}
else
{
tp1=polyint(polymult(fdv,{P*P,-2*P,1}));
iints[j].emplace_back(Cint,2*(P-Rrg)*abs(polyval(tp1,Rrg)-polyval(tp1,Lrg)));
}
}
}
}
double operator () ()
{
double val=Get(rnd)*Int;
int rg=upper_bound(x_F.begin(),x_F.end(),make_pair(0,val),cmp_sec)-x_F.begin()-1;
double L=x_F[rg].first,R=x_F[rg+1].first,Qint=x_F[rg+1].second-x_F[rg].second;
int x_l=0,r_x=0;
for(int i=0;i<re.size();i++)
{
if(re[i]<=L)
{
double Cint=rints[rg][i]*(L-re[i]);
if(Get(rnd)*Qint>=Cint)x_l++;
}
else
{
double Cint=rints[rg][i]*(re[i]-R);
if(Get(rnd)*Qint>=Cint)r_x++;
}
}
for(int i=0;i<im.size();i++)
{
double P=real(im[i]),Q2=imag(im[i])*imag(im[i]);
if(Get(rnd)*Qint<iints[rg][i].first*Q2)continue;
double Tint=Qint-iints[rg][i].first*Q2;
if(P<=L)
{
if(Get(rnd)*Tint>=(L-P)*iints[rg][i].second)
x_l++;
if(Get(rnd)*Tint>=(L-P)*iints[rg][i].second)
x_l++;
}
else
{
if(Get(rnd)*Tint>=(P-R)*iints[rg][i].second)
r_x++;
if(Get(rnd)*Tint>=(P-R)*iints[rg][i].second)
r_x++;
}
}
return genpwr(L,R,x_l,r_x);
}
};
//polygenerator A(0,1,{1},{0.5+0.86602540378443864676372317075294i});
//polygenerator A(0,1,{},{0.5});
//polygenerator A(0,1,{-1,1},{});
//polygenerator A(0,1,{},{0.5,-0.5+0.86602540378443864676372317075294i});
void Plot(int Tms,double D,int L,int R,int perstar)
{
int res[1005]={0};
for(int i=1;i<=Tms;i++)
{
double tmp=A();
assert(L<=tmp and tmp<R);
res[int((tmp-L)*D)]++;
}
for(int i=0;i<=D;i++)
{
printf("%.3f: %6d",i/D+L,res[i]);
for(int j=0;j<res[i]/perstar;j++)
printf("*");
puts("");
}
}
int main()
{
//Plot(10,10,0,1,20);
Plot(1000000,100,0,1,400);
}小练习 1 答案：
如果忽略
不难发现，半圆弧长度
在此基础上只要随一个极角即可。
综上，在
小练习 2 答案：