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 > > --- > > 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 > > + * > > + * 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 > > + > > +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