Random(...,1,2^80) test fails on SPARC Solaris

[dimpase]

Observed behaviour

testing: /datapool/dima/Sage/gap/tst/testinstall/intarith.tst
########> Diff in /datapool/dima/Sage/gap/tst/testinstall/intarith.tst:948
# Input is:
Random(mysource, 1, 2^80);
# Expected output:
1071135285340982180054653
# But found:
2461224813218428180

(the other tests there pass)

Expected behaviour

No diff, obviously

Copy and paste GAP banner (to tell us about your setup)

 *********   GAP 4.8.8-5202-gb88ad0f of today
 *  GAP  *   https://www.gap-system.org
 *********   Architecture: sparc-sun-solaris2.11-default64
 Configuration:  gmp 6.1.2
 Loading the library and packages ...
 Packages:   AClib 1.2, Alnuth 3.1.0, AtlasRep 1.5.1, AutPGrp 1.8, CRISP 1.4.4, Cryst 4.1.13, CrystCat 1.1.6, CTblLib 1.2.2, 
             FactInt 1.6.1, FGA 1.3.1, GAPDoc 1.6.1, IRREDSOL 1.4, LAGUNA 3.8.0, Polenta 1.3.8, Polycyclic 2.11, PrimGrp 3.3.0, 
             RadiRoot 2.7, ResClasses 4.7.1, smallgrp 1.2, Sophus 1.23, SpinSym 1.5, TomLib 1.2.6, TransGrp 2.0.2, utils 0.53

Built with vanilla gcc 7.2.

---

Note that as #2274 is merged, this and #2194 are the only things to be fixed for testinstall pass on this platform.

#116 makes it hard to understand what's going on, precisely...

[ChrisJefferson]

Would it be possible to bisect this issue if it only appeared in the last couple of weeks, or has it been around for a long time (I made some recentish random changes).

[dimpase]

I don't recall seeing it on 4.9. I will bisect this.

[dimpase]

It's still there on 4.9.0, but it does work on 4.8.6. I will try finding the exact commit.

It's also not in 4.8.10, so one has to chase the 4.9 branch.

[dimpase]

Could you give me a reasonable starting point on the 4.9 branch to start bisecting from?

[fingolfin]

@dimpase One thing you can try is to start bisecting after the big build system work was settled, that makes it much less painful.
From a quick look, perhaps something like 3e23772 would be a good starting point -- assuming it builds for you, and assuming the error does not occur there.

[fingolfin]

BTW, is that machine a big endian machine by chance? That would give an important hint: it would suggest that some of the kernel random code is not sufficiently endian safe (e.g. the mersenne twister). It then probably wasn't in GAP 4.8 either, but we also changed which random sources we use to be more consistent, which could have exposed the problem.

This is just a wild guess, though.

[dimpase]

Yes, the endianness is big, as is the case with SPARCs.

[fingolfin]

OK. (Note that SPARC v9 and later is bi-endian, so given the information in this issue, I was not able to able to tell the endianess).

[fingolfin]

A quick look at the kernel code, specifically FuncRandomIntegerMT, suggests that my guess may be correct: that function does not seem to take endianess into account; even if the underlying mersenne twister code does, still FuncRandomIntegerMT writes data 4 bytes at a time, which could lead to differing results on a 64bit big endian machine. But this is still guessing, I did not try to understand this code properly.

[dimpase]

Just to be sure, here is the test:

$ python
Python 2.7.9 (default, Apr 27 2017, 13:50:01) [C] on sunos5
Type "help", "copyright", "credits" or "license" for more information.
>>> from sys import byteorder
>>> print(byteorder)
big

[fingolfin]

Excellent!
You could test if the mersenne twister code is involved or not by comparing with this snippet:

gap> state:=InitRandomMT("abc");;
gap> RandomIntegerMT(state, 80);
358386947634599872691919

If it gives identical results, we should look elsewhere; if not we have found the bug (or at least part of it).

[dimpase]

test if the mersenne twister code is involved or not

Both GAP 4.8.6 and 4.9.0 on SPARC Solaris 11 return

1539244741327137841491211282754109

I guess it's hidden on 4.8.6 only because tests for this are not present.

[fingolfin]

So one needs to look at FuncRandomIntegerMT in src/integer.c, which does this:

Obj FuncRandomIntegerMT(Obj self, Obj mtstr, Obj nrbits)
{
  Obj res;
  Int i, n, q, r, qoff, len;
  UInt4 *mt, rand;
  UInt4 *pt;
...
  n = INT_INTOBJ(nrbits);
...
  else {
     /* large int case */
     q = n / 32;
     r = n - q * 32;
     /* qoff = number of 32 bit words we need */
     qoff = q + (r==0 ? 0:1);
     /* len = number of limbs we need (limbs currently are either 32 or 64 bit wide) */
     len = (qoff*4 +  sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
     res = NewBag( T_INTPOS, len*sizeof(mp_limb_t) );
     pt = (UInt4*) ADDR_INT(res);
     mt = (UInt4*) CHARS_STRING(mtstr);
     for (i = 0; i < qoff; i++, pt++) {
       rand = (UInt4) nextrandMT_int32(mt);
       *pt = rand;  // <--- PROBLEM IS AROUND HERE ...
     }
     if (r != 0) {
       /* we generated too many random bits -- chop of the extra bits */
       pt = (UInt4*) ADDR_INT(res);
       pt[qoff-1] = pt[qoff-1] & ((UInt4)(-1) >> (32-r));  // <--- ... AND HERE
     }
     /* shrink bag if necessary */
     res = GMP_NORMALIZE(res);
     /* convert result if small int */
     res = GMP_REDUCE(res);
  }

  return res;
}

I marked the two problematic spots (roughly): this code creates a GMP bigint 32bits at a time. This ensures that we get identical results regardless of whether we are on a 32 or 64 bit system.

I am not sure whether we get identical results on 32 bit LE vs BE systems. Perhaps you could try it? i.e. can you build a 32bit SPARC versions of GAP and see what it outputs? If that already is broken, we need to change the byte order when writing in those two places.

But in addition, for 64bit BE, also the 32bit words need to be shuffled. To do this, in the first loop, one should call nextrandMT_int32 twice and swap the write order. It might be easier, though, to simply use an #ifdef SYS_IS_64_BIT to split the code into a 32bit version, and a 64bit version; and then change pt to be of type mp_limb_t; and then write 64bit at a time as follows:

   rand = (UInt4) nextrandMT_int32(mt);
  *pt = ((UInt8) nextrandMT_int32(mt) << 32) | rand;

(except in the last iteration, if the total number of 32 bit blocks is odd, then of course the high bits should be set to 0 instead.

And of course the final masking must be done carefully, too.

[dimpase]

Why does GAP roll its own randoms here, instead of relying on GMP? Is it because GMP does not do a good job here?

[ChrisJefferson]

Until a couple of months ago, gmp was an optional dependency. Also, we have to be careful using gmp functionality because we need to be careful about memory manager interactions.