Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Lagged Fibonacci Generator #727

Closed
wants to merge 12 commits into from

2 participants

@monarchdodra
Collaborator

D2 test results for 764

Ported from boost:
http://www.boost.org/doc/libs/1_50_0/boost/random/linear_congruential.hpp

While boost provided two different classes for UInt/Real types, this implementation merges both into the same class.

My doubts on the implementation are:
*Aliases: Boost uses doubles as the defaults for the engine, and doesn't provide any aliases for the unsigned types. Should we follow the same scheme? Provided alliases for both?
*Default Alias: Not sure which config to use for the "LaggedFibonacci", so I chose LaggedFibonacci607, but documented it as "implementation chosen".
*Reference type & unseeded engine: The generator is a reference type forward range. In thoery, it can't be used until it has been seeded. I have put in asserts about this. Maybe it should be enforece instead?
*Boost provides a "generate" function, which is used to fill a range of integers. I did not write these, as they can be replaced by simply chosing the correct engine, an manually calling std.algo.fill.

Open to all suggestions/criticism, even the "minor" ones, such as documentation/licensing.

@monarchdodra
Collaborator

I will be away from August 6 to August 20. Apologies for not being able to answer comments during that time.

std/random.d
((52 lines not shown))
+ //********
+ // Begin Payload
+ // We use a payload for "cheap" copy/move
+ // Payload implements:
+ // seed(range)
+ // fill
+ // front
+ // popFront
+ // discard
+ //********
+ struct Payload
+ {
+ Type[longLag] x;
+ size_t i = 0;
+
+ void ThrowSeedException()
@andralex Owner
andralex added a note

lowercase t

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((64 lines not shown))
+ Type[longLag] x;
+ size_t i = 0;
+
+ void ThrowSeedException()
+ {
+ throw new Exception(text("LaggedFibonacciEngine.seed: Input range didn't provide enough"
+ " elements: Need ", longLag, " elements."));
+ }
+
+ void seed(Range)(Range range)
+ {
+ static if (isUnsigned!Type)
+ {
+ i = 0;
+ typeof(range.front) front;
+ for(size_t j = 0; j < longLag; j++)
@andralex Owner
andralex added a note

use foreach

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((67 lines not shown))
+ void ThrowSeedException()
+ {
+ throw new Exception(text("LaggedFibonacciEngine.seed: Input range didn't provide enough"
+ " elements: Need ", longLag, " elements."));
+ }
+
+ void seed(Range)(Range range)
+ {
+ static if (isUnsigned!Type)
+ {
+ i = 0;
+ typeof(range.front) front;
+ for(size_t j = 0; j < longLag; j++)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
@andralex Owner
andralex added a note

one statement per line

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((68 lines not shown))
+ {
+ throw new Exception(text("LaggedFibonacciEngine.seed: Input range didn't provide enough"
+ " elements: Need ", longLag, " elements."));
+ }
+
+ void seed(Range)(Range range)
+ {
+ static if (isUnsigned!Type)
+ {
+ i = 0;
+ typeof(range.front) front;
+ for(size_t j = 0; j < longLag; j++)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
+ UIntType val = cast(UIntType)front;
@andralex Owner
andralex added a note

don't stutter, use auto

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((69 lines not shown))
+ throw new Exception(text("LaggedFibonacciEngine.seed: Input range didn't provide enough"
+ " elements: Need ", longLag, " elements."));
+ }
+
+ void seed(Range)(Range range)
+ {
+ static if (isUnsigned!Type)
+ {
+ i = 0;
+ typeof(range.front) front;
+ for(size_t j = 0; j < longLag; j++)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
+ UIntType val = cast(UIntType)front;
+ for(size_t k = 1; k < (bits+31)/32; ++k)
@andralex Owner
andralex added a note

use foreach and precompute (bits+31)/32 outside the loop

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((72 lines not shown))
+
+ void seed(Range)(Range range)
+ {
+ static if (isUnsigned!Type)
+ {
+ i = 0;
+ typeof(range.front) front;
+ for(size_t j = 0; j < longLag; j++)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
+ UIntType val = cast(UIntType)front;
+ for(size_t k = 1; k < (bits+31)/32; ++k)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
@andralex Owner
andralex added a note

ditto

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((83 lines not shown))
+ UIntType val = cast(UIntType)front;
+ for(size_t k = 1; k < (bits+31)/32; ++k)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
+ val += cast(UIntType)front << 32*k;
+ }
+ x[j] = val & mask;
+ }
+ fill();
+ }
+ else //static if (isFloatingPoint!Type)
+ {
+ i = 0;
+ typeof(range.front) front;
+ for(size_t j = 0; j < longLag; ++j)
@andralex Owner
andralex added a note

foreach

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((82 lines not shown))
+ front = range.front; range.popFront();
+ UIntType val = cast(UIntType)front;
+ for(size_t k = 1; k < (bits+31)/32; ++k)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
+ val += cast(UIntType)front << 32*k;
+ }
+ x[j] = val & mask;
+ }
+ fill();
+ }
+ else //static if (isFloatingPoint!Type)
+ {
+ i = 0;
+ typeof(range.front) front;
@andralex Owner
andralex added a note

why is this defined here? could be defined at initialization.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((101 lines not shown))
+ RealType mult = divisor;
+ for(int k = 0; k < bits/32; ++k)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
+ val += cast(uint)front * mult;
+ mult *= two32;
+ }
+ static if(mask != 0)
+ {
+ if(range.empty) ThrowSeedException();
+ front = range.front; range.popFront();
+ val += (cast(uint)front & mask) * mult;
+ }
+ assert(val >= 0);
+ assert(val < 1);
@andralex Owner
andralex added a note

consolidate asserts, print text(val) if assert fails

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((114 lines not shown))
+ }
+ assert(val >= 0);
+ assert(val < 1);
+ x[j] = val;
+ }
+ fill();
+ }
+ }
+
+ @safe nothrow
+ void fill()
+ {
+ static if (isUnsigned!Type)
+ {
+ // two loops to avoid costly modulo operations
+ for(size_t j = 0; j < shortLag; ++j)
@andralex Owner
andralex added a note

foreach

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((116 lines not shown))
+ assert(val < 1);
+ x[j] = val;
+ }
+ fill();
+ }
+ }
+
+ @safe nothrow
+ void fill()
+ {
+ static if (isUnsigned!Type)
+ {
+ // two loops to avoid costly modulo operations
+ for(size_t j = 0; j < shortLag; ++j)
+ x[j] = (x[j] + x[j+(longLag-shortLag)]) & mask;
+ for(size_t j = shortLag; j < longLag; ++j)
@andralex Owner
andralex added a note

foreach

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((122 lines not shown))
+
+ @safe nothrow
+ void fill()
+ {
+ static if (isUnsigned!Type)
+ {
+ // two loops to avoid costly modulo operations
+ for(size_t j = 0; j < shortLag; ++j)
+ x[j] = (x[j] + x[j+(longLag-shortLag)]) & mask;
+ for(size_t j = shortLag; j < longLag; ++j)
+ x[j] = (x[j] + x[j-shortLag]) & mask;
+ }
+ else //static if (isFloatingPoint!Type)
+ {
+ // two loops to avoid costly modulo operations
+ for(size_t j = 0; j < shortLag; ++j)
@andralex Owner
andralex added a note

foreach

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((129 lines not shown))
+ for(size_t j = 0; j < shortLag; ++j)
+ x[j] = (x[j] + x[j+(longLag-shortLag)]) & mask;
+ for(size_t j = shortLag; j < longLag; ++j)
+ x[j] = (x[j] + x[j-shortLag]) & mask;
+ }
+ else //static if (isFloatingPoint!Type)
+ {
+ // two loops to avoid costly modulo operations
+ for(size_t j = 0; j < shortLag; ++j)
+ {
+ RealType t = x[j] + x[j+(longLag-shortLag)];
+ if(t >= one)
+ t -= one;
+ x[j] = t;
+ }
+ for(size_t j = shortLag; j < longLag; ++j)
@andralex Owner
andralex added a note

foreach

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((140 lines not shown))
+ if(t >= one)
+ t -= one;
+ x[j] = t;
+ }
+ for(size_t j = shortLag; j < longLag; ++j)
+ {
+ RealType t = x[j] + x[j-shortLag];
+ if(t >= one)
+ t -= one;
+ x[j] = t;
+ }
+ }
+ }
+
+ @property @safe nothrow const
+ Type front()
@andralex Owner
andralex added a note

ref?

@monarchdodra Collaborator

I don't think so actually:
-None of the random generators return by ref
-Boost doesn't return by ref
-The return type is just a built in integer/floating point

Unless you mean it as a way to "write to front and modify the payload"? I don't think we should leak this kind of control, but don't have the PRNG knowledge to certify it. I'd rather not be the one to open this functionality

Maybe return by const ref then? Your call.

@andralex Owner

OK - by value's fine.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((161 lines not shown))
+ void popFront()
+ {
+ ++i;
+ if(i == longLag)
+ {
+ fill();
+ i = 0;
+ }
+ }
+
+ @safe nothrow
+ void discard(size_t n)
+ {
+ while(n >= longLag)
+ {
+ fill();
@andralex Owner
andralex added a note

two space indent?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((185 lines not shown))
+ }
+ }
+ //********
+ // End Payload
+ //********
+ Payload* payload;
+
+ public:
+
+/**
+ Creates a new $(D LaggedFibonacciEngine) and calls $(D seed()).
+*/
+ static
+ This opCall()
+ {
+ This ret;
@andralex Owner
andralex added a note

whoa, a lot of two space indent around here

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((202 lines not shown))
+ return ret;
+ }
+
+/**
+ Creates a new $(D LaggedFibonacciEngine) and calls $(D seed(value)).
+*/
+ this(T)(T value)
+ if (isIntegral!T)
+ {seed(value);}
+
+/**
+ Creates a new $(D LaggedFibonacciEngine) and calls $(D seed(range)).
+ */
+ this(Range)(Range range)
+ if (isInputRange!Range)
+ {seed(range);}
@andralex Owner
andralex added a note

odd formatting

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@andralex
Owner

ping?

@andralex
Owner

Sorry, saw the notice only now. Let's get back to this after Aug 20. Enjoy vacations!

@andralex
Owner

due for review next week

@monarchdodra
Collaborator

Thanks a lot for the re-read. In my defense, I tried to stay close to the boost code, but I think your comments are pertinent, and the code should be more D-Like.

I have nothing to add to any of the comments, except the "ref?" one.

Will provide a new version in 2 days.

@monarchdodra
Collaborator

Ok, I did all the above fixes, plus a few internal tweaks. Also added the unittests I had forgotten to add.

Added complexity ratings, and merged a few docs with ditto's.

std/random.d
((191 lines not shown))
+ {
+ seed(value);
+ }
+
+/// ditto
+ this(Range)(Range range)
+ if (isInputRange!Range)
+ {
+ seed(range);
+ }
+
+/**
+ Tells if the generator has been seeded/initialized
+ */
+ @property @safe nothrow const
+ bool seeded()
@monarchdodra Collaborator

Is this a good name? Maybe "isSeeded" instead?

Property?

@andralex Owner

seeded is fine as long as it's a property

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
@@ -1008,6 +1008,443 @@ unittest
}
}
+/**
+ * Lagged Fibonacci generator.
+ *
+ * The single Engine can be started with either an Unsigned Integral Type, or a Real Type.
@andralex Owner

Why all of the capitalized names? If they represent reified concepts you'd need to define them. (I don't think that's necessary though.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((39 lines not shown))
+ enum UIntType mask = ~(~cast(UIntType)0 << bits); //A mask that is the size of "bits"
+ enum UIntType maxPrivate = mask; //maxPrivate is needed to have a single documented max
+ }
+ else //static if (isFloatingPoint!Type)
+ {
+ alias Type RealType;
+ enum RealType two32 = 2.0^^32;
+ enum RealType divisor = 1.0/(2.0^^bits);
+ enum size_t kMax = bits/32; //Amount of whole "u32" types that can fit in a type of size bits
+ enum uint mask = ~((~cast(uint)0) << (bits%32)); //A mask of the remaining bits
+ enum RealType maxPrivate = 1; //maxPrivate is needed to have a single documented max
+ }
+
+ private:
+ // We use a payload for "cheap" copy/move
+ struct Payload
@andralex Owner
static
@monarchdodra Collaborator

I did not know what static meant in this context, so I inquired over at digitalmars.D.learn.

It turns out that it is not clear whether or not the "outer" feature applies for nested structs or not :
it doesn't work, but it could just be "not yet implemented".

May I request your participation?
http://forum.dlang.org/thread/cwohbwozowasyfpkcoim@forum.dlang.org

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((50 lines not shown))
+ }
+
+ private:
+ // We use a payload for "cheap" copy/move
+ struct Payload
+ {
+ Type[longLag] x;
+ size_t i = 0;
+
+ void throwSeedException() const
+ {
+ throw new Exception(text("LaggedFibonacciEngine.seed: Input range didn't provide"
+ " enough elements: Need ", longLag, " elements."));
+ }
+
+ void seed(Range)(Range range)
@andralex Owner

constraint?

@monarchdodra Collaborator

This method is internal to the private Payload struct, the constraint should already be validated in the public "seed(*)" methods. It could warrant an internal static assert for sanity though

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@andralex andralex commented on the diff
std/random.d
((92 lines not shown))
+ range.popFront();
+ mult *= two32;
+ }
+ static if(mask != 0)
+ {
+ if(range.empty) throwSeedException();
+ val += (cast(uint)range.front & mask) * mult;
+ range.popFront();
+ }
+ assert(val >= 0.0, text("LaggedFibonacciEngine: Floating point value is smaller than 0: ", val));
+ assert(val < 1.0, text("LaggedFibonacciEngine: Floating point value is bigger than 1: ", val));
+ x[j] = val;
+ }
+ }
+ fill();
+ }
@andralex Owner

Since the bodies are so different it would make sense to define two functions with different constraints.

@monarchdodra Collaborator

Actually, it is really just an implementation detail about how the internal array is filled.
The constraints are the same, and they both finish with a call to "fill". I'd rather keep the static if inside.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((123 lines not shown))
+ foreach(j; 0..shortLag)
+ {
+ RealType t = x[j] + x[j+(longLag-shortLag)];
+ if(t >= 1.0)
+ t -= 1.0;
+ x[j] = t;
+ }
+ foreach(j; shortLag..longLag)
+ {
+ RealType t = x[j] + x[j-shortLag];
+ if(t >= 1.0)
+ t -= 1.0;
+ x[j] = t;
+ }
+ }
+ }
@andralex Owner

ditto

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((151 lines not shown))
+ {
+ fill();
+ i = 0;
+ }
+ }
+
+ @safe nothrow
+ void discard(size_t n)
+ {
+ while(n >= longLag)
+ {
+ fill();
+ n -= longLag;
+ }
+ i += n;
+ if(i >= longLag)
@andralex Owner

Not clear on the context, but shouldn't this be while?

@monarchdodra Collaborator

The while is already processed in "while(n >= longLag)".
After that while we are guaranteed that i + n < 2*longLag.

In theory, I could have just started with i+=n and looped on "i"... But there is an integer overflow risk there => Loop on "n", finish the last iteration on n + i.

I'll add an in-code comment for posterity.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@andralex andralex commented on the diff
std/random.d
((158 lines not shown))
+ void discard(size_t n)
+ {
+ while(n >= longLag)
+ {
+ fill();
+ n -= longLag;
+ }
+ i += n;
+ if(i >= longLag)
+ {
+ fill();
+ i -= longLag;
+ }
+ }
+ }
+ Payload* payload;
@andralex Owner

RefCounted!Payload?

@monarchdodra Collaborator

Isn't it simpler to just let the GC handle it? I don't think there is any reason to force RAII on the allocated space here.

@andralex Owner
andralex added a note

fine (this actually works great with the suggestion to make this a class)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((166 lines not shown))
+ if(i >= longLag)
+ {
+ fill();
+ i -= longLag;
+ }
+ }
+ }
+ Payload* payload;
+
+ public:
+
+/**
+ Creates a new $(D LaggedFibonacciEngine) and calls $(D seed()).
+*/
+ static
+ This opCall()
@andralex Owner

