On date Wednesday 2024-01-03 23:24:33 +0100, Michael Niedermayer wrote:
> This is the 64bit version of George Marsaglias KISS PRNG
>
> 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 for the 64bit version
> Also the initialization for our LFG is slower.
> OTOH the LFG is probably faster.
>
> Signed-off-by: Michael Niedermayer <mich...@niedermayer.cc>
> ---
> libavutil/eval.c | 24 ++--
> libavutil/kiss_prng.h | 64 +++++++++++
> tests/fate/libswresample.mak | 210 +++++++++++++++++------------------
> tests/ref/fate/eval | 2 +-
> 4 files changed, 184 insertions(+), 116 deletions(-)
> create mode 100644 libavutil/kiss_prng.h
>
> diff --git a/libavutil/eval.c b/libavutil/eval.c
> index 89c61ba4bf5..20b4719a75b 100644
> --- a/libavutil/eval.c
> +++ b/libavutil/eval.c
> @@ -33,6 +33,7 @@
> #include "eval.h"
> #include "ffmath.h"
> #include "internal.h"
> +#include "kiss_prng.h"
> #include "log.h"
> #include "mathematics.h"
> #include "time.h"
> @@ -55,7 +56,7 @@ typedef struct Parser {
> void *log_ctx;
> #define VARS 10
> double *var;
> - uint64_t *var_uint64;
> + KISS64State *prng_state;
> } Parser;
>
> static const AVClass eval_class = {
> @@ -174,7 +175,7 @@ struct AVExpr {
> } a;
> struct AVExpr *param[3];
> double *var;
> - uint64_t *var_uint64;
> + KISS64State *prng_state;
> };
>
> static double etime(double v)
> @@ -232,10 +233,13 @@ static double eval_expr(Parser *p, AVExpr *e)
> }
> case e_random:{
> int idx= av_clip(eval_expr(p, e->param[0]), 0, VARS-1);
> - uint64_t r= p->var_uint64[idx] ? p->var_uint64[idx] :
> (isnan(p->var[idx]) ? 0 : p->var[idx]);
> - r= r*1664525+1013904223;
> + KISS64State *s = p->prng_state + idx;
> + uint64_t r;
> +
> + if (!s->y)
> + init_random64(s, isnan(p->var[idx]) ? 0 : p->var[idx]);
> + r = get_random64(s);
needs to be rebased after the randomi introduction (also looks like
an incremental patch on top of your previous patch)
> p->var[idx]= r;
> - p->var_uint64[idx]= r;
> return e->value * (r * (1.0/UINT64_MAX));
> }
> case e_while: {
> @@ -324,7 +328,7 @@ static double eval_expr(Parser *p, AVExpr *e)
> case e_last:return e->value * d2;
> case e_st : {
> int index = av_clip(d, 0, VARS-1);
> - p->var_uint64[index] = 0;
> + p->prng_state[index].y = 0;
> return e->value * (p->var[index]= d2);
> }
> case e_hypot:return e->value * hypot(d, d2);
> @@ -346,7 +350,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->var_uint64);
> + av_freep(&e->prng_state);
> av_freep(&e);
> }
>
> @@ -732,8 +736,8 @@ int av_expr_parse(AVExpr **expr, const char *s,
> goto end;
> }
> e->var= av_mallocz(sizeof(double) *VARS);
> - e->var_uint64= av_mallocz(sizeof(uint64_t) *VARS);
> - if (!e->var || !e->var_uint64) {
> + e->prng_state = av_mallocz(sizeof(*e->prng_state) *VARS);
> + if (!e->var || !e->prng_state) {
> ret = AVERROR(ENOMEM);
> goto end;
> }
> @@ -775,7 +779,7 @@ double av_expr_eval(AVExpr *e, const double
> *const_values, void *opaque)
> {
> Parser p = { 0 };
> p.var= e->var;
> - p.var_uint64= e->var_uint64;
> + p.prng_state= e->prng_state;
>
> p.const_values = const_values;
> p.opaque = opaque;
> diff --git a/libavutil/kiss_prng.h b/libavutil/kiss_prng.h
> new file mode 100644
> index 00000000000..fb3941b5df5
> --- /dev/null
> +++ b/libavutil/kiss_prng.h
> @@ -0,0 +1,64 @@
> +/*
> + * Copyright (c) 2008,2024 Michael Niedermayer
> <michael-ffm...@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
> + *
> + * 64-bit PRNG by George Marsaglia
> https://www.thecodingforums.com/threads/64-bit-kiss-rngs.673657/
> + *
> + */
> +
> +/**
> + * @file
> + * simple Pseudo Random Number Generator
> + *
> + * The implementation was verified to match the value after 100 milllion
> iterations
> + */
> +
> +#ifndef AVUTIL_KISS_PRNG_H
> +#define AVUTIL_KISS_PRNG_H
> +
> +#include <inttypes.h>
> +
> +typedef struct KISS64State{
> + uint64_t x,y,z,c;
> +}KISS64State;
missing namespace here and in the functions below.
For this you could probably use AVKISS64PRNG (or AVKISS64), and use
av_kiss64prng_ or av_kiss64_ as function prefix (I'd probably go with
the shortest variant).
> +
> +static inline uint64_t get_random64(KISS64State *s){
av_kiss64_get() ?
> + uint64_t
> + t = (s->x << 58) + s->c;
> + s->c = s->x >> 6;
> + s->x+= t;
> + s->c+= s->x < t;
> +
> + s->y ^= s->y << 13;
> + s->y ^= s->y >> 17;
> + s->y ^= s->y << 43;
> +
> + s->z = 1234567 + 6906969069LL * s->z;
> +
> + return s->x + s->y + s->z;
> +}
> +
> +static inline void init_random64(KISS64State *s, uint64_t seed){
av_kiss64_init() ?
to make it consistent with the av_lfg API.
> + //Constants based on SHA512 of "FFmpeg" and Marsaglias values so seed=0
> matches his
> + s->x = 1234567890987654321ULL ^ seed;
> + s->y =( 362436362436362436ULL ^ (0xd255973df01e5086*seed)) | 4;
> + s->z = 1066149217761810ULL + seed;
> + s->c = 123456123456123456ULL - 0x33730c0524f137da*seed;
> +}
> +
> +#endif /* AVUTIL_KISS__PRNGH */
> diff --git a/tests/fate/libswresample.mak b/tests/fate/libswresample.mak
> index 0d29f760248..ea421d45fb6 100644
> --- a/tests/fate/libswresample.mak
> +++ b/tests/fate/libswresample.mak
> @@ -359,7 +359,7 @@ fate-swr-resample_nn-s16p-8000-44100: SIZE_TOLERANCE =
> 96000 - 20480
> define ARESAMPLE_ASYNC
> FATE_SWR_RESAMPLE += fate-swr-resample_async-$(3)-$(1)-$(2)
> fate-swr-resample_async-$(3)-$(1)-$(2): tests/data/asynth-$(1)-1.wav
> -fate-swr-resample_async-$(3)-$(1)-$(2): CMD = ffmpeg -i
> $(TARGET_PATH)/tests/data/asynth-$(1)-1.wav -af
> atrim=end_sample=10240,asetpts=PTS+random\(0\)*200-100,aresample=$(2):async=50:min_hard_comp=0.100000:first_pts=0:linear_interp=0:exact_rational=0:internal_sample_fmt=$(3),aformat=$(3),aresample=$(1):linear_interp=0:exact_rational=0:internal_sample_fmt=$(3)
> -f wav -c:a pcm_s16le -
> +fate-swr-resample_async-$(3)-$(1)-$(2): CMD = ffmpeg -i
> $(TARGET_PATH)/tests/data/asynth-$(1)-1.wav -af
> atrim=end_sample=10240,asetpts=PTS+random\(1\)*200-100,aresample=$(2):async=50:min_hard_comp=0.100000:first_pts=0:linear_interp=0:exact_rational=0:internal_sample_fmt=$(3),aformat=$(3),aresample=$(1):linear_interp=0:exact_rational=0:internal_sample_fmt=$(3)
> -f wav -c:a pcm_s16le -
any reason to change the random() index?
[...]
> Evaluating '7000000B*random(0)'
> -'7000000B*random(0)' -> 0.003078
> +'7000000B*random(0)' -> 27118453.055396
unrelated comment: we might provide a better example/reference here
_______________________________________________
ffmpeg-devel mailing list
ffmpeg-devel@ffmpeg.org
https://ffmpeg.org/mailman/listinfo/ffmpeg-devel
To unsubscribe, visit link above, or email
ffmpeg-devel-requ...@ffmpeg.org with subject "unsubscribe".
