[FFmpeg-devel] [PATCH] Box-Muller gaussian generator

[Stefano Sabatini]

On date Friday 2009-12-11 01:49:21 +0100, Michael Niedermayer encoded:
> On Fri, Dec 11, 2009 at 01:40:08AM +0100, Michael Niedermayer wrote:
> > On Sun, Dec 06, 2009 at 11:13:37PM +0100, Stefano Sabatini wrote:
> > > On date Sunday 2009-12-06 19:40:20 +0100, Vitor Sessak encoded:
> > > > Stefano Sabatini wrote:
> > > >> Hi,
> > > >>
> > > >> I'm going to use this for a noise filter, I believe this could be
> > > >> useful in other part of libav* though.
> > > >
> > > >> +double av_bmg_get(AVBMG *bmg)
> > > >> +{
> > > >> +    double x1, x2, w, y;
> > > >> +    AVLFG *lfg = &bmg->lfg;
> > > >> +
> > > >> +    do {
> > > >> +        x1 = 2.0 * ranf(lfg) - 1.0;
> > > >> +        x2 = 2.0 * ranf(lfg) - 1.0;
> > > >> +        w = x1*x1 + x2*x2;
> > > >> +    } while (w >= 1.0);
> > > >> +
> > > >> +    w = sqrt((-2.0 * log(w)) / w);
> > > >> +    y = x2 * w;
> > > >> +
> > > >> +    return bmg->mean + y * bmg->std_dev;
> > > >> +}
> > > >
> > > > This generates two numbers to use only one. Isn't it faster to store the  
> > > > second number for the next call or have a function to fill a buffer of  
> > > > random numbers?
> > > 
> > > Implemented the first idea.
> > 
> > void av_bmg_get(AVLFG *lfg, double out[2]);
> >     double x1, x2, w;
> > 
> >     do {
> >         x1 = 2.0/UINT_MAX*av_lfg_get(lfg) - 1.0;
> >         x2 = 2.0/UINT_MAX*av_lfg_get(lfg) - 1.0;
> >         w = x1*x1 + x2*x2;
> >     } while (w >= 1.0);
> >
> >     w = sqrt((-2.0 * log(w)) / w);

If av_lfg_get() returns UINT_MAX/2 two times in a row then we have:
x1 = x2 = w = 0
while this isn't very likely (I don't even know if av_lfg() can return
two times in a row the same number) we'll have a division-by-zero.
Is this the bug you mentioned in your first mail?
Also how it is fixed in this second variant?

> >     out[0] = x1 * w;
> >     out[1] = x2 * w;
> > }
> > 
> > comments welcome, if none ill commit above
> > 
> > libavutil is supposed to be clean and i wont tolerate it being filled up with
> 
> hmm, the word tolerate sounds a little harsh, i guess i could have said this
> more diplomatically but the point stands we dont need a new file and whole
> struct+API for something that could be a single 8 line function

I appreciate your wish to keep simple the API, well this slightly
complicates its usage which becomes:

int main(void)
{
    int i;
    AVLFG lfg;
    double mean   = 1000;
    double stddev = 53;

    av_lfg_init(&lfg, 42);

    for (i = 0; i < 1000; i+=2) {
        double bmg_out[2];
        av_bmg_get(&lfg, bmg_out);
        printf("%f\n%f\n",
               bmg_out[0] * stddev + mean,
               bmg_out[1] * stddev + mean);
    }

    return 0;
}

but I can easily live with it and maybe I even prefer your variant.

Regards.
-- 
FFmpeg = Fundamental and Forgiving Mastering Portable EnGine