[Bug libstdc++/51795] New: linear_congruential_engine doesn't work correctly

[Bug libstdc++/51795] New: linear_congruential_engine doesn't work correctly

[spoon.reloaded at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

             Bug #: 51795
           Summary: linear_congruential_engine doesn't work correctly
    Classification: Unclassified
           Product: gcc
           Version: 4.6.1
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libstdc++
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: spoon.reloaded@gmail.com


I am trying to use a simple linear congruential engine from the random library
in C++11, and I am not getting what I expect.

#include <random>
#include <iostream>

int main() {
  std::linear_congruential_engine<uint64_t, 1103515245, 12345, 2147483648>
foo(1103527590);
  std::cout << foo() << std::endl; // prints 17294696140058308839

  return 0;
}

This is obviously incorrect, as the result is far bigger than the modulus
(2147483648). I think the correct result should be (a * seed + c) mod m =
(1103515245 * 1103527590 + 12345) % 2147483648 = 377401575. Am I missing
something?

I am using the gcc and libstdc++ from Ubuntu oneiric, which is version 4.6.1. I
am using a 32-bit machine.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[3dw4rd at verizon dot net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #1 from Ed Smith-Rowland <3dw4rd at verizon dot net> 2012-01-27 04:48:39 UTC ---
I random.tcc in _Mod
__x = 1103527590
__m = 2147483648
__a = 1103515245
__c = 12345
          static const _Tp __q = __m / __a; // 1
          static const _Tp __r = __m % __a; // 1043968403

          _Tp __t1 = __a * (__x % __q);     // 0
          _Tp __t2 = __r * (__x / __q);     // 1152047935798738770
          if (__t1 >= __t2)
        __x = __t1 - __t2;              // 
          else
        __x = __m - __t2 + __t1;        // 17294696140058296494 !!!

      if (__c != 0)
        {
          const _Tp __d = __m - __x;        // 1152047935798738770
          if (__d > __c)
        __x += __c;                     // 17294696140058308839
          else
        __x = __c - __d;
        }

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |paolo.carlini at oracle dot
                   |                            |com

--- Comment #2 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 10:48:34 UTC ---
Thus Ed, is the fix already obvious to you? Or should I look into this? Thanks.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |ASSIGNED
   Last reconfirmed|                            |2012-01-27
         AssignedTo|unassigned at gcc dot       |paolo.carlini at oracle dot
                   |gnu.org                     |com
     Ever Confirmed|0                           |1

--- Comment #3 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 11:02:46 UTC ---
Thus looks like we have a very very long standing bug in the implementation of
Schrage's algorithm, which nobody noticed so far. Let's fix it.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #4 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 11:32:10 UTC ---
The problem is that isn't true that __r < __q, thus one of the preconditions
for the algorithm (as explained in eg, Numerical Recipes) is violated.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #5 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 11:41:28 UTC ---
Roughly, we have to fix either __calc or the generator itself to work fine also
when isn't true that __a << __m (which is the "normal" situation)

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|paolo.carlini at oracle dot |marc.glisse at normalesup
                   |com                         |dot org

--- Comment #6 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 11:45:29 UTC ---
Marc, can you see a "quick and dirty" (ie, suited for 4.7 too) way to fix this?

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #7 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 12:35:28 UTC ---
As a matter of fact, I don't see why here (and elsewhere, eg generate) we have
to use the Schrage trick to avoid integer overflow: we know that the integers
are unsigned and I don't see anything about that in the Standard.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #8 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 12:47:10 UTC ---
Eg:

Index: random.h
===================================================================
--- random.h    (revision 183615)
+++ random.h    (working copy)
@@ -263,7 +263,7 @@
       result_type
       operator()()
       {
-    _M_x = __detail::__mod<_UIntType, __m, __a, __c>(_M_x);
+    _M_x = (__a * _M_x + __c) % __m;
     return _M_x;
       }

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[marc.glisse at normalesup dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #9 from Marc Glisse <marc.glisse at normalesup dot org> 2012-01-27 12:53:38 UTC ---
(In reply to comment #6)
> Marc, can you see a "quick and dirty" (ie, suited for 4.7 too) way to fix this?

Note that I am not a specialist.
1) add an assertion "not implemented yet, try a smaller constant"
2) if you know how to multiply by small a, you can decompose large a into
a1*a2+a3 and perform the operations separately (variants might make sense,
including using powers of 2 and also splitting x)
3) do the double precision operation. longlong.h contains generic code for that
(but beware UDIV_NEEDS_NORMALIZATION).

About comment #7, if that's how you interpret it, it may be worth opening a DR.
I believe the formula has to be interpreted with mathematical, arbitrary size
integers (although of course you can optimize when m*m fits in _Tp).

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #10 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 13:01:23 UTC ---
Really, I don't think the intent here is involving arbitrary size integers.
This stuff must be quick, and all the integers are unsigned. I don't think we
have here something like the templates in time, really.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[marc.glisse at normalesup dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #11 from Marc Glisse <marc.glisse at normalesup dot org> 2012-01-27 13:12:05 UTC ---
(In reply to comment #10)
> Really, I don't think the intent here is involving arbitrary size integers.
> This stuff must be quick, and all the integers are unsigned. I don't think we
> have here something like the templates in time, really.

But if we round first modulo 2^64 then modulo m, it isn't the linear
congruential algorithm it is supposed to be. I think I'd prefer for the
compiler to reject the code when used with values that would make the generator
too slow.

Let me ask on the reflector, to be sure.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #12 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 13:17:22 UTC ---
I see. Before asking: 26.5/4 says that "all descriptions of calculations in
this subclause use mathematical real numbers". Thus should we use floats?!?

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #13 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 13:23:01 UTC ---
Note that besides the case of the linear_congruential (where we could
definitely static_assert, I don't consider such uses really important in
practice, given also that generators much better from the statistical point of
view are highly recommended (for example the last Numeric Recipes recommends
avoiding linear congruential and doesn't provide anymore code for it), we have
also to resolve the problem with the implementation of generate, where an
algorithm is described in detail involving multiplications and currently we
also use __mod for unsafe arguments.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[marc.glisse at normalesup dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #14 from Marc Glisse <marc.glisse at normalesup dot org> 2012-01-27 13:35:29 UTC ---
(In reply to comment #12)
> I see. Before asking: 26.5/4 says that "all descriptions of calculations in
> this subclause use mathematical real numbers".

Ok, no need to ask then.

> Thus should we use floats?!?

Actually, using double is a common trick (that wouldn't be sufficient here):
compute (double)a*(double)x/(double)m and round it to an integer q; then you
can compute a*x-m*q modulo 2^64.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #15 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 13:42:27 UTC ---
I'm not sure to understand: is using floating point quantities the way to go or
not? I'm particularly worried by generate, as I said, much more than linear
congruential. Note we can also use long double.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #16 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 13:54:44 UTC ---
Uhm, in the case of generate, however, "each operation is to be carried out
modulo 2^32", thus I guess it's safe to not use __mod at all for the
multiplication.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[marc.glisse at normalesup dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #17 from Marc Glisse <marc.glisse at normalesup dot org> 2012-01-27 13:57:43 UTC ---
(In reply to comment #15)
> I'm not sure to understand: is using floating point quantities the way to go or
> not?
> Note we can also use long double.

I don't think there is any guarantee that long double has enough precision.

> I'm particularly worried by generate, as I said, much more than linear
> congruential.

Then it may be safer to implement a slow but exact fallback implementation of
_Mod.

Sorry, I can't look at it in more details right now...

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #18 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 14:04:28 UTC ---
Okay, thanks, but for generate we are 'lucky'. Thus, for 4.7 I'm going to add
the static_assert in linear_congruential and tweak generate to not use __mod at
all. Then we'll see.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[spoon.reloaded at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #19 from spoon.reloaded at gmail dot com 2012-01-27 21:21:40 UTC ---
Paulo, in response to your suggestion to simply do multiplication and modulo in
#7 and #8, I don't think that would work in general. The example I gave
happened to have m = a power of 2 (namely 2^31), and so the truncation that we
would get from integer overflow (whether by 2^32 if we use uint32_t or 2^64 if
we use uint64_t) does not affect the result. However, if we choose any other
number as a modulo (e.g. 2^31 - 1) and say we use uint32_t as the type, it will
not work:

(1103515245 * 1103527590 + 12345) % 2147483647 = 944465040

but in uint32_t arithmetic:

(1103515245 * 1103527590 + 12345) % (1 << 32) % 2147483647 = 377401576

This is why I think we still need something like the Schrage's algorithm

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[spoon.reloaded at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #20 from spoon.reloaded at gmail dot com 2012-01-27 21:39:11 UTC ---
Someone posted a thread here:
http://objectmix.com/java/312426-extending-schrage-multiplication.html
about ways to overcome the preconditions of Shrage's algorithm. Maybe that can
help

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #21 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 21:42:53 UTC ---
Paolo, not Paulo, by the way. Sure, you are right, I didn't personally write
these bits of random and for years I simply trusted the documentation which
doesn't point to any restriction. Too bad. Anyway, for 4.7.0 we are simply
going to static_assert the requirements. In general, as I already said, to be
fully honest I don't consider the issue particularly serious, because LCG is a
pretty poor generator anyway and sensible values of the constants (per the
typedefs too) are already supported well. Thus for 4.8.0 I will welcome
patches, or at some point I will get around working on the issue but I don't
have it in my TODO list for the immediate future.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #22 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 21:43:50 UTC ---
Thanks for the link!

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[spoon.reloaded at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #23 from spoon.reloaded at gmail dot com 2012-01-27 21:52:48 UTC ---
By the way, the Boost library also has an implementation of this. Are there any
problems with the license to look at it?

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #24 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-27 22:01:14 UTC ---
There aren't if afterwards you don't submit a patch to GCC ;) Seriously, I'm
pretty sure that there are many implementations around and Boost would be by
and large ok, but normally it's much better and straightforward from the legal
point of view if you look only at GPL code, eg, the GNU Scientific Library
would be perfect.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[marc.glisse at normalesup dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #25 from Marc Glisse <marc.glisse at normalesup dot org> 2012-01-27 22:32:17 UTC ---
On the other hand, you can test the boost implementation without looking at the
code and notice the assertion failure if you ask for uint32_t (with uint64_t,
there is enough space to do the multiplication). So they are cleaner but don't
handle everything either.

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #26 from paolo at gcc dot gnu.org <paolo at gcc dot gnu.org> 2012-01-27 23:30:41 UTC ---
Author: paolo
Date: Fri Jan 27 23:30:28 2012
New Revision: 183655

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=183655
Log:
2012-01-27  Paolo Carlini  <paolo.carlini@oracle.com>

    PR libstdc++/51795
    * include/bits/random.h (linear_congruential_generator): Add
    static_assert preventing instantiation for values of 'a' and 'm'
    currently handled incorrectly but _Mod::__calc.
    * include/bits/random.tcc (seed_seq::generate): Avoid unsafe
    uses of _Mod::__calc.

Modified:
    trunk/libstdc++-v3/ChangeLog
    trunk/libstdc++-v3/include/bits/random.h
    trunk/libstdc++-v3/include/bits/random.tcc

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #27 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-28 10:51:12 UTC ---
It seems to me that something still safe to do now for 4.7.0 is special casing
the case x == m - 1 to x' = m - a, because, if I understand correctly what I'm
reading around, in that case too Schrage in general is not correct, ie, for the
constants proposed by Park & Miller it works fine only by chance. Can you
people confirm?

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #28 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-01-30 15:52:00 UTC ---
As a matter of fact, I'm not able to prove that things can go wrong with the
normal Schrage when x == m - 1, at least given our other conds (eg, a < m). I
guess better not fiddling with this, unless somebody has a testcase...

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #29 from paolo at gcc dot gnu.org <paolo at gcc dot gnu.org> 2012-02-01 11:10:35 UTC ---
Author: paolo
Date: Wed Feb  1 11:10:30 2012
New Revision: 183795

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=183795
Log:
2012-02-01  Paolo Carlini  <paolo.carlini@oracle.com>

    PR libstdc++/51795
    * include/bits/random.h (linear_congruential_generator): Add
    static_assert preventing instantiation for values of 'a' and 'm'
    currently handled incorrectly but _Mod::__calc.
    * include/bits/random.tcc (seed_seq::generate): Avoid unsafe
    uses of _Mod::__calc.

Modified:
    branches/gcc-4_6-branch/libstdc++-v3/ChangeLog
    branches/gcc-4_6-branch/libstdc++-v3/include/bits/random.h
    branches/gcc-4_6-branch/libstdc++-v3/include/bits/random.tcc

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

--- Comment #30 from paolo at gcc dot gnu.org <paolo at gcc dot gnu.org> 2012-04-29 23:36:15 UTC ---
Author: paolo
Date: Sun Apr 29 23:36:09 2012
New Revision: 186948

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=186948
Log:
2012-04-29  Marc Glisse  <marc.glisse@inria.fr>
        Paolo Carlini  <paolo.carlini@oracle.com>

    PR libstdc++/51795
    * include/bits/stl_algobase.h (__lg<>(_Size)): Remove.
    (__lg(int), __lg(unsigned), __lg(long), __lg(unsigned long),
    __lg(long long), __lg(unsigned long long)): Define constexpr.
    * include/bits/random.h (_Mod<>): Overcome Schrage's algorithm
    limitations.
    (__mod): Adjust.
    (linear_congruential): Remove FIXME static_assert.
    * include/bits/random.tcc (_Mod<>): Adjust.
    * testsuite/26_numerics/random/linear_congruential_engine/operators/
    51795.cc: New.

Added:
   
trunk/libstdc++-v3/testsuite/26_numerics/random/linear_congruential_engine/operators/51795.cc
Modified:
    trunk/libstdc++-v3/ChangeLog
    trunk/libstdc++-v3/include/bits/random.h
    trunk/libstdc++-v3/include/bits/random.tcc
    trunk/libstdc++-v3/include/bits/stl_algobase.h

[Bug libstdc++/51795] linear_congruential_engine doesn't work correctly

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51795

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|                            |FIXED
   Target Milestone|---                         |4.8.0

--- Comment #31 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-04-29 23:38:45 UTC ---
Completely fixed in mainline.