public static double lcgrand(int stream) {
        long zi, lowprd, hi31;

        zi = zrng[stream];
        lowprd = (zi & 65535) * MULT1;
        hi31 = (zi >> 16) * MULT1 + (lowprd >> 16);
        zi = ((lowprd & 65535) - MODLUS)
                + ((hi31 & 32767) << 16) + (hi31 >> 15);
        if (zi < 0) {
            zi += MODLUS;
        }
        lowprd = (zi & 65535) * MULT2;
        hi31 = (zi >> 16) * MULT2 + (lowprd >> 16);
        zi = ((lowprd & 65535) - MODLUS)
                + ((hi31 & 32767) << 16) + (hi31 >> 15);
        if (zi < 0) {
            zi += MODLUS;
        }
        zrng[stream] = zi;
        return (zi >> 7 | 1) / 16777216.0;
    }