As far as I am aware, there is no production quality native PL/I compiler
available for Linux. There is VisualAge PL/I for Windows, which IBM makes
available through its Scholars Program to those who qualify or which may be
purchased (at significant cost) as part of the Rational Developer for
System Z product.

Not  that I want to  plug a competing product, but Micro Focus Open PL/I
runs on Linux.  (http://www.microfocus.com/products/studio/open-pli.aspx)
What's your definition of "production quality?"  I haven't tried this,
but it at least sounds good.

Of course Micro Focus invites this anonymity, since their marketing is
probably second only to IBM's in its badness: no pricing, no demo, no
advertising or promotion, etc.

If you want something to try *now*, look at Iron Spring
(http://www.iron-spring.com/) - possibly not yet "production quality"
but getting there.  If it doesn't do what you want, just ask.

It produced identical results for me when I tried it in C using this
code:

typedef unsigned int Uint32;
Uint32 sign(Uint32 value) {return value>>31;}
Uint32 index, carry, data[4691];
Uint32 mwc4691() {
        index = (index < 4690) ? index + 1 : 0;
        Uint32 x = data[index];
        Uint32 t = (x << 13) + carry + x;
        if (sign((x<<13)+carry) == sign(x))
                carry = sign(x) + (x>>19);
        else carry = 1 - sign(t) + (x>>19);
        data[index] = t;
        return t;

From:
http://docs.sun.com/source/819-0492/4_f95.html

"4.5 Unsigned Integers
The Fortran 95 compiler accepts a new data type, UNSIGNED, as an
extension to the language. Four KIND parameter values are accepted with
UNSIGNED: 1, 2, 4, and 8, corresponding to 1-, 2-, 4-, and 8-byte
unsigned integers, respectively.

The form of an unsigned integer constant is a digit-string followed by
the upper or lower case letter U, optionally followed by an underscore
and kind parameter. The following examples show the maximum values for
unsigned integer constants:

                 255u_1
                 65535u_2
                 4294967295U_4
                 18446744073709551615U_8

Expressed without a kind parameter (12345U), the default is the same as
for default integer. This is U_4 but can be changed by the -xtypemap
option, which will change the kind type for default unsigned integers.

Declare an unsigned integer variable or array with the UNSIGNED type
specifier:

Aside from those qualifiers, the results are independent of whether
the numbers are considered to be signed or unsigned.  I probably
should have mentioned #4, but since I believe Fortran treats all right-
shifts as unsigned, it should not matter in Fortran.

Some C code which produces identical results and illustrates that:

typedef unsigned int Uint32;
typedef signed int Sint32;
Sint32 index, carry, data[4691];
Sint32 rshift(Sint32 value, int bits) {return Uint32(value)>>bits;}
Sint32 sign(Sint32 value) {return rshift(value,31);}
Sint32 mwc4691() {
   index = (index < 4690) ? index + 1 : 0;
   Sint32 x = data[index];
   Sint32 t = (x << 13) + carry + x;
   if (sign((x<<13)+carry) == sign(x)) carry = sign(x) + rshift(x,19);
   else carry = 1 - sign(t) + rshift(x,19);
   data[index] = t;
   return t;

(snip)

(Many Fortran rules haven't been changed for the same reason.)

There is no guarantee that twos complement overflow gives
the low bits of the full result, but it usually does.
(Fortran allows for almost anything to happen.)

(Hopefully ishft(1,-31) is a compile-time constant.  Otherwise,
substitute the appropriate constant.)

MWC () is supposed to implement multiplication of a huge number by
8193, with the huge number split into 32 bit items.
We could implement this by calculating (Q[i] * 8193 + carry) with 64
bit arithmetic, storing 32 bits, moving the upper 32 bits into carry.
In 32 bit arithmetic, we calculate the 32 bits to be stored as (Q[i] +
(Q[i]<<13) + carry), that part is fine.
The new carry should be (Q[i] >> 19), plus the carry in adding (Q[i] +
(Q[i]<<13) + carry) in 32 bit arithmetic.

The carry is found by adding all three items, and checking whether the
sum is less than Q [i].
That is incorrect when (Q[i]<<13) + carry produces a carry in 32 bit
arithmetic.
And that happens when we start with carry = 0x800 and Q [i] =
0xffffffff.

if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then
   if (t < x) then
     tLTx = 1
   else
     tLTx = 0
   endif
elseif (x < 0) then
   tLTx = 1
elseif (t < 0) then
   tLTx = 0
endif
c = tLTx + SHIFTR(x,19)

The following Fortran, modulo 2**32, gives the right answers:
The C code has been left in as comments.  The before and after
comments removed, but one could add them back in from the original
post.

As mentioned previously, this assumes 32 bit, twos complement, and
that on overflow the proper low order bits are returned.

I haven't done any timing, or time comparisons to C.

      integer function MWC()
! when did zero argument functions get added to Fortran?
      integer q
      common q(0:4690)
      integer t,x,c,j
      save c
      data j/4691/
! unsigned long t,x,i; static c=0,j=4691;
      j=j+1
      if(j>4690) j=0
!  j=(j<4690)? j+1:0;
      x=q(j)
!  x=Q[j];
      t=ishft(x,13)+c+x
      c=ishft(x,-19)
      if(ieor(t,ishft(1,31)).lt.ieor(x,ishft(1,31))) c=c+1
!  t=(x<<13)+c+x; c=(t<x)+(x>>19);
      q(j)=t
      MWC=t
      return
      end
!  return (Q[j]=t);

! int main()
      program main
      integer xs,xcng
      data xs,xcng/521288629,362436069/
! static unsigned long xs=521288629,xcng=362436069,Q[4691];
      integer i,x
      integer q
      common q(0:4690)
! {unsigned long i,x;
      do i=0,4690
         xcng=69069*xcng+123
! #define CNG ( xcng=69069*xcng+123 )
         xs=ieor(xs,ishft(xs, 13))
         xs=ieor(xs,ishft(xs,-17))
         xs=ieor(xs,ishft(xs,  5))
! #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
         q(i)=xcng+xs
      enddo
!  for(i=0;i<4691;i++) Q[i]=CNG+XS;
      do i=1,1000000000
         x=MWC()
      enddo
!  for(i=0;i<1000000000;i++) x=MWC();
      print *," MWC result=3740121002 (or -554846295) ?",x
!  printf(" MWC result=3740121002 ?\n%22u\n",x);
!  for(i=0;i<1000000000;i++) x=KISS;
      do i=1,1000000000
         xcng=69069*xcng+123
! #define CNG ( xcng=69069*xcng+123 )
         xs=ieor(xs,ishft(xs, 13))
         xs=ieor(xs,ishft(xs,-17))
         xs=ieor(xs,ishft(xs,  5))
! #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
         x=MWC()+xcng+xs
      enddo

! #define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
      print *,"KISS result=2224631993 (or -2070335303) ?",x;
!  printf("KISS result=2224631993 ?\n%22u\n",x);
! }
      end

new x = a*x+c mod b,  new c = integerpartof((a*x+c)/b)

and the new c can never reach, or exceed, a.

Formally, we may consider the set Z of a*b^4169 "vectors" z:

  z=[c;x_1,x_2,...,x_{4169}],  0 <= c < a and 0 <= x_i < 2^32

and a function f(z):  with t=a*x+c,

f([c;x_1,x_2,...,x_{4169}])=
   [trunc(t/b);x_2,x_3,...,x_{4169},t mod b],

If k is the order of b mod p, the directed graph:z -> f(z)
has (p-1)/k loops of size k and two "loops" of size 1:

f(z)=z for z=[0;0,0,...,0] and z=[b-1,a-1,a-1,...,a-1].

The concatenated sequence made by taking the rightmost x
from each z forms, in reverse order,
the base-b expansion of j/p for some 0 <= j <= p.

When z=[0;0,0,...,0]
we get 0/p=.0000000...
When z=[a-1;d,d,d,...,d] with d=b-1,
we get p/p=.dddd...

In our case, we have four possible loops,
two of size (p-1)/2, two of size 1,
for a total of  2(p-1)/2 + 2 = p+1 = ab^4169,
related by the mappings
c <-> a-1-c, x <-> b-1-x, j/p <-> (p-j)/p.

The function f has an inverse: with s=c*b+x_{4169},

f^{-1}(z)=[s mod a;trunc(s/a),x_1,x_2,...,x_{4168}],

so that with g(z)=f^{-1}(z),
the trailing element in each of the vectors
 g(z),g(g(z)),g(g(g(z))),...
forms, in normal order, the base-b expansion of some j/p.
But arithmetic in g(z) is not well-suited for computer
implementation, while that of f(z) is eminently so, the
reason I created it---the method, not the arithmetic.

By finding the prime p=ab^n-1 with
b=2^32, a=8193=2^13+1 and n=4169,
I as able to ensure that the necessary arithmetic
could be simply---and quickly--carried out with
the base b=2^32 arithmetic available to most of us.

In my original post I neglected to point out that the
initial choice of c should satisfy 0 <= c < a=8193
and that, choosing
c=0 and all the Q's zero, or
c=a-1 and all the Q's b-1
would makes the MWC generator have period 1
and the resulting KISS period a mere 2^64-2^32
rather than one exceeding 10^40000.

Please forgive me.

the problem is not that the carry could match or exceed the multiplier
8193.
The problem is that the carry can be as large as 8192.

In your code you write

  t=(x<<13)+c+x; c=(t<x)+(x>>19);

To make it a bit more obvious what happens, I can change this to the
100% equivalent

  t = (x << 13) + c;
  c = x >> 19;
  t += x; if (t < x) ++c;

The last line is the usual trick to add numbers and check for carry: t
will be less than x exactly if the addition t += x wrapped around and
produced a carry.

The problem is that (x << 13) + c can also produce a carry, and that
goes unnoticed. If c = 8192, and the lowest 19 bits in x are all set,
then (x << 13) has the highest 19 bits set and 13 bits zero. Adding c
"overflows" and gives a result of 0, and adding x leaves x. The
comparison (t < x) is false and produces zero, even though (x<<13)+c+x
should have produced a carry.

This is very, very rare. It doesn't actually happen in your original
source code with a total of 2 billion random numbers. If you produce
four billion random numbers, then it happens at i = 3596309491 and
again at i = 3931531774.

I am more interested in 64 bit performance, so I just made c, t, and x
64 bit numbers and wrote

  uint64_t x = Q [i];
  uint64_t t = (x << 13) + x + c;
  c = (t >> 32);

   MWC result=3740121002 (or -554846295) ?  -554846294
  KISS result=2224631993 (or -2070335303) ? -2070335303

integer :: u = 3740121002
write(*,*) u

"christian.bau" <christian....@cbau.wanadoo.co.uk> ha scritto nel messaggio
news:9309cd00-acf9-4908-9f88-d8e2ffad4e70@x21g2000yqa.googlegroups.com...

---------------------------------------------
[mesA: D$ 0
mesB: D$ 0

StatC: D$ 0
StatJ: D$ 4691
xs: D$ 521288629
xcng: D$ 362436069
Q: D$ 0 #2346
secondQ: D$ 0 #2346 ; 2346*2=4692
; 0..4690

IMWCIIrisultatoI3740121002I: B$ "MWC risultato 3740121002?", 0
IKISSIrisultatoI2224631993I: B$ "KISS risultato 2224631993?", 0
msg: B$ "On each window press ok and wait for the next window", 0 , 0

]
[ArrForAll: D$ ? #2048
]

mwc:
push  esi |push   edi
mov   eax, D$StatJ  |cmp   eax,  4690 |jae  @1 |inc   eax |jmp   @2
@1:  xor   eax,  eax
@2:  mov  D$StatJ ,  eax
mov   edx, D$Q+eax*4  |mov   ecx, D$StatC  |mov   edi,  edx |xor   esi,  esi
shl   edx,  13 |add   edx,  ecx |adc   esi,  0
add   edx,  edi |adc   esi,  0
shr   edi,  19 |lea   ecx, D$edi+esi
mov  D$StatC ,  ecx
mov  D$Q+eax*4 ,  edx
mov   eax,  edx
pop    edi |pop   esi
ret

CNGmXS:
mov   eax, D$xcng  |mov   edx,  69069 |mul edx |add   eax,  123 |mov  D$xcng ,
eax
mov   eax, D$xs  |mov   ecx, D$xs  |shl   eax,  13 |xor   ecx,  eax
mov   eax,  ecx |shr   eax,  17 |xor   ecx,  eax
mov   eax,  ecx |shl   eax,  5 |xor   ecx,  eax
mov  D$xs ,  ecx |mov   eax, D$xcng  |add   eax,  ecx
ret

KISS:
push  esi |push   edi
mov   eax, D$StatJ  |cmp   eax,  4690 |jae  @1 |inc   eax |jmp   @2
@1:  xor   eax,  eax
@2:  mov  D$StatJ ,  eax
mov   edx, D$Q+eax*4  |mov   ecx, D$StatC  |mov   edi,  edx |xor   esi,  esi
shl   edx,  13 |add   edx,  ecx |adc   esi,  0
add   edx,  edi |adc   esi,  0
shr   edi,  19 |lea   ecx, D$edi+esi
mov  D$StatC ,  ecx
mov  D$Q+eax*4 ,  edx
pop    edi |pop   esi
push  edx
mov   eax, D$xcng  |mov   edx,  69069 |mul edx |add   eax,  123 |mov  D$xcng ,
eax
mov   eax, D$xs  |mov   ecx, D$xs  |shl   eax,  13 |xor   ecx,  eax
mov   eax,  ecx |shr   eax,  17 |xor   ecx,  eax
mov   eax,  ecx |shl   eax,  5 |xor   ecx,  eax
mov  D$xs ,  ecx |mov   eax,  ecx |mov   eax, D$xcng  |add   eax,  ecx
pop   edx
add   eax,  edx
ret

Main:
pushad
;iint3
push  &MB_OK |push   msg |push   msg |push  0 |call   'user32.MessageBoxA'
mov   esi,  0
@1:  call  CNGmXS |mov  D$Q+esi*4 ,  eax |inc   esi |cmp   esi,  4691 |jb   @1
;iint3
mov   esi,  0
@2:  call  mwc |inc   esi |cmp   esi,  1000000000 |jb   @2
mov  D$mesA ,  eax
mov   edx,  esp |mov   D$esp ,   "%u"  |push   eax |push   edx |push  ArrForAll
call   'user32.wsprintfA' |add   esp,  12
push  &MB_OK |push   IMWCIIrisultatoI3740121002I  |push   ArrForAll |push  0
call   'user32.MessageBoxA'

;iint3
mov   esi,  0
@3:  call  KISS |inc   esi |cmp   esi,  1000000000 |jb   @3
mov  D$mesB ,  eax
mov   edx,  esp |mov   D$esp ,   "%u"  |push   eax |push   edx |push  ArrForAll
call   'user32.wsprintfA' |add   esp,  12
push  &MB_OK |push   IKISSIrisultatoI2224631993I  |push   ArrForAll |push  0
call   'user32.MessageBoxA'
;iint3
popad
mov   eax,  0
ret
------------------End RosAsm

-------bigin the language i use
PREPARSE MACRO2D

section DATA

%define  MsgBoxA             'user32.MessageBoxA'
%define  wsprintfA           'user32.wsprintfA'

mesA    dd           0
mesB    dd           0

StatC   dd           0
StatJ   dd        4691
xs        dd   521288629
xcng      dd   362436069
Q:         times  2346    dd  0
secondQ    times  2346    dd  0   ; 2346*2=4692
; 0..4690

"MWC  risultato 3740121002?"   db  "MWC risultato 3740121002?" , 0
"KISS risultato 2224631993?"   db  "KISS risultato 2224631993?", 0
msg  db "On each window press ok and wait for the next window", 0, 0

section       BSS
ArrForAll:    resd   2048
section TEXT

mwc:
<i,j
    a=*StatJ|a<4690!#.1|++a|#.2
.1:            a^=a
.2: *StatJ=a
    r=*Q+a*4|c=*StatC|j=r|i^=i
    r<<=13  |r+=c   |i++=0
    r+=j            |i++=0
    j>>=19  |c=&*j+i
    *StatC=c
    *Q+a*4=r
    a=r

CNGmXS:
    a=*xcng|r=69069|mul  r|a+=123|*xcng=a
    a=*xs|c=*xs|a<<=13|c^=a
    a=c        |a>>=17|c^=a
    a=c        |a<<=5 |c^=a
    *xs=c|a=*xcng|a+=c
ret

KISS:
<i,j
    a=*StatJ|a<4690!#.1|++a|#.2
.1:            a^=a
.2: *StatJ=a
    r=*Q+a*4|c=*StatC|j=r|i^=i
    r<<=13  |r+=c   |i++=0
    r+=j            |i++=0
    j>>=19  |c=&*j+i
    *StatC=c
    *Q+a*4=r

Main:
pushad
    ;iint3
    MsgBoxA(0, msg, msg, &MB_OK )
    i=0
.1: CNGmXS()|*Q+i*4=a|++i|i<4691#.1
    ;iint3
    i=0
.2: mwc()   |++i|i< 1000000000#.2
    *mesA=a
    r=s|D*s="%u"|wsprintfA<(ArrForAll,  r, a)
    MsgBoxA(0, ArrForAll, "MWC  risultato 3740121002?" , &MB_OK )

The entire MWC() routine should have that
(t<x) to (t<=x) correction:

unsigned long MWC(void)
{ unsigned long t,x,i; static c=0,j=4691;
    j=(j<4690)? j+1:0;
    x=Q[j];
    t=(x<<13)+c+x; c=(t<=x)+(x>>19);
    return (Q[j]=t);

Keeping that (t<x) may be viewed as, rather than making a
circular transition through the 8193*2^133407-1 base-b digits
of the expansion of some j/p, we jump---after an average of
b=2^32 steps---to another point in the immense circle.

It could be that the period is increased rather than
decreased, but it remains a curiosity. Perhaps light could
be shed with primes p=ab^n-1 for which the actual distribution
of such jumps, averaging every b steps, could be examined.

Or perhaps it could remain a mystery
and be further fodder for those who tend to
equate lack of understanding to "true" randomness.