Re: [FFmpeg-devel] [PATCH] avutil/eval: Use even better PRNG

[Michael Niedermayer <michael@niedermayer.cc>]

[-- Attachment #1.1: Type: text/plain, Size: 8288 bytes --]

On Sun, Jan 14, 2024 at 03:14:23PM +0100, Stefano Sabatini wrote:
> On date Saturday 2024-01-13 04:51:06 +0100, Michael Niedermayer wrote:
> > This is the 64bit version of Chris Doty-Humphreys SFC64
> > 
> > Compared to the LCGs these produce much better quality numbers.
> > Compared to LFGs this needs less state. (our LFG has 224 byte
> > state for its 32bit version) this has 32byte state
> > Also the initialization for our LFG is slower.
> > This is also much faster than KISS or PCG.
> > 
> > This commit replaces the broken LCG used before.
> > (broken as it had only a period ~200M due to being put in a double)
> > 
> > This changes the output from random() which is why libswresample.mak
> > is updated, update was done using the command in libswresample.mak
> > 
> > Signed-off-by: Michael Niedermayer <michael@niedermayer.cc>
> > ---
> >  libavutil/eval.c             |  24 +++-
> >  libavutil/sfc64.h            |  85 ++++++++++++++
> >  tests/fate/libswresample.mak | 208 +++++++++++++++++------------------
> >  tests/ref/fate/eval          |   2 +-
> >  4 files changed, 210 insertions(+), 109 deletions(-)
> >  create mode 100644 libavutil/sfc64.h
> > 
> > diff --git a/libavutil/eval.c b/libavutil/eval.c
> > index dc6b3697bc2..349015d4fa3 100644
> > --- a/libavutil/eval.c
> > +++ b/libavutil/eval.c
> > @@ -35,6 +35,7 @@
> >  #include "internal.h"
> >  #include "log.h"
> >  #include "mathematics.h"
> > +#include "sfc64.h"
> >  #include "time.h"
> >  #include "avstring.h"
> >  #include "timer.h"
> > @@ -55,6 +56,7 @@ typedef struct Parser {
> >      void *log_ctx;
> >  #define VARS 10
> >      double *var;
> > +    FFSFC64 *prng_state;
> >  } Parser;
> >  
> >  static const AVClass eval_class = {
> > @@ -173,6 +175,7 @@ struct AVExpr {
> >      } a;
> >      struct AVExpr *param[3];
> >      double *var;
> > +    FFSFC64 *prng_state;
> >  };
> >  
> >  static double etime(double v)
> > @@ -231,8 +234,14 @@ static double eval_expr(Parser *p, AVExpr *e)
> >  
> >  #define COMPUTE_NEXT_RANDOM()                                        \
> >              int idx = av_clip(eval_expr(p, e->param[0]), 0, VARS-1); \
> > -            uint64_t r = isnan(p->var[idx]) ? 0 : p->var[idx];       \
> > -            r = r * 1664525 + 1013904223;                            \
> > +            FFSFC64 *s = p->prng_state + idx;                        \
> > +            uint64_t r;                                              \
> > +                                                                     \
> > +            if (!s->counter) {                                       \
> > +                r = isnan(p->var[idx]) ? 0 : p->var[idx];            \
> > +                ff_sfc64_init(s, r, r, r, 12);                       \
> > +            }                                                        \
> > +            r = ff_sfc64_get(s);                                     \
> >              p->var[idx] = r;                                         \
> >  
> >          case e_random: {
> > @@ -329,7 +338,11 @@ static double eval_expr(Parser *p, AVExpr *e)
> >                  case e_div: return e->value * (d2 ? (d / d2) : d * INFINITY);
> >                  case e_add: return e->value * (d + d2);
> >                  case e_last:return e->value * d2;
> > -                case e_st : return e->value * (p->var[av_clip(d, 0, VARS-1)]= d2);
> > +                case e_st :  {
> > +                    int index = av_clip(d, 0, VARS-1);
> > +                    p->prng_state[index].counter = 0;
> > +                    return e->value * (p->var[index]= d2);
> > +                }
> >                  case e_hypot:return e->value * hypot(d, d2);
> >                  case e_atan2:return e->value * atan2(d, d2);
> >                  case e_bitand: return isnan(d) || isnan(d2) ? NAN : e->value * ((long int)d & (long int)d2);
> > @@ -349,6 +362,7 @@ void av_expr_free(AVExpr *e)
> >      av_expr_free(e->param[1]);
> >      av_expr_free(e->param[2]);
> >      av_freep(&e->var);
> > +    av_freep(&e->prng_state);
> >      av_freep(&e);
> >  }
> >  
> > @@ -736,7 +750,8 @@ int av_expr_parse(AVExpr **expr, const char *s,
> >          goto end;
> >      }
> >      e->var= av_mallocz(sizeof(double) *VARS);
> > -    if (!e->var) {
> > +    e->prng_state = av_mallocz(sizeof(*e->prng_state) *VARS);
> > +    if (!e->var || !e->prng_state) {
> >          ret = AVERROR(ENOMEM);
> >          goto end;
> >      }
> > @@ -778,6 +793,7 @@ double av_expr_eval(AVExpr *e, const double *const_values, void *opaque)
> >  {
> >      Parser p = { 0 };
> >      p.var= e->var;
> > +    p.prng_state= e->prng_state;
> >  
> >      p.const_values = const_values;
> >      p.opaque     = opaque;
> > diff --git a/libavutil/sfc64.h b/libavutil/sfc64.h
> > new file mode 100644
> > index 00000000000..05f1e84cc68
> > --- /dev/null
> > +++ b/libavutil/sfc64.h
> > @@ -0,0 +1,85 @@
> > +/*
> > + * Copyright (c) 2024 Michael Niedermayer <michael-ffmpeg@niedermayer.cc>
> > + *
> > + * This file is part of FFmpeg.
> > + *
> > + * FFmpeg is free software; you can redistribute it and/or
> > + * modify it under the terms of the GNU Lesser General Public
> > + * License as published by the Free Software Foundation; either
> > + * version 2.1 of the License, or (at your option) any later version.
> > + *
> > + * FFmpeg is distributed in the hope that it will be useful,
> > + * but WITHOUT ANY WARRANTY; without even the implied warranty of
> > + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> > + * Lesser General Public License for more details.
> > + *
> > + * You should have received a copy of the GNU Lesser General Public
> > + * License along with FFmpeg; if not, write to the Free Software
> > + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
> > + *
> 
> > + * This is a implementation of SFC64, a 64-bit PRNG by Chris Doty-Humphrey.
> > + *
> > + * This Generator is much faster (0m1.872s) than 64bit KISS (0m3.823s) and PCG-XSH-RR-64/32 (0m2.700s)
> > + * And passes testu01 and practrand test suits.
> > + *
> 
> nit: possibly better to put this in @file

will move


> 
> > + */
> > +
> > +/**
> > + * @file
> > + * simple Pseudo Random Number Generator
> > + *
> > + */
> > +
> > +#ifndef AVUTIL_SFC64_H
> > +#define AVUTIL_SFC64_H
> > +
> > +#include <inttypes.h>
> > +
> > +typedef struct FFSFC64 {
> > +    uint64_t a,b,c,counter;
> > +} FFSFC64;
> > +
> > +static inline uint64_t ff_sfc64_get(FFSFC64 *s) {
> > +    uint64_t tmp = s->a + s->b + s->counter++;
> > +    s->a = s->b ^ (s->b >> 11);
> > +    s->b = s->c + (s->c << 3); // This is a multiply by 9
> > +    s->c = (s->c << 24 | s->c >> 40) + tmp;
> > +    return tmp;
> > +}
> > +
> 
> > +/**
> > + * returns the previous random value, and steps the generator backward.
> > + *
> > + * Its safe to take values before the first, but such values can be highly
> > + * correlated to the seeds.
> 
> Return ..., and step the generator...
> 
> It is safe to take values before the first, but such values can be highly
> correlated to the seeds.

will change

> 
> > + */
> > +static inline uint64_t ff_sfc64_reverse_get(FFSFC64 *s) {
> > +    uint64_t prev_c = s->b * 0x8E38E38E38E38E39;
> > +    uint64_t tmp = s->c - (prev_c << 24 | prev_c >> 40);
> > +    s->b = s->a ^ (s->a >> 11);
> > +    s->b ^= s->b >> 22;
> > +    s->b ^= s->b >> 44;
> > +
> > +    s->a = tmp - s->b - --s->counter;
> > +    s->c = prev_c;
> > +
> > +    return tmp;
> > +}
> > +
> > +/**
> > + * Initialize sfc64 with up to 3 seeds.
> > + *
> 
> > + * @param rounds number of rounds mixing up state during init. generally 8-18, larger numbers will help with bad quality seeds.
> > + *               12 is a good choice if all 3 seeds are equal
> 
> uppercase Generally after the dot.
> 
> Looks good to me apart for the minor nits.

ok, will apply with these changes

thx

[...]
-- 
Michael     GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB

Everything should be made as simple as possible, but not simpler.
-- Albert Einstein

[-- Attachment #1.2: signature.asc --]
[-- Type: application/pgp-signature, Size: 195 bytes --]

[-- Attachment #2: Type: text/plain, Size: 251 bytes --]

_______________________________________________
ffmpeg-devel mailing list
ffmpeg-devel@ffmpeg.org
https://ffmpeg.org/mailman/listinfo/ffmpeg-devel

To unsubscribe, visit link above, or email
ffmpeg-devel-request@ffmpeg.org with subject "unsubscribe".