|
| 1 | +/** |
| 2 | + * Copyright 2026 Alexander Herzog |
| 3 | + * |
| 4 | + * Licensed under the Apache License, Version 2.0 (the "License"); |
| 5 | + * you may not use this file except in compliance with the License. |
| 6 | + * You may obtain a copy of the License at |
| 7 | + * |
| 8 | + * http://www.apache.org/licenses/LICENSE-2.0 |
| 9 | + * |
| 10 | + * Unless required by applicable law or agreed to in writing, software |
| 11 | + * distributed under the License is distributed on an "AS IS" BASIS, |
| 12 | + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 13 | + * See the License for the specific language governing permissions and |
| 14 | + * limitations under the License. |
| 15 | + */ |
| 16 | +package mathtools.distribution.tools; |
| 17 | + |
| 18 | +/** |
| 19 | + * Minimal, compatible implementation of C's drand48/srand48 and erand48 in Java. |
| 20 | + * |
| 21 | + * Uses the POSIX/standard 48-bit LCG: |
| 22 | + * X_{n+1} = (a * X_n + c) mod 2^48 |
| 23 | + * with a = 0x5DEECE66DL, c = 0xB. |
| 24 | + * |
| 25 | + * - srand48(seedval) sets the internal 48-bit state to ((seedval << 16) + 0x330E) & mask |
| 26 | + * - drand48() advances state and returns state / 2^48 as a double in [0,1) |
| 27 | + * - erand48(xsubi) advances the state stored in the provided 3-element short[] and returns the double |
| 28 | + * |
| 29 | + * Notes about Java types: C uses unsigned short for xsubi; Java shorts are signed, |
| 30 | + * so we mask with 0xFFFF when converting to/from the short[]. |
| 31 | + */ |
| 32 | +public final class Drand48Mix implements Drand48Interface { |
| 33 | + /** LCG parameter (48-bit) MULTIPLIER */ |
| 34 | + private static final long MULTIPLIER = 0x5DEECE66DL; |
| 35 | + /** LCG parameter (48-bit) ADDEND */ |
| 36 | + private static final long ADDEND = 0xBL; |
| 37 | + /** LCG parameter (48-bit) MASK48 */ |
| 38 | + private static final long MASK48 = (1L << 48) - 1; |
| 39 | + |
| 40 | + /** internal 48-bit state (only low 48 bits used) */ |
| 41 | + private long state; |
| 42 | + |
| 43 | + /** |
| 44 | + * Create and seed the generator with given seedval (like srand48). |
| 45 | + * @param seedval seed value (any long; only low bits used similarly to C) |
| 46 | + */ |
| 47 | + public Drand48Mix(final long seedval) { |
| 48 | + srand48(seedval); |
| 49 | + } |
| 50 | + |
| 51 | + /** |
| 52 | + * Create and seed the generator with given seedval (like srand48). |
| 53 | + */ |
| 54 | + public Drand48Mix() { |
| 55 | + srand48(0); |
| 56 | + } |
| 57 | + |
| 58 | + /** |
| 59 | + * 128-bittige Mix-Funktion |
| 60 | + * @param x Eingangswert |
| 61 | + * @return Ausgangswert nach Mix |
| 62 | + * @see L32X64Mix#lea32(int) |
| 63 | + */ |
| 64 | + public static long mixLong(final long x) { |
| 65 | + int low = (int) (x & 0xFFFFFFFFL); /* untere 32 Bit */ |
| 66 | + int high = (int) ((x >>> 32) & 0xFFFFFFFFL); /* obere 32 Bit */ |
| 67 | + |
| 68 | + /* Schritt 2: Mix-Funktion auf beide Teile anwenden */ |
| 69 | + int mixedLow = L32X64Mix.lea32(low); |
| 70 | + int mixedHigh = L32X64Mix.lea32(high); |
| 71 | + |
| 72 | + /* Schritt 3: Wieder zusammenführen zu einem long */ |
| 73 | + /* Beachte: mixedHigh ist jetzt 32 Bit, aber als int (signed). |
| 74 | + * Wir müssen sicherstellen, dass es als unsigned 32-Bit interpretiert wird. */ |
| 75 | + return ((long) mixedHigh << 32) | (mixedLow & 0xFFFFFFFFL); |
| 76 | + } |
| 77 | + |
| 78 | + /** |
| 79 | + * Seed the generator (like C's srand48). |
| 80 | + * state := ((seedval << 16) + 0x330E) & ((1<<48)-1) |
| 81 | + * @param seedval seed value |
| 82 | + */ |
| 83 | + @Override |
| 84 | + public void srand48(final long seedval) { |
| 85 | + state = ((seedval << 16) + 0x330EL) & MASK48; |
| 86 | + } |
| 87 | + |
| 88 | + /** |
| 89 | + * Advance the internal state and return a double in [0.0, 1.0). |
| 90 | + * Equivalent to C's drand48(). |
| 91 | + * @return uniformly distributed double in [0,1) |
| 92 | + */ |
| 93 | + @Override |
| 94 | + public double drand48() { |
| 95 | + state = (state * MULTIPLIER + ADDEND) & MASK48; |
| 96 | + return (state / (double) (1L << 48)); |
| 97 | + } |
| 98 | + |
| 99 | + /** |
| 100 | + * Return a 31-bit non-negative integer like C's lrand48 (optional helper). |
| 101 | + * @return int in [0, 2^31) |
| 102 | + */ |
| 103 | + @Override |
| 104 | + public int lrand48() { |
| 105 | + // advance state and return high-order 31 bits |
| 106 | + state = (state * MULTIPLIER + ADDEND) & MASK48; |
| 107 | + return (int) (mixLong(state) >>> (48 - 31)); |
| 108 | + } |
| 109 | + |
| 110 | + /** |
| 111 | + * erand48: updates the provided seed array in-place and returns a double in [0,1). |
| 112 | + * The input array must have at least 3 elements. Each element represents an |
| 113 | + * unsigned 16-bit word; Java's short is signed, so mask with 0xFFFF when reading. |
| 114 | + * |
| 115 | + * Array layout: |
| 116 | + * xsubi[0] = low 16 bits |
| 117 | + * xsubi[1] = middle 16 bits |
| 118 | + * xsubi[2] = high 16 bits (most significant) |
| 119 | + * |
| 120 | + * This mirrors the C POSIX definition. |
| 121 | + * |
| 122 | + * @param xsubi short[3] (modified in-place) |
| 123 | + * @return double in [0,1) |
| 124 | + * @throws IllegalArgumentException if xsubi.length < 3 |
| 125 | + */ |
| 126 | + public static double erand48(final short[] xsubi) { |
| 127 | + if (xsubi == null || xsubi.length < 3) { |
| 128 | + throw new IllegalArgumentException("xsubi must be a short[3] (or longer)"); |
| 129 | + } |
| 130 | + // assemble 48-bit state from 3 unsigned 16-bit words (xsubi[0] = least significant) |
| 131 | + long s0 = xsubi[0] & 0xFFFFL; |
| 132 | + long s1 = xsubi[1] & 0xFFFFL; |
| 133 | + long s2 = xsubi[2] & 0xFFFFL; |
| 134 | + long state = (s2 << 32) | (s1 << 16) | s0; |
| 135 | + // update |
| 136 | + state = (state * MULTIPLIER + ADDEND) & MASK48; |
| 137 | + // write back into xsubi (low..high) |
| 138 | + xsubi[0] = (short) (state & 0xFFFFL); |
| 139 | + xsubi[1] = (short) ((state >>> 16) & 0xFFFFL); |
| 140 | + xsubi[2] = (short) ((state >>> 32) & 0xFFFFL); |
| 141 | + return mixLong(state) / (double) (1L << 48); |
| 142 | + } |
| 143 | + |
| 144 | + /** |
| 145 | + * Convenience overload that accepts int[] - each int must be in 0..65535. |
| 146 | + * Useful because Java's short is signed. |
| 147 | + * @param xsubi int[3] (modified in-place) |
| 148 | + * @return double in [0,1) |
| 149 | + * @throws IllegalArgumentException if xsubi.length < 3 |
| 150 | + */ |
| 151 | + public static double erand48FromInts(final int[] xsubi) { |
| 152 | + if (xsubi == null || xsubi.length < 3) { |
| 153 | + throw new IllegalArgumentException("xsubi must be an int[3] (or longer)"); |
| 154 | + } |
| 155 | + long s0 = xsubi[0] & 0xFFFFL; |
| 156 | + long s1 = xsubi[1] & 0xFFFFL; |
| 157 | + long s2 = xsubi[2] & 0xFFFFL; |
| 158 | + long state = (s2 << 32) | (s1 << 16) | s0; |
| 159 | + state = (state * MULTIPLIER + ADDEND) & MASK48; |
| 160 | + xsubi[0] = (int) (state & 0xFFFFL); |
| 161 | + xsubi[1] = (int) ((state >>> 16) & 0xFFFFL); |
| 162 | + xsubi[2] = (int) ((state >>> 32) & 0xFFFFL); |
| 163 | + return mixLong(state) / (double) (1L << 48); |
| 164 | + } |
| 165 | +} |
| 166 | + |
0 commit comments