Unfortunately this is a bit confusing...

LaggedFibonacciEngine!(stuff) x; // not seeded
auto y = LaggedFibonacciEngine!(stuff)(); // seeded

Not sure whether this idiom will make it or not in the long run.

@monarchdodra Collaborator

Written like that, I agree, but like this:

LaggedFibonacci607 x; // not seeded (LaggedFibonacci607.init)
LaggedFibonacci607 y = LaggedFibonacci607(); // seeded
LaggedFibonacci607 y = LaggedFibonacci607(31); // seeded

This goes back to my comment that it would really be helpful if D allowed for "Constructor with no argument" (Not to be confused with "Default Constructor"). This would help in particular for all structs that have payload idioms: RefCounted, Array, etc...


But this is not the place to discuss this. Regarding Fibonacci, I see 3 options:
*1:Leave it as is.
*2:Remove ALL constuctors, and request an explicit call to seed.
*3:Make it a class. //Meh

Also, thinking a bit harder about it: Is it really the implementation's job to force reference semantics onto the generator? Shouldn't we just make the generator a BIG struct, but strongly hint that it should be allocated on the heap?

auto gen = new LaggedFibonacci607;
gen.seed();
gen.front;
...

I think it is a better approach here.


I don't have the experience/hindsight to answer these questions myself, so I'm deferring to your best judgement.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((215 lines not shown))
+ None: Uses a default seed.
+ T value: Uses a $(D MinstdRand0) generator seeded with value
+ InputRange range: Uses range.
+
+ Throws:
+ $(D Exception) if the InputRange didn't provide enough elements to seed the generator.
+ The number of elements required is the 'longLag' template parameter of the LaggedFibonacciEngine struct.
+*/
+ void seed()()
+ {
+ seed(defaultSeed);
+ }
+
+/// ditto
+ void seed(T)(T value)
+ if (isIntegral!T)
@andralex Owner

