Skip to content

Add JumpAhead / Advance #4

@JordanWeberAU

Description

@JordanWeberAU

Just leaving this here, in case anyone wants to use the jump ahead/rewind or distance functionality:
_increment is also known as "alternate stream"

//Jump Ahead: https://github.com/imneme/pcg-c/blob/master/src/pcg-advance-64.c
        /// <summary>
        /// If Jump(4), then it would fast simulate calling NextUInt32() 4 times
        /// </summary>
        /// <param name="jumpCount"></param>
        public void Jump(UInt64 jumpCount)
        {
            _state = JumpAhead(_state, jumpCount, Multiplier, _increment);
        }

        /// <summary>
        /// For small numbers, this takes the 'long way around'
        /// Rough estimates showing 100x slower than Jump forward by small numbers
        /// </summary>
        public void Rewind(UInt64 rewindCount)
        {
            _state = JumpAhead(_state , UInt64.MaxValue - rewindCount, Multiplier, _increment);
        }

        public void SkipOne()
        {
            _state = _state * Multiplier + _increment;
        }

        /* Multi-step advance functions (jump-ahead, jump-back)
         *
         * The method used here is based on Brown, "Random Number Generation
         * with Arbitrary Stride,", Transactions of the American Nuclear
         * Society (Nov. 1994).  The algorithm is very similar to fast
         * exponentiation.
         *
         * Even though delta is an unsigned integer, we can pass a
         * signed integer to go backwards, it just goes "the long way round".
         * 
         * Fast because of delta /= 2; cutting jump calculations in half each iteration
         */
        UInt64 JumpAhead(UInt64 state, UInt64 delta, UInt64 cur_mult,
                                    UInt64 cur_plus)
        {
            UInt64 acc_mult = 1;
            UInt64 acc_plus = 0;
            while (delta > 0)
            {
                //$"d{delta} - m{acc_mult} - p{acc_plus} - cm{cur_mult} - cp{cur_plus}";
                if (delta % 2 == 1) //0 or 1,  and only run if 1, because of delta /=2
                {
                    acc_mult *= cur_mult;
                    acc_plus = acc_plus * cur_mult + cur_plus;
                }
                cur_plus = (cur_mult + 1) * cur_plus;
                cur_mult *= cur_mult;
                delta /= 2;
            }
            return acc_mult * state + acc_plus;
        }

        public UInt64 Distance(PcgRng originalState)
        {
            if (originalState._increment != _increment) // || originalState.Multiplier != Multiplier)
                throw new Exception("Distance can only be calculated on the same streams");
            return Distance(originalState._state , _state , Multiplier, _increment);
        }
        public static UInt64 operator -(PcgRng lhs, PcgRng rhs) => lhs.Distance(rhs);

        UInt64 Distance(UInt64 cur_state, UInt64 newstate, UInt64 cur_mult,
                                    UInt64 cur_plus, UInt64 mask = UInt64.MaxValue)
        {
            const bool is_mcg = false;
            UInt64 the_bit = is_mcg ? 4ul : 1ul;
            UInt64 distance = 0ul;
            while ((cur_state & mask) != (newstate & mask))
            {
                if ((cur_state & the_bit) != (newstate & the_bit))
                {
                    cur_state = cur_state * cur_mult + cur_plus;
                    distance |= the_bit;
                }
                //Assert(cur_state % the_bit == 1) == (newstate % the_bit == 1)
                the_bit <<= 1;
                cur_plus = (cur_mult + 1) * cur_plus;
                cur_mult *= cur_mult;
            }
            return is_mcg ? distance >> 2 : distance;
        }

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions