smile.math.random.UniversalGenerator Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of openchemlib Show documentation
Open Source Chemistry Library
There is a newer version: 2024.11.2
/*******************************************************************************
 * Copyright (c) 2010 Haifeng Li
 *   
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *  
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *******************************************************************************/

package smile.math.random;

/**
 * The so called "Universal Generator" based on multiplicative congruential
 * method, which originally appeared in "Toward a Universal Random Number
 * Generator" by Marsaglia, Zaman and Tsang. It was later modified by F. James
 * in "A Review of Pseudo-random Number Generators". It passes ALL of the tests
 * for random number generators and has a period of 2¹⁴⁴. It is
 * completely portable (gives bit identical results on all machines with at
 * least 24-bit mantissas in the floating point representation).
 * 
 * @author Haifeng Li
 */
public class UniversalGenerator implements RandomNumberGenerator {

    /**
     * Default seed.  DEFAULT_RANDOM_SEED=54217137
     */
    private static final int DEFAULT_RANDOM_SEED = 54217137;
    /**
     * The 46,009,220nd prime number,
     * The largest prime less than 9*10⁸.  Used as a modulus
     * because this version of random() needs a seed between 0
     * and 9*10⁸ and BIG_PRIME isn't commensurate
     * with any regular period.
     * BIG_PRIME = 899999963
     */
    static final int BIG_PRIME = 899999963;
    private double c, cd, cm, u[];
    private int i97, j97;

    /**
     * Initialize Random with default seed.
     */
    public UniversalGenerator() {
        setSeed(DEFAULT_RANDOM_SEED);
    }

    /**
     * Initialize Random with a specified integer seed
     */
    public UniversalGenerator(int seed) {
        setSeed(seed);
    }

    /**
     * Initialize Random with a specified long seed
     */
    public UniversalGenerator(long seed) {
        setSeed(seed);
    }

    @Override
    public void setSeed(long seed) {
        u = new double[97];

        int ijkl = Math.abs((int) (seed % BIG_PRIME));
        int ij = ijkl / 30082;
        int kl = ijkl % 30082;

        // Handle the seed range errors
        // First random number seed must be between 0 and 31328
        // Second seed must have a value between 0 and 30081
        if (ij < 0 || ij > 31328 || kl < 0 || kl > 30081) {
            ij = ij % 31329;
            kl = kl % 30082;
        }

        int i = ((ij / 177) % 177) + 2;
        int j = (ij % 177) + 2;
        int k = ((kl / 169) % 178) + 1;
        int l = kl % 169;

        int m;
        double s, t;
        for (int ii = 0; ii < 97; ii++) {
            s = 0.0;
            t = 0.5;
            for (int jj = 0; jj < 24; jj++) {
                m = (((i * j) % 179) * k) % 179;
                i = j;
                j = k;
                k = m;
                l = (53 * l + 1) % 169;
                if (((l * m) % 64) >= 32) {
                    s += t;
                }
                t *= 0.5;
            }
            u[ii] = s;
        }

        c = 362436.0 / 16777216.0;
        cd = 7654321.0 / 16777216.0;
        cm = 16777213.0 / 16777216.0;
        i97 = 96;
        j97 = 32;
    }

    @Override
    public double nextDouble() {
        double uni;

        uni = u[i97] - u[j97];
        if (uni < 0.0) {
            uni += 1.0;
        }

        u[i97] = uni;
        if (--i97 < 0) {
            i97 = 96;
        }

        if (--j97 < 0) {
            j97 = 96;
        }

        c -= cd;
        if (c < 0.0) {
            c += cm;
        }

        uni -= c;
        if (uni < 0.0) {
            uni++;
        }

        return uni;
    }
    
    @Override
    public void nextDoubles(double[] d) {
        int n = d.length;

        double uni;

        for (int i = 0; i < n; i++) {
            uni = u[i97] - u[j97];
            if (uni < 0.0) {
                uni += 1.0;
            }
            u[i97] = uni;
            if (--i97 < 0) {
                i97 = 96;
            }
            if (--j97 < 0) {
                j97 = 96;
            }
            c -= cd;
            if (c < 0.0) {
                c += cm;
            }
            uni -= c;
            if (uni < 0.0) {
                uni += 1.0;
            }
            d[i] = uni;
        }
    }

    @Override
    public int next(int numbits) {
        return nextInt() >>> (32 - numbits);
    }
    
    @Override
    public int nextInt() {
        return (int) Math.floor(Integer.MAX_VALUE * (2 * nextDouble() - 1.0));        
    }
    
    @Override
    public int nextInt(int n) {
        if (n <= 0) {
            throw new IllegalArgumentException("n must be positive");
        }
        
        return (int) (nextDouble() * n);
    }
    
    @Override
    public long nextLong() {
        return (long) Math.floor(Long.MAX_VALUE * (2 * nextDouble() - 1.0));        
    }
}