Let's add here the constraint that sizeof(value) <= 4

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((257 lines not shown))
+ }
+
+/// dito
+ alias dup save;
+
+/**
+ Returns true if the two generators will produce identical
+ sequences of outputs.
+
+ $(BIGOH 1) average case complexity.
+ $(BIGOH longLag) worst case complexity.
+ */
+ @safe nothrow const
+ bool opEquals(const(This) b)
+ {
+ return (payload == b.payload) ||
@andralex Owner

is?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@andralex
Owner

Getting there. Nice work!

@monarchdodra
Collaborator

Thank you again for the re-read. I did not have time to work on it this week, but I needed your feedback on the constructor/payload issue before continuing anyways. Regarding everything else, I think one more week of work, and everything should be final.

D is still a learning process for me, so I don't want to rush this.

@andralex
Owner

Just do what the Mersenne Twister does (i.e. ctor with seed and also seed method). We can always add the opCall trick if we decide it's worth it. For now let's aim for conservative uniformity. Thanks!

@monarchdodra
Collaborator

Fixed what needed to be fixed. Made everything as conservative as and straight forward as possible. I kept the reference semantics. It is kind of discussed in this thread:
http://forum.dlang.org/thread/spjfsyaqrcpboypgfrsd@forum.dlang.org
And the rationale in this post:
http://forum.dlang.org/thread/spjfsyaqrcpboypgfrsd@forum.dlang.org?page=3#post-bloacgckmnejwhbnjdlp:40forum.dlang.org

Long story short: It is much more convenient for the type itself to provide reference semantics, then for the developers to do it themselves. As a reminder, LaggedFibonacci44497 is 300K in size (!)

Anyways, I don't think there is much else to say here on my end.

I was still curious about nested structs and have a hidden "outer" pointer:
According to http://dlang.org/struct.html, they do.
According to TDPL, they don't.

Could we have your input in this discussion? Or here and I'll relay.
http://forum.dlang.org/thread/cwohbwozowasyfpkcoim@forum.dlang.org

std/random.d
@@ -1008,6 +1008,462 @@ unittest
}
}
+/**
+* Lagged Fibonacci generator.
+*
+* The single engine can be started with either an unsigned integral type,
+* or a floating point type. Unsigned Integral Type engine will generate
+* numbers in the range [0 .. 2^^bits). Real Type will generate
+* floating points in the range [0 .. 1), with a step of 1/2^^bits.
+*
+* The engine uses reference semantics, and models a ForwardRange.
@andralex Owner
andralex added a note

s/ForwardRange/$(D ForwardRange)/

@andralex Owner
andralex added a note

The fact that it uses reference semantics must be emphasize and discussed with examples. Documentation should also mention that this semantics is different from all other types in this module.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
@@ -1008,6 +1008,462 @@ unittest
}
}
+/**
+* Lagged Fibonacci generator.
+*
+* The single engine can be started with either an unsigned integral type,
+* or a floating point type. Unsigned Integral Type engine will generate
+* numbers in the range [0 .. 2^^bits). Real Type will generate
+* floating points in the range [0 .. 1), with a step of 1/2^^bits.
+*
+* The engine uses reference semantics, and models a ForwardRange.
+*
+* The engine must be seeded prior to use.
+*
+* Ported from boost 1.50 on July 30th 2012
+*/
+struct LaggedFibonacciEngine(Type, size_t bits, size_t longLag, size_t shortLag)
@andralex Owner
andralex added a note

