anchor.diff:      end.
   
   %% =====================================================================
  +%% exrop PRNG: Xoroshiro116+
  +%%
  +%% Reference URL: http://xorshift.di.unimi.it/
  +%%
  +%% 58 bits fits into an immediate on 64bits Erlang and is thus much faster.
  +%% In fact, an immediate number is 60 bits signed in Erlang so you can
  +%% add two positive 58 bit numbers and get a 59 bit number that still is
  +%% a positive immediate, which is a property we utilize here...
  +%%
  +%% Modification of the original Xororhiro128+ algorithm to 116 bits
  +%% by Sebastiano Vigna.  A lot of thanks for his help and work.
  +%% =====================================================================
  +%% (a, b, c) = (24, 2, 35)
  +%% JUMP Polynomial = 0x9863200f83fcd4a11293241fcb12a (116 bit)
  +%%
  +%% From http://xoroshiro.di.unimi.it/xoroshiro116plus.c:
  +%% ---------------------------------------------------------------------
  +%% /* Written in 2017 by Sebastiano Vigna (vigna@acm.org).
  +%%
  +%% To the extent possible under law, the author has dedicated all copyright
  +%% and related and neighboring rights to this software to the public domain
  +%% worldwide. This software is distributed without any warranty.
  +%%
  +%% See <http://creativecommons.org/publicdomain/zero/1.0/>. */
  +%%
  +%% #include <stdint.h>
  +%%
  +%% #define UINT58MASK (uint64_t)((UINT64_C(1) << 58) - 1)
  +%%
  +%% uint64_t s[2];
  +%%
  +%% static inline uint64_t rotl58(const uint64_t x, int k) {
  +%%     return (x << k) & UINT58MASK | (x >> (58 - k));
  +%% }
  +%% 
  +%% uint64_t next(void) {
  +%%     uint64_t s1 = s[1];
  +%%     const uint64_t s0 = s[0];
  +%%     const uint64_t result = (s0 + s1) & UINT58MASK;
  +%%
  +%%     s1 ^= s0;
  +%%     s[0] = rotl58(s0, 24) ^ s1 ^ ((s1 << 2) & UINT58MASK); // a, b
  +%%     s[1] = rotl58(s1, 35); // c
  +%%     return result;
  +%% }
  +%%
  +%% void jump(void) {
  +%%     static const uint64_t JUMP[] =
  +%%         { 0x4a11293241fcb12a, 0x0009863200f83fcd };
  +%%
  +%%     uint64_t s0 = 0;
  +%%     uint64_t s1 = 0;
  +%%     for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
  +%%         for(int b = 0; b < 64; b++) {
  +%%             if (JUMP[i] & UINT64_C(1) << b) {
  +%%                 s0 ^= s[0];
  +%% 	           s1 ^= s[1];
  +%% 	       }
  +%% 	       next();
  +%% 	   }
  +%%     s[0] = s0;
  +%%     s[1] = s1;
  +%% }
  +
  +-opaque exrop_state() :: nonempty_improper_list(uint58(), uint58()).
  +
  +-dialyzer({no_improper_lists, exrop_seed/1}).
  +exrop_seed({A1, A2, A3}) ->
  +    [_|S1] =
  +        exrop_next_s(
  +          ?MASK(58, (A1 * 4294967197) + 1),
  +          ?MASK(58, (A2 * 4294967231) + 1)),
  +    exrop_next_s(?MASK(58, (A3 * 4294967279) + 1), S1).
  +
  +-dialyzer({no_improper_lists, exrop_next_s/2}).
  +%% Advance xoroshiro116+ state one step
  +%% [a, b, c] = [24, 2, 35]
  +-define(
  +   exrop_next_s(S0, S1, S1_a),
  +   begin
  +       S1_a = S1 bxor S0,
  +       [?ROTL(58, S0, 24) bxor S1_a bxor ?BSL(58, S1_a, 2)| % a, b
  +        ?ROTL(58, S1_a, 35)] % c
  +   end).
  +exrop_next_s(S0, S1) ->
  +    ?exrop_next_s(S0, S1, S1_a).
  +
  +-dialyzer({no_improper_lists, exrop_next/1}).
  +%% Advance xoroshiro116+ state one step, generate 58 bit unsigned integer,
  +%% and waste the lowest bit since it is of lower randomness quality
  +exrop_next([S0|S1]) ->
  +    {?MASK(58, S0 + S1), ?exrop_next_s(S0, S1, S1_a)}.
  +
  +exrop_uniform({Alg, R}) ->
  +    {V, R1} = exrop_next(R),
  +    %% Waste the lowest bit since it is of lower
  +    %% randomness quality than the others
  +    {(V bsr (58-53)) * ?TWO_POW_MINUS53, {Alg, R1}}.
  +
  +exrop_uniform(Range, {Alg, R}) ->
  +    {V, R1} = exrop_next(R),
  +    MaxMinusRange = ?BIT(58) - Range,
  +    ?uniform_range(Range, Alg, R1, V, MaxMinusRange, I).
  +
  +%% Split a 116 bit constant into two 58 bit words,
  +%% a top '1' marks the end of the low word.
  +-define(
  +   JUMP_116(Jump),
  +   [?BIT(58) bor ?MASK(58, (Jump)),(Jump) bsr 58]).
  +%%
  +exrop_jump({Alg,S}) ->
  +    [J|Js] = ?JUMP_116(16#9863200f83fcd4a11293241fcb12a),
  +    {Alg, exrop_jump(S, 0, 0, J, Js)}.
  +%%
  +-dialyzer({no_improper_lists, exrop_jump/5}).
  +exrop_jump(_S, S0, S1, 0, []) -> % End of jump constant
  +    [S0|S1];
  +exrop_jump(S, S0, S1, 1, [J|Js]) -> % End of word
  +    exrop_jump(S, S0, S1, J, Js);
  +exrop_jump([S__0|S__1] = _S, S0, S1, J, Js) ->
  +    case ?MASK(1, J) of
  +        1 ->
  +            NewS = exrop_next_s(S__0, S__1),
  +            exrop_jump(NewS, S0 bxor S__0, S1 bxor S__1, J bsr 1, Js);
  +        0 ->
  +            NewS = exrop_next_s(S__0, S__1),
  +            exrop_jump(NewS, S0, S1, J bsr 1, Js)
  +    end.
  +
  +%% =====================================================================
   %% Ziggurat cont
   %% =====================================================================
   -define(NOR_R, 3.6541528853610087963519472518).
   -define(NOR_INV_R, 1/?NOR_R).
   
   %% return a {sign, Random51bits, State}
  +get_52({Alg=#{bits:=Bits, next:=Next}, S0}) ->
  +    %% Use the high bits
  +    {Int,S1} = Next(S0),
  +    {?BIT(Bits - 51 - 1) band Int, Int bsr (Bits - 51), {Alg, S1}};
   get_52({Alg=#{next:=Next}, S0}) ->
       {Int,S1} = Next(S0),
  -    {((1 bsl 51) band Int), Int band ((1 bsl 51)-1), {Alg, S1}}.
  +    {?BIT(51) band Int, ?MASK(51, Int), {Alg, S1}}.
   
   %% Slow path
   normal_s(0, Sign, X0, State0) ->
  
positive.code: %%
  %% %CopyrightBegin%
  %%
  %% Copyright Ericsson AB 2015-2016. All Rights Reserved.
  %%
  %% Licensed under the Apache License, Version 2.0 (the "License");
  %% you may not use this file except in compliance with the License.
  
negative.code: %% Modification of the original Xorshift128+ algorithm to 116
  %% by Sebastiano Vigna, a lot of thanks for his help and work.
  %% =====================================================================
  -opaque exsplus_state() :: nonempty_improper_list(uint58(), uint58()).
  
  -dialyzer({no_improper_lists, exsplus_seed/1}).
  
  exsplus_seed({A1, A2, A3}) ->
      {_, R1} = exsplus_next(
                  [?MASK(58, (A1 * 4294967197) + 1)|
                   ?MASK(58, (A2 * 4294967231) + 1)]),
      {_, R2} = exsplus_next(
                  [?MASK(58, (A3 * 4294967279) + 1)|
                   tl(R1)]),
      R2.
  
  -dialyzer({no_improper_lists, exsplus_next/1}).
  