You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by lu...@free.fr on 2010/09/20 12:30:48 UTC
[math] new pseudo random number generators
Hi all,
As those subscribed to the commit list may have noticed, I have added several new pseudo-random number generators to commons-math.
I should not have add such a feature without discussing on the list before. It was a wrong move and I apologize about it.
So I would like to start a discussion about this now.
Up to 2.1, we have few pseudo random number generators. We have an interface RandomGenerator implemented by three classes:
- JDKRandomGenerator that extends the JDK provided generator
- AbstractRandomGenerator as a helper for users generators
- BitStreamGenerator which in turn is extended only by MersenneTwister
The JDK provided generator is a simple one that can be used only for very simple needs. The Mersenne Twister is on the other hand a fast generator with good properties well suited for Monte-Carlo simulation. It is equidistributed for generating vectors up to dimension 623 and has a huge period: 2^19937 - 1.
Since Mersenne-Twister inception in 1997, some new generators have been created, retaining the good properties of Mersenne twister but removing some of its (few) drawbacks. The main one is that if initialized with a bits pool containing lots of zeroes, the pool will take a very long time time to stabilize with a roughly balanced number of zeros and ones.
I would like to add such generators (well, I already did but can withdraw my commit). The ones I want to add are the WELL generators (Well Equidistributed Long period Linear) created by François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. They are described in their 2006 paper: Improved
* Long-Period Generators Based on Linear Recurrences Modulo 2, ransactions on Mathematical Software, 32, 1 (2006) which is available at <http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf>.
The paper describes several generators from one family. The most interesting ones are WELL1024 (a small one), WELL19937 and WELL44497 (two large ones with huge periods that are both Mersenne primes). The two large ones exist in a basic version that is not Maximally Equidistributed (i.e. equidistributed in all possible dimensions and number of bits up to some threshold) and also in an extended version that add some tempering to the output to get an ME generator.
The paper does not put any restriction on the algorithm. The reference implementation in C on the other hand is limited to non-commercial use only. In my implementation, I did not refer to this C implementation at all. I only compiled it on my personal computer to generate reference values and copied the values in the test cases. So this implementation is completely independent (and in fact the code is rather different since I used an abstract class for the global algorithm and one derived class for each specific generator. I rely on the JVM optimizer to do all the inlining and constants simplification that they did manually in the reference implementation.
Another difference with the reference implementation is that our BitStreamGenerator abstract class requires the generator to provide only bits blocks and not double, and it computes the double by itself. This computation uses the full IEEE754 width for the doubles, i.e. it uses 53 bits. The generators are guaranteed to be equidistributed on large vectors of 32 bits (up to dimension 1390 for the 44497 version, since 1390 * 32 < 44497). I guess this should work well also using chunks of 53 bits (but of course with smaller dimensions), but it is not mathematically proven.
So what do other think about this ?
Should this really be included (in which case I will open a JIRA issue and resolve it immediately) or should it be removed from commons-math ? Is the independent implementation and the off-line reference data generation sufficient to say this code is unencumbered ? Is the 53 bits use for double generation a good thing or should we reduce it to 32 bits to stay within the proven behavior ?
Once again, sorry for this too late discussion, I really should have started it before.
Luc
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Ted Dunning <te...@gmail.com>.
On Mon, Sep 20, 2010 at 6:16 AM, Gilles Sadowski <
gilles@harfang.homelinux.org> wrote:
> >
> > Sounds OK, but we should verify somehow that the algorithm itself is
> > not patented. [...]
>
> This is probably not the place to start a discussion on the difference
> between patent and copyright :-}
> Copyright can apply to source code; neither patent nor copyright should
> apply to algorithms (even though American patent law has been subverted
> in order to include them).
>
> Does the Apache foundation give in to this "IP" bull<snip> ?
Apache generally tries to avoid trouble by settling the questions at
development time. That isn't the same as "giving in", but you could read it
that way. I prefer to think of it as having a different goal.
Re: [math] new pseudo random number generators
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
> [...]
>
> Sounds OK, but we should verify somehow that the algorithm itself is
> not patented. [...]
This is probably not the place to start a discussion on the difference
between patent and copyright :-}
Copyright can apply to source code; neither patent nor copyright should
apply to algorithms (even though American patent law has been subverted
in order to include them).
Does the Apache foundation give in to this "IP" bull<snip> ?
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
> [...]
> >
> > They do not use this kind of tricks in their implementation. By using
> > them in ours, it would make more clear our code is really an original
> > one and not derived from their GPL code.
> It sounds reasonable, but I'm thinking that it might be better to add
> the current code, and then optimise it further afterwards? [...]
+1
[Not being sure that the speed gain provided by the "optimized/obfuscated"
code will not be dwarfed by e.g. the next version of the JVM.]
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Phil Steitz <ph...@gmail.com>.
On 9/21/10 6:40 AM, sebb wrote:
> On 21 September 2010 09:40, Mikkel Meyer Andersen<mi...@mikl.dk> wrote:
>> 2010/9/21 Luc Maisonobe<Lu...@free.fr>:
>>> Le 21/09/2010 10:08, Mikkel Meyer Andersen a écrit :
>>>> 2010/9/21 Luc Maisonobe<Lu...@free.fr>:
>>>>> Le 21/09/2010 09:26, Mikkel Meyer Andersen a écrit :
>>>>>>> Here is an extract from the answer from Pierre L'Ecuyer:
>>>>>>>
>>>>>>> Our code can be released under either a GPL or a commercial license.
>>>>>> Well, what about the Apache License then?
>>>>>
>>>>> It is their code and I did not use it. The main point is we
>>>>> reimplemented it. Of course, relying on their code that is GPLed would
>>>>> be a clear no-go as GPL is a category X license (see
>>>>> <http://www.apache.org/legal/resolved.html#category-x>).
>>>> Exactly. But as I understood, you used some of their optimization
>>>> principles - but that is not an issue or?
>>>
>>> I finally inlined everything, yes but I doubt it would count as a
>>> specific optimization found only in this original code. It's a common
>>> way and especially when the code to be inlined is one or two
>>> instructions only. The elementary transforms here are very simple
>>> (things like "v0 ^ (v0<< 25)" or "((z2& 0x00020000) != 0) ? (z2Prime ^
>>> 0xb729fcec) : z2Prime"). Inlining this is common sense.
>> Agree.
>>>
>>>>>
>>>>> Perhaps I should give a go to still other optimizations just to make
>>>>> more clear our code is really different ?
>>>> I'm not sure what you mean?
>>>
>>> Some ways to optimize further could be to process several numbers in raw
>>> to benefit from some of the values being already in register rather than
>>> putting them in a memory array to read them back just the next iteration
>>> (one Well iteration uses 5 values from the bits pool array, two being
>>> close together and the three other being spread in the array).
>>> Basically, it would be packing two or more iterations in one call and
>>> cache the results. Then, when the user asks for 1 million numbers, some
>>> calls would result in simply retrieving one value from the cache and
>>> only a few calls would mean cherry picking in the array and performing
>>> the bits twiddling.
>>>
>>> They do not use this kind of tricks in their implementation. By using
>>> them in ours, it would make more clear our code is really an original
>>> one and not derived from their GPL code.
>> It sounds reasonable, but I'm thinking that it might be better to add
>> the current code, and then optimise it further afterwards? As I see
>> it, that would provide some advantages (e.g. testing, performance etc.
>> for everybody).
>
> +1
>
> Also, it would help future maintenance if the optimisations were
> documented as comments in the code (rather than being buried in SVN
> history).
+1
Phil
>
>> Mikkel.
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>> For additional commands, e-mail: dev-help@commons.apache.org
>>
>>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by sebb <se...@gmail.com>.
On 21 September 2010 09:40, Mikkel Meyer Andersen <mi...@mikl.dk> wrote:
> 2010/9/21 Luc Maisonobe <Lu...@free.fr>:
>> Le 21/09/2010 10:08, Mikkel Meyer Andersen a écrit :
>>> 2010/9/21 Luc Maisonobe <Lu...@free.fr>:
>>>> Le 21/09/2010 09:26, Mikkel Meyer Andersen a écrit :
>>>>>> Here is an extract from the answer from Pierre L'Ecuyer:
>>>>>>
>>>>>> Our code can be released under either a GPL or a commercial license.
>>>>> Well, what about the Apache License then?
>>>>
>>>> It is their code and I did not use it. The main point is we
>>>> reimplemented it. Of course, relying on their code that is GPLed would
>>>> be a clear no-go as GPL is a category X license (see
>>>> <http://www.apache.org/legal/resolved.html#category-x>).
>>> Exactly. But as I understood, you used some of their optimization
>>> principles - but that is not an issue or?
>>
>> I finally inlined everything, yes but I doubt it would count as a
>> specific optimization found only in this original code. It's a common
>> way and especially when the code to be inlined is one or two
>> instructions only. The elementary transforms here are very simple
>> (things like "v0 ^ (v0 << 25)" or "((z2 & 0x00020000) != 0) ? (z2Prime ^
>> 0xb729fcec) : z2Prime"). Inlining this is common sense.
> Agree.
>>
>>>>
>>>> Perhaps I should give a go to still other optimizations just to make
>>>> more clear our code is really different ?
>>> I'm not sure what you mean?
>>
>> Some ways to optimize further could be to process several numbers in raw
>> to benefit from some of the values being already in register rather than
>> putting them in a memory array to read them back just the next iteration
>> (one Well iteration uses 5 values from the bits pool array, two being
>> close together and the three other being spread in the array).
>> Basically, it would be packing two or more iterations in one call and
>> cache the results. Then, when the user asks for 1 million numbers, some
>> calls would result in simply retrieving one value from the cache and
>> only a few calls would mean cherry picking in the array and performing
>> the bits twiddling.
>>
>> They do not use this kind of tricks in their implementation. By using
>> them in ours, it would make more clear our code is really an original
>> one and not derived from their GPL code.
> It sounds reasonable, but I'm thinking that it might be better to add
> the current code, and then optimise it further afterwards? As I see
> it, that would provide some advantages (e.g. testing, performance etc.
> for everybody).
+1
Also, it would help future maintenance if the optimisations were
documented as comments in the code (rather than being buried in SVN
history).
> Mikkel.
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Mikkel Meyer Andersen <mi...@mikl.dk>.
2010/9/21 Luc Maisonobe <Lu...@free.fr>:
> Le 21/09/2010 10:08, Mikkel Meyer Andersen a écrit :
>> 2010/9/21 Luc Maisonobe <Lu...@free.fr>:
>>> Le 21/09/2010 09:26, Mikkel Meyer Andersen a écrit :
>>>>> Here is an extract from the answer from Pierre L'Ecuyer:
>>>>>
>>>>> Our code can be released under either a GPL or a commercial license.
>>>> Well, what about the Apache License then?
>>>
>>> It is their code and I did not use it. The main point is we
>>> reimplemented it. Of course, relying on their code that is GPLed would
>>> be a clear no-go as GPL is a category X license (see
>>> <http://www.apache.org/legal/resolved.html#category-x>).
>> Exactly. But as I understood, you used some of their optimization
>> principles - but that is not an issue or?
>
> I finally inlined everything, yes but I doubt it would count as a
> specific optimization found only in this original code. It's a common
> way and especially when the code to be inlined is one or two
> instructions only. The elementary transforms here are very simple
> (things like "v0 ^ (v0 << 25)" or "((z2 & 0x00020000) != 0) ? (z2Prime ^
> 0xb729fcec) : z2Prime"). Inlining this is common sense.
Agree.
>
>>>
>>> Perhaps I should give a go to still other optimizations just to make
>>> more clear our code is really different ?
>> I'm not sure what you mean?
>
> Some ways to optimize further could be to process several numbers in raw
> to benefit from some of the values being already in register rather than
> putting them in a memory array to read them back just the next iteration
> (one Well iteration uses 5 values from the bits pool array, two being
> close together and the three other being spread in the array).
> Basically, it would be packing two or more iterations in one call and
> cache the results. Then, when the user asks for 1 million numbers, some
> calls would result in simply retrieving one value from the cache and
> only a few calls would mean cherry picking in the array and performing
> the bits twiddling.
>
> They do not use this kind of tricks in their implementation. By using
> them in ours, it would make more clear our code is really an original
> one and not derived from their GPL code.
It sounds reasonable, but I'm thinking that it might be better to add
the current code, and then optimise it further afterwards? As I see
it, that would provide some advantages (e.g. testing, performance etc.
for everybody).
Mikkel.
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Luc Maisonobe <Lu...@free.fr>.
Le 21/09/2010 10:08, Mikkel Meyer Andersen a écrit :
> 2010/9/21 Luc Maisonobe <Lu...@free.fr>:
>> Le 21/09/2010 09:26, Mikkel Meyer Andersen a écrit :
>>>> Here is an extract from the answer from Pierre L'Ecuyer:
>>>>
>>>> Our code can be released under either a GPL or a commercial license.
>>> Well, what about the Apache License then?
>>
>> It is their code and I did not use it. The main point is we
>> reimplemented it. Of course, relying on their code that is GPLed would
>> be a clear no-go as GPL is a category X license (see
>> <http://www.apache.org/legal/resolved.html#category-x>).
> Exactly. But as I understood, you used some of their optimization
> principles - but that is not an issue or?
I finally inlined everything, yes but I doubt it would count as a
specific optimization found only in this original code. It's a common
way and especially when the code to be inlined is one or two
instructions only. The elementary transforms here are very simple
(things like "v0 ^ (v0 << 25)" or "((z2 & 0x00020000) != 0) ? (z2Prime ^
0xb729fcec) : z2Prime"). Inlining this is common sense.
>>
>> Perhaps I should give a go to still other optimizations just to make
>> more clear our code is really different ?
> I'm not sure what you mean?
Some ways to optimize further could be to process several numbers in raw
to benefit from some of the values being already in register rather than
putting them in a memory array to read them back just the next iteration
(one Well iteration uses 5 values from the bits pool array, two being
close together and the three other being spread in the array).
Basically, it would be packing two or more iterations in one call and
cache the results. Then, when the user asks for 1 million numbers, some
calls would result in simply retrieving one value from the cache and
only a few calls would mean cherry picking in the array and performing
the bits twiddling.
They do not use this kind of tricks in their implementation. By using
them in ours, it would make more clear our code is really an original
one and not derived from their GPL code.
Luc
>>
>> Luc
>>
>>>
>>>> There is also a Java implementation with multiple streams and
>>>> substreams in SSJ: see the package rng:
>>>> http://www.iro.umontreal.ca/~simardr/ssj/indexe.html
>>>> If you reimplement the code, this is a gray zone, but we do not have
>>>> a patent on the algorithm.
>>>>
>>>> So as I reimplemented from the paper itself (taking the erratas in <
>>>> http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt>
>>>> into account) and use different optimization techniques (precomputed
>>>> indices tables), I would consider it is safe to publish this home-grown
>>>> code.
>>>
>>> ---------------------------------------------------------------------
>>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>>> For additional commands, e-mail: dev-help@commons.apache.org
>>>
>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>> For additional commands, e-mail: dev-help@commons.apache.org
>>
>>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Mikkel Meyer Andersen <mi...@mikl.dk>.
2010/9/21 Luc Maisonobe <Lu...@free.fr>:
> Le 21/09/2010 09:26, Mikkel Meyer Andersen a écrit :
>>> Here is an extract from the answer from Pierre L'Ecuyer:
>>>
>>> Our code can be released under either a GPL or a commercial license.
>> Well, what about the Apache License then?
>
> It is their code and I did not use it. The main point is we
> reimplemented it. Of course, relying on their code that is GPLed would
> be a clear no-go as GPL is a category X license (see
> <http://www.apache.org/legal/resolved.html#category-x>).
Exactly. But as I understood, you used some of their optimization
principles - but that is not an issue or?
>
> Perhaps I should give a go to still other optimizations just to make
> more clear our code is really different ?
I'm not sure what you mean?
>
> Luc
>
>>
>>> There is also a Java implementation with multiple streams and
>>> substreams in SSJ: see the package rng:
>>> http://www.iro.umontreal.ca/~simardr/ssj/indexe.html
>>> If you reimplement the code, this is a gray zone, but we do not have
>>> a patent on the algorithm.
>>>
>>> So as I reimplemented from the paper itself (taking the erratas in <
>>> http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt>
>>> into account) and use different optimization techniques (precomputed
>>> indices tables), I would consider it is safe to publish this home-grown
>>> code.
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>> For additional commands, e-mail: dev-help@commons.apache.org
>>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Phil Steitz <ph...@gmail.com>.
On 9/21/10 3:34 AM, Luc Maisonobe wrote:
> Le 21/09/2010 09:26, Mikkel Meyer Andersen a écrit :
>>> Here is an extract from the answer from Pierre L'Ecuyer:
>>>
>>> Our code can be released under either a GPL or a commercial license.
>> Well, what about the Apache License then?
>
> It is their code and I did not use it. The main point is we
> reimplemented it. Of course, relying on their code that is GPLed would
> be a clear no-go as GPL is a category X license (see
> <http://www.apache.org/legal/resolved.html#category-x>).
>
> Perhaps I should give a go to still other optimizations just to make
> more clear our code is really different ?
>
Thanks, Luc! We are probably OK, but best to check on legal-discuss.
Phil
> Luc
>
>>
>>> There is also a Java implementation with multiple streams and
>>> substreams in SSJ: see the package rng:
>>> http://www.iro.umontreal.ca/~simardr/ssj/indexe.html
>>> If you reimplement the code, this is a gray zone, but we do not have
>>> a patent on the algorithm.
>>>
>>> So as I reimplemented from the paper itself (taking the erratas in<
>>> http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt>
>>> into account) and use different optimization techniques (precomputed
>>> indices tables), I would consider it is safe to publish this home-grown
>>> code.
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>> For additional commands, e-mail: dev-help@commons.apache.org
>>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Luc Maisonobe <Lu...@free.fr>.
Le 21/09/2010 09:26, Mikkel Meyer Andersen a écrit :
>> Here is an extract from the answer from Pierre L'Ecuyer:
>>
>> Our code can be released under either a GPL or a commercial license.
> Well, what about the Apache License then?
It is their code and I did not use it. The main point is we
reimplemented it. Of course, relying on their code that is GPLed would
be a clear no-go as GPL is a category X license (see
<http://www.apache.org/legal/resolved.html#category-x>).
Perhaps I should give a go to still other optimizations just to make
more clear our code is really different ?
Luc
>
>> There is also a Java implementation with multiple streams and
>> substreams in SSJ: see the package rng:
>> http://www.iro.umontreal.ca/~simardr/ssj/indexe.html
>> If you reimplement the code, this is a gray zone, but we do not have
>> a patent on the algorithm.
>>
>> So as I reimplemented from the paper itself (taking the erratas in <
>> http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt>
>> into account) and use different optimization techniques (precomputed
>> indices tables), I would consider it is safe to publish this home-grown
>> code.
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Mikkel Meyer Andersen <mi...@mikl.dk>.
> Here is an extract from the answer from Pierre L'Ecuyer:
>
> Our code can be released under either a GPL or a commercial license.
Well, what about the Apache License then?
> There is also a Java implementation with multiple streams and
> substreams in SSJ: see the package rng:
> http://www.iro.umontreal.ca/~simardr/ssj/indexe.html
> If you reimplement the code, this is a gray zone, but we do not have
> a patent on the algorithm.
>
> So as I reimplemented from the paper itself (taking the erratas in <
> http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt>
> into account) and use different optimization techniques (precomputed
> indices tables), I would consider it is safe to publish this home-grown
> code.
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Luc Maisonobe <Lu...@free.fr>.
Le 20/09/2010 14:40, Phil Steitz a écrit :
> On 9/20/10 6:30 AM, luc.maisonobe@free.fr wrote:
>> The paper does not put any restriction on the algorithm. The
>> reference implementation in C on the other hand is limited to
>> non-commercial use only. In my implementation, I did not refer to
>> this C implementation at all. I only compiled it on my personal
>> computer to generate reference values and copied the values in
>> the test cases. So this implementation is completely independent
>> (and in fact the code is rather different since I used an
>> abstract class for the global algorithm and one derived class for
>> each specific generator. I rely on the JVM optimizer to do all
>> the inlining and constants simplification that they did manually
>> in the reference implementation.
>
> Sounds OK, but we should verify somehow that the algorithm itself is not
> patented. Does the C source say anything about that? Might be good to
> ask about this on legal-discuss or ask one of the authors. I know this
> could apply to any of the algorithms we implement, but the ones that are
>>100 yrs old are a little safer ;)
Here is an extract from the answer from Pierre L'Ecuyer:
Our code can be released under either a GPL or a commercial license.
There is also a Java implementation with multiple streams and
substreams in SSJ: see the package rng:
http://www.iro.umontreal.ca/~simardr/ssj/indexe.html
If you reimplement the code, this is a gray zone, but we do not have
a patent on the algorithm.
So as I reimplemented from the paper itself (taking the erratas in <
http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt>
into account) and use different optimization techniques (precomputed
indices tables), I would consider it is safe to publish this home-grown
code.
>
>>
>> Another difference with the reference implementation is that our
>> BitStreamGenerator abstract class requires the generator to
>> provide only bits blocks and not double, and it computes the
>> double by itself. This computation uses the full IEEE754 width
>> for the doubles, i.e. it uses 53 bits. The generators are
>> guaranteed to be equidistributed on large vectors of 32 bits (up
>> to dimension 1390 for the 44497 version, since 1390 * 32<
>> 44497). I guess this should work well also using chunks of 53
>> bits (but of course with smaller dimensions), but it is not
>> mathematically proven.
Pierre L'Ecuyer confirmed me we could do that despite equidistribution
is not mathematically proven in this case. In fact, he wrote me they do
the same ...
Luc
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by lu...@free.fr.
----- "Phil Steitz" <ph...@gmail.com> a écrit :
> On 9/20/10 6:30 AM, luc.maisonobe@free.fr wrote:
> > Hi all,
> >
> > As those subscribed to the commit list may have noticed, I have
> > added several new pseudo-random number generators to
> > commons-math.
> >
> > I should not have add such a feature without discussing on the
> > list before. It was a wrong move and I apologize about it. So I
> > would like to start a discussion about this now.
>
> Chalk it up to network latency ;)
> >
> > Up to 2.1, we have few pseudo random number generators. We have
> > an interface RandomGenerator implemented by three classes: -
> > JDKRandomGenerator that extends the JDK provided generator -
> > AbstractRandomGenerator as a helper for users generators -
> > BitStreamGenerator which in turn is extended only by
> > MersenneTwister
> >
> > The JDK provided generator is a simple one that can be used only
> > for very simple needs. The Mersenne Twister is on the other hand
> > a fast generator with good properties well suited for Monte-Carlo
> > simulation. It is equidistributed for generating vectors up to
> > dimension 623 and has a huge period: 2^19937 - 1.
> >
> > Since Mersenne-Twister inception in 1997, some new generators
> > have been created, retaining the good properties of Mersenne
> > twister but removing some of its (few) drawbacks. The main one is
> > that if initialized with a bits pool containing lots of zeroes,
> > the pool will take a very long time time to stabilize with a
> > roughly balanced number of zeros and ones.
> >
> > I would like to add such generators (well, I already did but can
> > withdraw my commit). The ones I want to add are the WELL
> > generators (Well Equidistributed Long period Linear) created by
> > François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. They are
> > described in their 2006 paper: Improved * Long-Period Generators
> > Based on Linear Recurrences Modulo 2, ransactions on Mathematical
> > Software, 32, 1 (2006) which is available
> > at<http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf>.
> >
> > The paper describes several generators from one family. The most
> > interesting ones are WELL1024 (a small one), WELL19937 and
> > WELL44497 (two large ones with huge periods that are both
> > Mersenne primes). The two large ones exist in a basic version
> > that is not Maximally Equidistributed (i.e. equidistributed in
> > all possible dimensions and number of bits up to some threshold)
> > and also in an extended version that add some tempering to the
> > output to get an ME generator.
> >
> > The paper does not put any restriction on the algorithm. The
> > reference implementation in C on the other hand is limited to
> > non-commercial use only. In my implementation, I did not refer to
> > this C implementation at all. I only compiled it on my personal
> > computer to generate reference values and copied the values in
> > the test cases. So this implementation is completely independent
> > (and in fact the code is rather different since I used an
> > abstract class for the global algorithm and one derived class for
> > each specific generator. I rely on the JVM optimizer to do all
> > the inlining and constants simplification that they did manually
> > in the reference implementation.
>
> Sounds OK, but we should verify somehow that the algorithm itself is
> not patented. Does the C source say anything about that? Might be
> good to ask about this on legal-discuss or ask one of the authors.
> I know this could apply to any of the algorithms we implement, but
> the ones that are >100 yrs old are a little safer ;)
OK, I'll write to the authors off-list and ask them. I also needed to ask them a question about some of the constants since there seem to be some typos in the tables and I want to be sure to use the ones that have the good properties we need.
Luc
>
> >
> > Another difference with the reference implementation is that our
> > BitStreamGenerator abstract class requires the generator to
> > provide only bits blocks and not double, and it computes the
> > double by itself. This computation uses the full IEEE754 width
> > for the doubles, i.e. it uses 53 bits. The generators are
> > guaranteed to be equidistributed on large vectors of 32 bits (up
> > to dimension 1390 for the 44497 version, since 1390 * 32<
> > 44497). I guess this should work well also using chunks of 53
> > bits (but of course with smaller dimensions), but it is not
> > mathematically proven.
> >
> > So what do other think about this ? Should this really be
> > included (in which case I will open a JIRA issue and resolve it
> > immediately) or should it be removed from commons-math ? Is the
> > independent implementation and the off-line reference data
> > generation sufficient to say this code is unencumbered ? Is the
> > 53 bits use for double generation a good thing or should we
> > reduce it to 32 bits to stay within the proven behavior ?
> >
>
> +1 for including assuming all is OK with IP
>
> Phil
>
>
>
> > Once again, sorry for this too late discussion, I really should
> > have started it before. Luc
> >
> >
> ---------------------------------------------------------------------
> >
> >
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > For additional commands, e-mail: dev-help@commons.apache.org
> >
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Phil Steitz <ph...@gmail.com>.
On 9/20/10 6:30 AM, luc.maisonobe@free.fr wrote:
> Hi all,
>
> As those subscribed to the commit list may have noticed, I have
> added several new pseudo-random number generators to
> commons-math.
>
> I should not have add such a feature without discussing on the
> list before. It was a wrong move and I apologize about it. So I
> would like to start a discussion about this now.
Chalk it up to network latency ;)
>
> Up to 2.1, we have few pseudo random number generators. We have
> an interface RandomGenerator implemented by three classes: -
> JDKRandomGenerator that extends the JDK provided generator -
> AbstractRandomGenerator as a helper for users generators -
> BitStreamGenerator which in turn is extended only by
> MersenneTwister
>
> The JDK provided generator is a simple one that can be used only
> for very simple needs. The Mersenne Twister is on the other hand
> a fast generator with good properties well suited for Monte-Carlo
> simulation. It is equidistributed for generating vectors up to
> dimension 623 and has a huge period: 2^19937 - 1.
>
> Since Mersenne-Twister inception in 1997, some new generators
> have been created, retaining the good properties of Mersenne
> twister but removing some of its (few) drawbacks. The main one is
> that if initialized with a bits pool containing lots of zeroes,
> the pool will take a very long time time to stabilize with a
> roughly balanced number of zeros and ones.
>
> I would like to add such generators (well, I already did but can
> withdraw my commit). The ones I want to add are the WELL
> generators (Well Equidistributed Long period Linear) created by
> François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. They are
> described in their 2006 paper: Improved * Long-Period Generators
> Based on Linear Recurrences Modulo 2, ransactions on Mathematical
> Software, 32, 1 (2006) which is available
> at<http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf>.
>
> The paper describes several generators from one family. The most
> interesting ones are WELL1024 (a small one), WELL19937 and
> WELL44497 (two large ones with huge periods that are both
> Mersenne primes). The two large ones exist in a basic version
> that is not Maximally Equidistributed (i.e. equidistributed in
> all possible dimensions and number of bits up to some threshold)
> and also in an extended version that add some tempering to the
> output to get an ME generator.
>
> The paper does not put any restriction on the algorithm. The
> reference implementation in C on the other hand is limited to
> non-commercial use only. In my implementation, I did not refer to
> this C implementation at all. I only compiled it on my personal
> computer to generate reference values and copied the values in
> the test cases. So this implementation is completely independent
> (and in fact the code is rather different since I used an
> abstract class for the global algorithm and one derived class for
> each specific generator. I rely on the JVM optimizer to do all
> the inlining and constants simplification that they did manually
> in the reference implementation.
Sounds OK, but we should verify somehow that the algorithm itself is
not patented. Does the C source say anything about that? Might be
good to ask about this on legal-discuss or ask one of the authors.
I know this could apply to any of the algorithms we implement, but
the ones that are >100 yrs old are a little safer ;)
>
> Another difference with the reference implementation is that our
> BitStreamGenerator abstract class requires the generator to
> provide only bits blocks and not double, and it computes the
> double by itself. This computation uses the full IEEE754 width
> for the doubles, i.e. it uses 53 bits. The generators are
> guaranteed to be equidistributed on large vectors of 32 bits (up
> to dimension 1390 for the 44497 version, since 1390 * 32<
> 44497). I guess this should work well also using chunks of 53
> bits (but of course with smaller dimensions), but it is not
> mathematically proven.
>
> So what do other think about this ? Should this really be
> included (in which case I will open a JIRA issue and resolve it
> immediately) or should it be removed from commons-math ? Is the
> independent implementation and the off-line reference data
> generation sufficient to say this code is unencumbered ? Is the
> 53 bits use for double generation a good thing or should we
> reduce it to 32 bits to stay within the proven behavior ?
>
+1 for including assuming all is OK with IP
Phil
> Once again, sorry for this too late discussion, I really should
> have started it before. Luc
>
> ---------------------------------------------------------------------
>
>
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by lu...@free.fr.
----- "Mikkel Meyer Andersen" <mi...@mikl.dk> a écrit :
> Hi,
>
> I think it sounds interesting.
>
> A few comments:
> 1) Does it support paralleled generation?
No, it updates its bits pool using an iterative formula.
> 2) In regards to the user, which should be the default for e.g.
> sampling from probability distributions?
For simple sampling, I guess WELL1024a could be sufficient. For Monte-Carlo simulations with more than 32 independent variables (i.e. dimensions), it would probably be insufficient so WELL19937c or WELL44497b would be better. See also answer to the following point.
> 3) How is the performance compared to the existing algorithms in
> Commons Math?
I did not check yet, but the paper shows benchmark with performances of the same order of magnitude as the MersenneTwister. All these are really fast generators (between 30 and 40 seconds to generate one billion numbers on a 2006 computer). These benchmark also show all generators ("short" periods like 2^1024-1 and long periods like 2^44497-1) have similar performances in term of generation rate. There is a difference in memory consumption, but this size is not a real concern unless you use hundreds of thousands of generators at the same time. The pool size of WELL44497b is 44497 bits stored in a 1391 elements integer array, so less than 6 kilobytes ...
Beware my implementation is really different and relies a lot on JVM optimizations (mainly inlining since there is a two-levels indirection for each of the 7 basic transforms t0 to t6). We may have a huge performance drop compared to the reference implementation if I did something wrong!
Although I forgot to point one important thing: these generators, just like our existing MersenneTwister are NOT suitable for cryptography. They are devoted to simulation, and to generate very long series with strong properties on the series as a whole (equidistribution, no correlation ...). They do not attempt to create small series but with very strong properties of unpredictability (does this word exist ?). So by adding these generators to our library we do not change its status with respect to exportation laws that apply in the US (well, and I developed them in France too).
Luc
>
> Cheers, Mikkel.
>
> 2010/9/20 <lu...@free.fr>:
> > Hi all,
> >
> > As those subscribed to the commit list may have noticed, I have
> added several new pseudo-random number generators to commons-math.
> >
> > I should not have add such a feature without discussing on the list
> before. It was a wrong move and I apologize about it.
> > So I would like to start a discussion about this now.
> >
> > Up to 2.1, we have few pseudo random number generators. We have an
> interface RandomGenerator implemented by three classes:
> > - JDKRandomGenerator that extends the JDK provided generator
> > - AbstractRandomGenerator as a helper for users generators
> > - BitStreamGenerator which in turn is extended only by
> MersenneTwister
> >
> > The JDK provided generator is a simple one that can be used only for
> very simple needs. The Mersenne Twister is on the other hand a fast
> generator with good properties well suited for Monte-Carlo simulation.
> It is equidistributed for generating vectors up to dimension 623 and
> has a huge period: 2^19937 - 1.
> >
> > Since Mersenne-Twister inception in 1997, some new generators have
> been created, retaining the good properties of Mersenne twister but
> removing some of its (few) drawbacks. The main one is that if
> initialized with a bits pool containing lots of zeroes, the pool will
> take a very long time time to stabilize with a roughly balanced number
> of zeros and ones.
> >
> > I would like to add such generators (well, I already did but can
> withdraw my commit). The ones I want to add are the WELL generators
> (Well Equidistributed Long period Linear) created by François
> Panneton, Pierre L'Ecuyer and Makoto Matsumoto. They are described in
> their 2006 paper: Improved
> > * Long-Period Generators Based on Linear Recurrences Modulo 2,
> ransactions on Mathematical Software, 32, 1 (2006) which is available
> at <http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf>.
> >
> > The paper describes several generators from one family. The most
> interesting ones are WELL1024 (a small one), WELL19937 and WELL44497
> (two large ones with huge periods that are both Mersenne primes). The
> two large ones exist in a basic version that is not Maximally
> Equidistributed (i.e. equidistributed in all possible dimensions and
> number of bits up to some threshold) and also in an extended version
> that add some tempering to the output to get an ME generator.
> >
> > The paper does not put any restriction on the algorithm. The
> reference implementation in C on the other hand is limited to
> non-commercial use only. In my implementation, I did not refer to this
> C implementation at all. I only compiled it on my personal computer to
> generate reference values and copied the values in the test cases. So
> this implementation is completely independent (and in fact the code is
> rather different since I used an abstract class for the global
> algorithm and one derived class for each specific generator. I rely on
> the JVM optimizer to do all the inlining and constants simplification
> that they did manually in the reference implementation.
> >
> > Another difference with the reference implementation is that our
> BitStreamGenerator abstract class requires the generator to provide
> only bits blocks and not double, and it computes the double by itself.
> This computation uses the full IEEE754 width for the doubles, i.e. it
> uses 53 bits. The generators are guaranteed to be equidistributed on
> large vectors of 32 bits (up to dimension 1390 for the 44497 version,
> since 1390 * 32 < 44497). I guess this should work well also using
> chunks of 53 bits (but of course with smaller dimensions), but it is
> not mathematically proven.
> >
> > So what do other think about this ?
> > Should this really be included (in which case I will open a JIRA
> issue and resolve it immediately) or should it be removed from
> commons-math ? Is the independent implementation and the off-line
> reference data generation sufficient to say this code is unencumbered
> ? Is the 53 bits use for double generation a good thing or should we
> reduce it to 32 bits to stay within the proven behavior ?
> >
> > Once again, sorry for this too late discussion, I really should have
> started it before.
> > Luc
> >
> >
> ---------------------------------------------------------------------
> > To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > For additional commands, e-mail: dev-help@commons.apache.org
> >
> >
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Phil Steitz <ph...@gmail.com>.
On 9/20/10 6:58 AM, Gilles Sadowski wrote:
>
>> [...]
>> 2) In regards to the user, which should be the default for e.g.
>> sampling from probability distributions?
>
> I have the same question for possibly replacing the default RNG in class
> "UnitSphereRandomVectorGenerator".
>
>
Very good point. We made all of the PRNG stuff pluggable so users
could use whatever they want; but we should probably not just always
default to the JDK generator.
Phil
> Gilles
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
> [...]
> 2) In regards to the user, which should be the default for e.g.
> sampling from probability distributions?
I have the same question for possibly replacing the default RNG in class
"UnitSphereRandomVectorGenerator".
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Mikkel Meyer Andersen <mi...@mikl.dk>.
Hi,
I think it sounds interesting.
A few comments:
1) Does it support paralleled generation?
2) In regards to the user, which should be the default for e.g.
sampling from probability distributions?
3) How is the performance compared to the existing algorithms in Commons Math?
Cheers, Mikkel.
2010/9/20 <lu...@free.fr>:
> Hi all,
>
> As those subscribed to the commit list may have noticed, I have added several new pseudo-random number generators to commons-math.
>
> I should not have add such a feature without discussing on the list before. It was a wrong move and I apologize about it.
> So I would like to start a discussion about this now.
>
> Up to 2.1, we have few pseudo random number generators. We have an interface RandomGenerator implemented by three classes:
> - JDKRandomGenerator that extends the JDK provided generator
> - AbstractRandomGenerator as a helper for users generators
> - BitStreamGenerator which in turn is extended only by MersenneTwister
>
> The JDK provided generator is a simple one that can be used only for very simple needs. The Mersenne Twister is on the other hand a fast generator with good properties well suited for Monte-Carlo simulation. It is equidistributed for generating vectors up to dimension 623 and has a huge period: 2^19937 - 1.
>
> Since Mersenne-Twister inception in 1997, some new generators have been created, retaining the good properties of Mersenne twister but removing some of its (few) drawbacks. The main one is that if initialized with a bits pool containing lots of zeroes, the pool will take a very long time time to stabilize with a roughly balanced number of zeros and ones.
>
> I would like to add such generators (well, I already did but can withdraw my commit). The ones I want to add are the WELL generators (Well Equidistributed Long period Linear) created by François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. They are described in their 2006 paper: Improved
> * Long-Period Generators Based on Linear Recurrences Modulo 2, ransactions on Mathematical Software, 32, 1 (2006) which is available at <http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf>.
>
> The paper describes several generators from one family. The most interesting ones are WELL1024 (a small one), WELL19937 and WELL44497 (two large ones with huge periods that are both Mersenne primes). The two large ones exist in a basic version that is not Maximally Equidistributed (i.e. equidistributed in all possible dimensions and number of bits up to some threshold) and also in an extended version that add some tempering to the output to get an ME generator.
>
> The paper does not put any restriction on the algorithm. The reference implementation in C on the other hand is limited to non-commercial use only. In my implementation, I did not refer to this C implementation at all. I only compiled it on my personal computer to generate reference values and copied the values in the test cases. So this implementation is completely independent (and in fact the code is rather different since I used an abstract class for the global algorithm and one derived class for each specific generator. I rely on the JVM optimizer to do all the inlining and constants simplification that they did manually in the reference implementation.
>
> Another difference with the reference implementation is that our BitStreamGenerator abstract class requires the generator to provide only bits blocks and not double, and it computes the double by itself. This computation uses the full IEEE754 width for the doubles, i.e. it uses 53 bits. The generators are guaranteed to be equidistributed on large vectors of 32 bits (up to dimension 1390 for the 44497 version, since 1390 * 32 < 44497). I guess this should work well also using chunks of 53 bits (but of course with smaller dimensions), but it is not mathematically proven.
>
> So what do other think about this ?
> Should this really be included (in which case I will open a JIRA issue and resolve it immediately) or should it be removed from commons-math ? Is the independent implementation and the off-line reference data generation sufficient to say this code is unencumbered ? Is the 53 bits use for double generation a good thing or should we reduce it to 32 bits to stay within the proven behavior ?
>
> Once again, sorry for this too late discussion, I really should have started it before.
> Luc
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] new pseudo random number generators
Posted by Gilles Sadowski <gi...@harfang.homelinux.org>.
> [...]
> So what do other think about this ?
> Should this really be included (in which case I will open a JIRA issue and resolve it immediately)
+1
Gilles
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org