Since the working set of LaggedFibonacciEngine is so large, it may be worthwhile just making it a class. That would simplify implementation considerably. Make methods final.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@andralex andralex commented on the diff
std/random.d
((18 lines not shown))
+struct LaggedFibonacciEngine(Type, size_t bits, size_t longLag, size_t shortLag)
+ if (isUnsigned!Type || isFloatingPoint!Type)
+{
+ public:
+ /// Mark this as a Rng
+ enum bool isUniformRandom = true;
+ /// Always $(D false) (random generators are infinite ranges).
+ enum empty = false;
+ /// Returns the smallest value that the generator can produce.
+ enum Type min = 0;
+ /// Returns the largest value that the generator can produce.
+ enum Type max = maxPrivate;
+
+ private:
+ enum uint defaultSeed = 331u;
+ alias typeof(this) This; //A *non-qualified* alias for typof(this).
@andralex Owner
andralex added a note

typo :o)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
std/random.d
((106 lines not shown))
+ assert(val >= 0.0, text("LaggedFibonacciEngine: Floating point value is smaller than 0: ", val));
+ assert(val < 1.0, text("LaggedFibonacciEngine: Floating point value is bigger than 1: ", val));
+ x[j] = val;
+ }
+ }
+ fill();
+ }
+
+ @safe nothrow
+ void fill()
+ {
+ static if (isUnsigned!Type)
+ {
+ // two loops to avoid costly modulo operations
+ foreach(j; 0..shortLag)
+ x[j] = (x[j] + x[j+(longLag-shortLag)]) & mask;
@andralex Owner
andralex added a note

I'm seeing inconsistent use of operator spacing, even in the same line. I'll use this as a pretext to require spacing around all operators, like in the rest of phobos.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@monarchdodra
Collaborator

Closed for development

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
Showing with 457 additions and 0 deletions.
  1. +457 −0 std/random.d
View
457 std/random.d
@@ -1008,6 +1008,463 @@ unittest
}
}
+/**
+* Lagged Fibonacci generator.
+*
+* The engine can be instanciated with either an unsigned integral type,
+* or a floating point type. An unsigned integral engine will generate
+* numbers in the range [0 .. 2^^bits). A floating point engine will generate
+* floating points in the range [0 .. 1), with a step of 1/2^^bits.
+*
+* A $(D LaggedFibonacciEngine) must maintain a very large set of data
+* to function. Because of this, and contrary to all other generators
+* in this module, the engine is allocated on the heap, and uses
+* reference semantics.
+*
+* The engine must be seeded prior to use. Use of an un-seeded
+* (un-initialized) engine will generate an assertion error.
+*
+* The $(D LaggedFibonacciEngine) models an $(D Infinite)
+* $(D ForwardRange).
+*
+* Ported from boost 1.50 on July 30th 2012
+*/
+struct LaggedFibonacciEngine(Type, size_t bits, size_t longLag, size_t shortLag)
+ if (isUnsigned!Type || isFloatingPoint!Type)
+{
+ public:
+ /// Mark this as a Rng
+ enum bool isUniformRandom = true;
+ /// Always $(D false) (random generators are infinite ranges).
+ enum empty = false;
+ /// Returns the smallest value that the generator can produce.
+ enum Type min = 0;
+ /// Returns the largest value that the generator can produce.
+ enum Type max = maxPrivate;
+
+ private:
+ enum uint defaultSeed = 331u;
+ alias typeof(this) This; //A *non-qualified* alias for typof(this).
@andralex Owner
andralex added a note

typo :o)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+ //Conditional Aliases/enums/asserts
+ static if (isUnsigned!Type)
+ {
+ static assert(UIntType.sizeof*8 >= bits, text("Provided UIntType not large enough to hold ", bits, "-bit integers."));
+
+ alias Type UIntType;
+ enum size_t kMax = (bits+31)/32; //Total amount of "u32" types needed to fill a type of size bits
+ enum UIntType mask = cast(UIntType)~(~cast(UIntType)0 << bits); //A mask that is the size of "bits"
+ enum UIntType maxPrivate = mask; //maxPrivate is needed to have a single documented max
+ }
+ else //static if (isFloatingPoint!Type)
+ {
+ alias Type RealType;
+ enum RealType two32 = 2.0^^32;
+ enum RealType divisor = 1.0/(2.0^^bits);
+ enum size_t kMax = bits/32; //Amount of whole "u32" types that can fit in a type of size bits
+ enum uint mask = ~((~cast(uint)0) << (bits%32)); //A mask of the remaining bits
+ enum RealType maxPrivate = 1; //maxPrivate is needed to have a single documented max
+ }
+
+ private:
+ //Payload; Reference semantics
+ static struct Payload
+ {
+ Type[longLag] x;
+ size_t i = 0;
+
+ void throwSeedException() const
+ {
+ throw new Exception(
+ text("LaggedFibonacciEngine.seed: Input range didn't provide enough elements: Need ",
+ longLag, " elements."));
+ }
+
+ void seed(Range)(Range range)
+ //Range conditions were validated in public interface
+ {
+ i = 0;
+ static if (isUnsigned!Type)
+ {
+ foreach(j; 0..longLag)
+ {
+ UIntType val = 0;
+ foreach(k; 0..kMax)
+ {
+ if(range.empty) throwSeedException();
+ val += cast(UIntType)range.front << 32*k;
+ range.popFront();
+ }
+ x[j] = val & mask;
+ }
+ }
+ else //static if (isFloatingPoint!Type)
+ {
+ foreach(j; 0..longLag)
+ {
+ RealType val = 0;
+ RealType mult = divisor;
+ foreach(k; 0..kMax)
+ {
+ if(range.empty) throwSeedException();
+ val += cast(uint)range.front * mult;
+ range.popFront();
+ mult *= two32;
+ }
+ static if(mask != 0)
+ {
+ if(range.empty) throwSeedException();
+ val += (cast(uint)range.front & mask) * mult;
+ range.popFront();
+ }
+ assert(val >= 0.0, text("LaggedFibonacciEngine: Floating point value is smaller than 0: ", val));
+ assert(val < 1.0, text("LaggedFibonacciEngine: Floating point value is bigger than 1: ", val));
+ x[j] = val;
+ }
+ }
+ fill();
+ }
@andralex Owner

Since the bodies are so different it would make sense to define two functions with different constraints.

@monarchdodra Collaborator

Actually, it is really just an implementation detail about how the internal array is filled.
The constraints are the same, and they both finish with a call to "fill". I'd rather keep the static if inside.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+ @safe nothrow
+ void fill()
+ {
+ //The basic algorithm might not be visible at first.
+ //
+ //It is the same for both integer/float:
+ //foreach(i; 0..longLag)
+ //{
+ // x[i] += x[i - shortLag]; //*with modulo wrap around if negative*
+ // x[i].ApplyMask
+ //}
+ //Where ApplyMask is
+ //*Integer: "x[i] &= mask"
+ //*Floats: "if(x[i]>1) x[i] -= 1"
+ //
+ //We exploit vectorial operations for perfomance.
+ //Data is iterated on in strips of shortLag length.
+ //Entire "stips at once" are processed in addAndMask.
+ //There is a special "first iteration" case to avoid the costly modulo operation
+ //There is a special "final iteration" for the last (shorter) strip
+ @safe nothrow void addAndMask(Type[] r1, Type[] r2)
+ {
+ try //@@@ 8651
+ {
+ r1[] += r2[];
+ static if (isUnsigned!Type) r1[] &= mask;
+ else foreach(ref t; r1[]) if(t >= 1.0) t -= 1.0;
+ }
+ catch(Exception) assert(0, "Unexpected exception in LaggedFibonacciEngine.Payload.Fill");
+ }
+
+ addAndMask(x[0 .. shortLag], x[longLag - shortLag .. longLag]);
+ size_t low = shortLag;
+ for( ; low + shortLag <= longLag ; low += shortLag) //Iteration necessarry to avoid overlap operations
+ addAndMask(x[low .. low + shortLag], x[low - shortLag .. low]);
+ addAndMask(x[low .. longLag], x[low - shortLag .. longLag - shortLag]);
+ }
+
+ @property @safe nothrow const
+ Type front()
+ {
+ return x[i];
+ }
+
+ @safe nothrow
+ void popFront()
+ {
+ ++i;
+ if(i == longLag)
+ {
+ fill();
+ i = 0;
+ }
+ }
+
+ @safe nothrow
+ void discard(size_t n)
+ {
+ //Do NOT replace with n+=i, or risk an integer overflow
+ while(n >= longLag)
+ {
+ fill();
+ n -= longLag;
+ }
+ //Now that n is smaller, we can take i into account and do the last iteration
+ i += n;
+ if(i >= longLag)
+ {
+ fill();
+ i -= longLag;
+ }
+ }
+ }
+ Payload* payload;
@andralex Owner

RefCounted!Payload?

@monarchdodra Collaborator

Isn't it simpler to just let the GC handle it? I don't think there is any reason to force RAII on the allocated space here.

@andralex Owner
andralex added a note

fine (this actually works great with the suggestion to make this a class)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+ public:
+/**
+Tells if the generator has been seeded/initialized
+*/
+ @property @safe nothrow const
+ bool seeded()
+ {
+ return payload != null;
+ }
+
+/**
+Seeds the $(D LaggedFibonacciEngine)
+Parameters:
+None: Uses a default seed.
+T value: Uses a $(D MinstdRand0) generator seeded with value
+InputRange range: Uses range.
+Throws:
+$(D Exception) if the InputRange didn't provide enough elements to seed the generator.
+The number of elements required is the 'longLag' template parameter of the LaggedFibonacciEngine struct.
+*/
+ void seed()()
+ {
+ seed(defaultSeed);
+ }
+
+/// ditto
+ void seed()(uint value)
+ {
+ auto gen = MinstdRand0(value);
+ seed(gen);
+ }
+
+/// ditto
+ void seed(InputRange)(InputRange range)
+ if ( isInputRange!InputRange && isIntegral!(ElementType!InputRange) &&
+ ( ElementType!InputRange.sizeof >= 4 || ElementType!InputRange.sizeof >= Type.sizeof ) )
+ {
+ if(!payload) payload = new Payload;
+ payload.seed(range);
+ }
+
+/**
+Creates a deep copy of the generator; Saves the range.
+$(BIGOH longLag) complexity.
+*/
+ @property const
+ auto dup()
+ {
+ This ret; //non qualified typeof(this)
+ if(payload)
+ {
+ ret.payload = new Payload;
+ *ret.payload = *payload;
+ }
+ return ret;
+ }
+
+/// dito
+ alias dup save;
+
+/**
+Returns true if the two generators will produce identical
+sequences of outputs.
+$(BIGOH 1) average case complexity.
+$(BIGOH longLag) worst case complexity.
+*/
+ @safe nothrow const
+ bool opEquals(const(This) b)
+ {
+ return (payload is b.payload) ||
+ (payload && b.payload && *payload == *b.payload);
+ }
+
+ private enum payloadErrorText = "LaggedFibonacciEngine: Attempt to use an un-initialized engine";
+
+/**
+Returns the current value of the generator. Constant complexity.
+*/
+ @property @safe nothrow const
+ Type front()
+ {
+ assert(payload, payloadErrorText);
+ return payload.front;
+ }
+
+/**
+Advances to the next value of the generator.
+
+Amortized constant time: $(BIGOH longLag) in case of generator refill,
+$(BIGOH 1) otherwise. The generator is refilled every time the generator
+is advanced by $(D longLag).
+*/
+ @safe nothrow
+ void popFront()
+ {
+ assert(payload, payloadErrorText);
+ payload.popFront();
+ }
+
+/**
+Advances the state of the generator by $(D z). More efficient than calling popFront $(D z) times.
+
+Amortized $(BIGOH z) time: advancing the get pointer is constant time, but each refill
+(every time the engine is advanced by $(D longLag)) will still require $(D longLag) time.
+*/
+ @safe nothrow
+ void discard(size_t z)
+ {
+ assert(payload, payloadErrorText);
+ payload.discard(z);
+ }
+}
+
+/**
+Define $(D_PARAM LaggedFibonacciEngine) generators with well-chosen
+parameters.
+$(D LaggedFibonacci) is an implementation chosen LaggedFibonacci engine.
+*/
+alias LaggedFibonacciEngine!(double, 48, 607, 273) LaggedFibonacci607;
+alias LaggedFibonacciEngine!(double, 48, 1279, 418) LaggedFibonacci1279; /// ditto
+alias LaggedFibonacciEngine!(double, 48, 2281, 1252) LaggedFibonacci2281; /// ditto
+alias LaggedFibonacciEngine!(double, 48, 3217, 576) LaggedFibonacci3217; /// ditto
+alias LaggedFibonacciEngine!(double, 48, 4423, 2098) LaggedFibonacci4423; /// ditto
+alias LaggedFibonacciEngine!(double, 48, 9689, 5502) laggedFibonacci9689; /// ditto
+alias LaggedFibonacciEngine!(double, 48, 19937, 9842) LaggedFibonacci19937; /// ditto
+alias LaggedFibonacciEngine!(double, 48, 23209, 13470) LaggedFibonacci23209; /// ditto
+alias LaggedFibonacciEngine!(double, 48, 44497, 21034) LaggedFibonacci44497; /// ditto
+alias LaggedFibonacci607 LaggedFibonacci;
+
+unittest
+{
+ //Data aquired from boost implementation
+
+ //First 5 numbers of the different implementations
+ immutable(long[][]) referenceInt = [
+ [111357645752581, 223546595107190, 276701472505536, 183354198929810, 259724424734707], // 607
+ [ 18516896918599, 266714042254566, 214677200432103, 238166017619873, 253427522088996], // 1279
+ [209321651636266, 234403323170313, 44940711525953, 270585935924073, 182951237378271], // 2281
+ [102247219917095, 111280010135238, 8455294141221, 268442649885962, 181826007198542], // 3217
+ [212931708246959, 67779908108345, 109786160258021, 197789232895085, 173746813492116], // 4423
+ [136413310818494, 11236510562599, 52970671921490, 193128799152157, 216674819310303], // 9689
+ [ 7126909551884, 86072826243616, 80675133038411, 251970503587899, 124423863440468], //19937
+ [146111616638474, 136512418917669, 259280540566850, 250350614477568, 143029744160684], //23209
+ [263558154538665, 124374555248232, 92765545546811, 244707388756702, 173997194527172], //44497
+ ];
+
+ immutable(double[][]) referenceFloat = [
+ [ 0.395622, 0.794197, 0.983041, 0.651405, 0.922727], // 607
+ [0.0657852, 0.947559, 0.762687, 0.846136, 0.900355], // 1279
+ [ 0.74366, 0.832768, 0.159661, 0.961314, 0.649973], // 2281
+ [ 0.363255, 0.395346, 0.0300392, 0.9537, 0.645976], // 3217
+ [ 0.756485, 0.240803, 0.390039, 0.702689, 0.617273], // 4423
+ [ 0.484637, 0.0399201, 0.18819, 0.686131, 0.769784], // 9689
+ [0.0253199, 0.305792, 0.286616, 0.895179, 0.442042], //19937
+ [ 0.519093, 0.48499, 0.92115, 0.889424, 0.508144], //23209
+ [ 0.936347, 0.441867, 0.329569, 0.869375, 0.618162], //44497
+ ];
+
+ //Numbers [5000-5005)
+ immutable(long[][]) referenceInt5000 = [
+ [266635669194050, 70515619676375, 83207904227719, 281373199363211, 116250967224354], // 607
+ [110135446738590, 36692226031177, 209483500649828, 268031999359600, 248891650590185], // 1279
+ [217776584741215, 215117240100973, 52618792839041, 22635584224685, 198265981921064], // 2281
+ [ 93750178568681, 165145932004335, 78808838225286, 232759037338400, 61858783708817], // 3217
+ [275211622257559, 20961152411805, 16486400488039, 239571949617653, 119147348429321], // 4423
+ [103616703846189, 134586041061516, 4244911899451, 5940954177811, 10074550564657], // 9689
+ [255367856135618, 277538191528067, 178667831440893, 22734204862614, 227099404643056], //19937
+ [121591596418969, 167750783695320, 32717268361943, 281360144972641, 255183279435442], //23209
+ [279851845097169, 249217869577016, 152361400411229, 261604351718748, 76994486632824], //44497
+ ];
+
+ immutable(double[][]) referenceFloat5000 = [
+ [ 0.94728, 0.250522, 0.295614, 0.999638, 0.413006], // 607
+ [ 0.39128, 0.130357, 0.744235, 0.952241, 0.884241], // 1279
+ [0.773698, 0.76425, 0.18694, 0.0804177, 0.704382], // 2281
+ [0.333068, 0.586716, 0.279985, 0.826926, 0.219767], // 3217
+ [0.977748, 0.074469, 0.0585715, 0.851131, 0.423296], // 4423
+ [ 0.36812, 0.478146, 0.015081, 0.0211065, 0.035792], // 9689
+ [0.907249, 0.986014, 0.634756, 0.0807681, 0.806819], //19937
+ [ 0.43198, 0.595971, 0.116235, 0.999592, 0.906593], //23209
+ [0.994233, 0.8854, 0.541296, 0.929405, 0.273539], //44497
+ ];
+
+ void doTests(T, V)(V values, V values5000)
+ {
+ foreach (I, Type; TypeTuple!(
+ LaggedFibonacciEngine!(T, 48, 607, 273),
+ LaggedFibonacciEngine!(T, 48, 1279, 418),
+ LaggedFibonacciEngine!(T, 48, 2281, 1252),
+ LaggedFibonacciEngine!(T, 48, 3217, 576),
+ LaggedFibonacciEngine!(T, 48, 4423, 2098),
+ LaggedFibonacciEngine!(T, 48, 9689, 5502),
+ LaggedFibonacciEngine!(T, 48, 19937, 9842),
+ LaggedFibonacciEngine!(T, 48, 23209, 13470),
+ LaggedFibonacciEngine!(T, 48, 44497, 21034)
+ )
+ )
+ {
+ auto r = Type();
+ assert(!r.seeded);
+ r.seed();
+ assert(r.seeded);
+
+ assert(equal!approxEqual(r.take(5), values[I]));
+ r.discard(4995);
+ assert(equal!approxEqual(r.take(5), values5000[I]));
+
+ r.seed(); //check re-seed and save
+ assert(equal!approxEqual(r.save.take(5), values[I]));
+ r.discard(5000);
+ assert(equal!approxEqual(r.save.take(5), values5000[I]));
+
+ r.discard(5000);
+ assert(r == r);
+ assert(equal!approxEqual(r.save.take(5), r.save.take(5)));
+ auto ra = Type();
+ auto rb = Type();
+ ra.seed();
+ rb.seed();
+ assert(ra !is rb);
+ assert(ra == rb);
+ ra.popFront();
+ assert(ra != rb);
+ }
+
+ //Verify validity of output "steps"
+ auto fourGenerator = LaggedFibonacciEngine!(T, 4, 607, 273)();
+ fourGenerator.seed(unpredictableSeed());
+ T[] sixteen = [ 0, 1, 2, 3,
+ 4, 5, 6, 7,
+ 8, 9, 10, 11,
+ 12, 13, 14, 15];
+ static if(is(T == double))
+ {
+ sixteen[] /= 16.0;
+ }
+
+ foreach(v; fourGenerator.take(10))
+ {
+ sixteen.canFind!approxEqual(v);
+ }
+ }
+
+ doTests!(ulong)(referenceInt, referenceInt5000);
+ doTests!(double)(referenceFloat, referenceFloat5000);
+}
+unittest
+{
+ //Test everything works for an 8 bit generator
+ auto a = LaggedFibonacciEngine!(ushort, 8, 607, 273)();
+ ushort b = 5;
+ ulong c = 5;
+ a.seed(b); //Can be built seeded from a ushort
+ a.seed(cast(uint)c); //ulong needs a cast
+ a.seed(unpredictableSeed()); //unperdictable seed is fine
+}
+unittest
+{
+ auto a = LaggedFibonacci();
+ assertThrown(a.seed([1, 2, 3])); //Input range too short
+}
+unittest
+{
+ auto p = new LaggedFibonacci();
+}
/**
A "good" seed for initializing random number engines. Initializing
Something went wrong with that request. Please try again.