add p345 solution

* p334:
  add working Python solution for p334
  add invariant version of bean game (bazillion times faster)
  add bean game simulator
  add script to simulate game
  add script to calculate sequence of inputs

java/Problem029.java

@@ -1,4 +1,4 @@
-iimport java.math.BigInteger;
+import java.math.BigInteger;
 import java.util.HashSet;
 import java.util.Set;
 import java.util.concurrent.TimeUnit; // Optional: for more readable time output
@@ -8,7 +8,7 @@ import java.util.concurrent.TimeUnit; // Optional: for more readable time output
  * for 2 <= a <= 100 and 2 <= b <= 100.
  * This implementation runs on a single core/thread.
  */
-public class DistinctPowers {
+public class Problem029 {
 
     // Define the inclusive range for base 'a'
     private static final int MIN_A = 2;

java/Problem323.java

@@ -0,0 +1,94 @@
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.math.RoundingMode;
+
+/**
+ * Calculates the expected number of steps N until a 32-bit integer,
+ * initially zero, becomes all ones (2^32 - 1) by repeatedly performing
+ * a bitwise OR with a sequence of random 32-bit integers.
+ *
+ * The expected value is calculated using the formula:
+ * E[N] = sum_{j=0 to infinity} (1 - (1 - (1/2)^j)^M)
+ * where M = 32.
+ *
+ * This program approximates the sum by summing a fixed number of terms
+ * using high-precision BigDecimal arithmetic.
+ */
+public class Problem323 {
+
+    public static void main(String[] args) {
+
+        // --- Configuration ---
+        final int M = 32; // Exponent (number of bits)
+        // Number of terms in the sum (j=0 to ITERATIONS-1).
+        // Chosen based on analysis to provide sufficient precision.
+        final int ITERATIONS = 70;
+        // Precision for BigDecimal calculations (number of significant digits)
+        // Chosen to be >> 10
+        final int PRECISION = 50;
+        // Final scale for the result (decimal places)
+        final int FINAL_SCALE = 10;
+        // Rounding mode for calculations and final result
+        final RoundingMode ROUNDING_MODE = RoundingMode.HALF_UP;
+
+        // --- Setup High Precision ---
+        // MathContext defines precision and rounding rules for operations
+        MathContext mc = new MathContext(PRECISION, ROUNDING_MODE);
+
+        // BigDecimal constants for efficiency and clarity
+        final BigDecimal ZERO = BigDecimal.ZERO;
+        final BigDecimal ONE = BigDecimal.ONE;
+        final BigDecimal TWO = new BigDecimal("2");
+
+        // --- Calculation ---
+        BigDecimal expectedValueSum = ZERO;
+        // Initialize halfPowerJ to (1/2)^0 = 1
+        BigDecimal halfPowerJ = ONE;
+
+        System.out.println("Calculating E[N] using " + ITERATIONS + " terms with precision " + PRECISION + "...");
+
+        for (int j = 0; j < ITERATIONS; j++) {
+            // Calculate the term Tj = 1 - (1 - (1/2)^j)^M
+
+            // 1. Calculate base: (1 - (1/2)^j)
+            BigDecimal termBase = ONE.subtract(halfPowerJ);
+
+            // 2. Calculate base^M: (1 - (1/2)^j)^M
+            // Use pow(int n, MathContext mc) for controlled precision during exponentiation
+            BigDecimal termBasePowM = termBase.pow(M, mc);
+
+            // 3. Calculate Tj = 1 - base^M
+            BigDecimal termTj = ONE.subtract(termBasePowM);
+
+            // 4. Add Tj to the running sum
+            expectedValueSum = expectedValueSum.add(termTj, mc);
+
+            // 5. Prepare (1/2)^j for the *next* iteration (j+1)
+            // (1/2)^(j+1) = (1/2)^j / 2
+            halfPowerJ = halfPowerJ.divide(TWO, mc);
+
+            // Optional: Print progress every 10 iterations
+            // if ((j + 1) % 10 == 0) {
+            //     System.out.println("j=" + j + ", Current Sum=" + expectedValueSum.round(new MathContext(FINAL_SCALE + 2)));
+            // }
+        }
+
+        // System.out.println("Calculation complete.");
+
+        // --- Round final result ---
+        BigDecimal finalResult = expectedValueSum.setScale(FINAL_SCALE, ROUNDING_MODE);
+
+        // --- Output ---
+        System.out.println("\nCalculated Expected Value E[N]:");
+        // Use toPlainString() to avoid potential scientific notation
+        System.out.println(finalResult.toPlainString());
+
+        // --- Sanity Check (Approximate value based on theoretical analysis) ---
+        double log2_M = Math.log(M) / Math.log(2.0);
+        double gamma = 0.57721566490153286060; // Euler-Mascheroni constant
+        double ln2 = Math.log(2.0);
+        double approxExpectedValue = log2_M + gamma / ln2 + 0.5;
+        System.out.println("\nApproximate theoretical value (for comparison): " + String.format("%.10f", approxExpectedValue));
+
+    }
+}

java/Problem345.java

@@ -0,0 +1,546 @@
+import java.util.Arrays;
+import java.util.List;
+import java.util.ArrayList;
+
+public class Problem345 {
+
+    // Method to calculate the maximum sum using the Hungarian Algorithm
+    public static double calculateMaxMatrixSum(double[][] matrix) {
+        int N = matrix.length;
+        if (N == 0 || matrix[0].length != N) {
+            throw new IllegalArgumentException("Matrix must be square and non-empty.");
+        }
+
+        // 1. Create the cost matrix (negated values)
+        double[][] costMatrix = new double[N][N];
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                costMatrix[i][j] = -matrix[i][j];
+            }
+        }
+
+        // 2. Solve the assignment problem (minimization on cost matrix)
+        // This part requires an actual Hungarian Algorithm implementation.
+        // Let's assume we have a class 'HungarianSolver'
+        HungarianSolver solver = new HungarianSolver(costMatrix);
+        int[] assignment = solver.execute(); // Returns assignment: assignment[row] = col
+
+        // 3. Calculate the minimum cost from the assignment
+        double minCost = 0;
+        for (int i = 0; i < N; i++) {
+            if (assignment[i] != -1) { // Check if row i has an assignment
+                 minCost += costMatrix[i][assignment[i]];
+            } else {
+                 // Handle case where no assignment is found (should not happen for square matrix if algo is correct)
+                 System.err.println("Warning: No assignment found for row " + i);
+            }
+           
+        }
+         // Or, if the solver directly provides the min cost:
+         // double minCost = solver.getOptimalCost(); 
+
+        // 4. The maximum sum is the negative of the minimum cost
+        return -minCost;
+    }
+
+    public static void main(String[] args) {
+        // The 15x15 matrix from the image
+        double[][] matrix = {
+            {  7,  53, 183, 439, 863, 497, 383, 563,  79, 973, 287,  63, 343, 169, 583},
+            {627, 343, 773, 959, 943, 767, 473, 103, 699, 303, 957, 703, 583, 639, 913},
+            {447, 283, 463,  29,  23, 487, 463, 993, 119, 883, 327, 493, 423, 159, 743},
+            {217, 623,   3, 399, 853, 407, 103, 983,  89, 463, 290, 516, 212, 462, 350},
+            {960, 376, 682, 962, 300, 780, 486, 502, 912, 800, 250, 346, 172, 812, 350},
+            {870, 456, 192, 162, 593, 473, 915,  45, 989, 873, 823, 965, 425, 329, 803},
+            {973, 965, 905, 919, 133, 673, 665, 235, 509, 613, 673, 815, 165, 992, 326},
+            {322, 148, 972, 962, 286, 255, 941, 541, 265, 323, 925, 281, 601,  95, 973},
+            {445, 721,  11, 525, 473,  65, 511, 164, 138, 672,  18, 428, 154, 448, 848},
+            {414, 456, 310, 312, 798, 104, 566, 520, 302, 248, 694, 976, 430, 392, 198},
+            {184, 829, 373, 181, 631, 101, 969, 613, 840, 740, 778, 458, 284, 760, 390},
+            {821, 461, 843, 513,  17, 901, 711, 993, 293, 157, 274,  94, 192, 156, 574},
+            { 34, 124,   4, 878, 450, 476, 712, 914, 838, 669, 875, 299, 823, 329, 699},
+            {815, 559, 813, 459, 522, 788, 168, 586, 966, 232, 308, 833, 251, 631, 107},
+            {813, 883, 451, 509, 615,  77, 281, 613, 459, 205, 380, 274, 302,  35, 805}
+        };
+
+        try {
+            double maxMatrixSum = calculateMaxMatrixSum(matrix);
+            // Use long for the result if expecting large integer sums
+            System.out.println("The Matrix Sum is: " + (long)maxMatrixSum); 
+        } catch (Exception e) {
+            System.err.println("Error calculating Matrix Sum: " + e.getMessage());
+            e.printStackTrace();
+        }
+    }
+}
+
+
+
+class HungarianSolver {
+    private double[][] costMatrix;
+    private int rows, cols;
+    private int[] assignment; // assignment[row] = assigned column
+    private double minCost;
+    private boolean[] rowCovered;
+    private boolean[] colCovered;
+    private double[][] originalCostMatrix;
+    private int[][] marked; // 0 = unmarked, 1 = starred, 2 = primed
+
+    // For numerical stability when comparing doubles
+    private static final double EPSILON = 1e-10;
+    private int originalRows, originalCols;
+
+    public HungarianSolver(double[][] matrix) {
+        // Store original dimensions
+        this.originalRows = matrix.length;
+        this.originalCols = matrix[0].length;
+
+        // Clone original matrix
+        this.originalCostMatrix = new double[originalRows][originalCols];
+        for (int i = 0; i < originalRows; i++) {
+            System.arraycopy(matrix[i], 0, originalCostMatrix[i], 0, originalCols);
+        }
+
+        // Make matrix square - use max dimension
+        int maxDim = Math.max(originalRows, originalCols);
+        this.rows = maxDim;
+        this.cols = maxDim;
+
+        this.costMatrix = new double[maxDim][maxDim];
+
+        // Fill with large values (effectively infinity for the algorithm)
+        for (int i = 0; i < maxDim; i++) {
+            for (int j = 0; j < maxDim; j++) {
+                costMatrix[i][j] = Double.MAX_VALUE / 2; // Avoid overflow in calculations
+            }
+        }
+
+        // Copy actual values
+        for (int i = 0; i < originalRows; i++) {
+            for (int j = 0; j < originalCols; j++) {
+                costMatrix[i][j] = matrix[i][j];
+            }
+        }
+
+        this.assignment = new int[maxDim];
+        this.rowCovered = new boolean[maxDim];
+        this.colCovered = new boolean[maxDim];
+        this.marked = new int[maxDim][maxDim];
+    }
+
+    public int[] execute() {
+        // Initialize assignments
+        Arrays.fill(assignment, -1);
+
+        // Clear markings
+        for (int i = 0; i < rows; i++) {
+            Arrays.fill(marked[i], 0);
+        }
+
+        // Step 1: Reduce the matrix
+        reduceMatrix();
+
+        // Step 2: Initial star zeros
+        initialStarZeros();
+
+        // Main loop
+        boolean done = false;
+        while (!done) {
+            // Step 3: Cover starred columns
+            coverColumns();
+
+            // Check if we're done
+            int coveredColumns = countCoveredColumns();
+            if (coveredColumns >= Math.min(rows, cols)) {
+                done = true;
+                continue;
+            }
+
+            // Step 4 & 5: Prime zeros and augment path
+            boolean pathFound = findAugmentingPath();
+
+            if (pathFound) {
+                // Step 5: Augment path was successful, go back to step 3
+                clearCovers();
+                clearPrimes();
+            } else {
+                // Step 6: Adjust matrix
+                adjustMatrix();
+            }
+        }
+
+        // Extract the assignment from the starred zeros
+        int[] result = new int[originalRows];
+        Arrays.fill(result, -1);
+
+        for (int i = 0; i < originalRows; i++) {
+            for (int j = 0; j < originalCols; j++) {
+                if (marked[i][j] == 1) {
+                    result[i] = j;
+                    break;
+                }
+            }
+        }
+
+        // Calculate optimal cost
+        calculateMinCost(result);
+
+        return result;
+    }
+
+    private void reduceMatrix() {
+        // Subtract row minimums
+        for (int i = 0; i < rows; i++) {
+            double minValue = Double.MAX_VALUE;
+            for (int j = 0; j < cols; j++) {
+                if (costMatrix[i][j] < minValue) {
+                    minValue = costMatrix[i][j];
+                }
+            }
+
+            // Only subtract if minimum is finite
+            if (minValue < Double.MAX_VALUE / 2) {
+                for (int j = 0; j < cols; j++) {
+                    if (costMatrix[i][j] < Double.MAX_VALUE / 2) {
+                        costMatrix[i][j] -= minValue;
+                    }
+                }
+            }
+        }
+
+        // Subtract column minimums
+        for (int j = 0; j < cols; j++) {
+            double minValue = Double.MAX_VALUE;
+            for (int i = 0; i < rows; i++) {
+                if (costMatrix[i][j] < minValue) {
+                    minValue = costMatrix[i][j];
+                }
+            }
+
+            // Only subtract if minimum is finite
+            if (minValue < Double.MAX_VALUE / 2) {
+                for (int i = 0; i < rows; i++) {
+                    if (costMatrix[i][j] < Double.MAX_VALUE / 2) {
+                        costMatrix[i][j] -= minValue;
+                    }
+                }
+            }
+        }
+    }
+
+    private void initialStarZeros() {
+        // Track if a row or column already has a star
+        boolean[] rowStarred = new boolean[rows];
+        boolean[] colStarred = new boolean[cols];
+
+        // First pass: try to star as many zeros as possible
+        for (int i = 0; i < rows; i++) {
+            for (int j = 0; j < cols; j++) {
+                if (isZero(costMatrix[i][j]) && !rowStarred[i] && !colStarred[j]) {
+                    marked[i][j] = 1; // Star this zero
+                    rowStarred[i] = true;
+                    colStarred[j] = true;
+                }
+            }
+        }
+    }
+
+    private void coverColumns() {
+        // Reset all covers
+        Arrays.fill(rowCovered, false);
+        Arrays.fill(colCovered, false);
+
+        // Cover all columns with starred zeros
+        for (int j = 0; j < cols; j++) {
+            for (int i = 0; i < rows; i++) {
+                if (marked[i][j] == 1) {
+                    colCovered[j] = true;
+                    break;
+                }
+            }
+        }
+    }
+
+    private int countCoveredColumns() {
+        int count = 0;
+        for (boolean covered : colCovered) {
+            if (covered) count++;
+        }
+        return count;
+    }
+
+    private boolean findAugmentingPath() {
+        int row = -1, col = -1;
+
+        // Find an uncovered zero
+        int[] zeroPos = findUncoveredZero();
+        row = zeroPos[0];
+        col = zeroPos[1];
+
+        if (row == -1) {
+            // No uncovered zeros left
+            return false;
+        }
+
+        // Mark the uncovered zero as primed
+        marked[row][col] = 2;
+
+        // Find a starred zero in the row
+        int starCol = findStarInRow(row);
+
+        if (starCol == -1) {
+            // No starred zero in this row
+            // Start the augmenting path with this zero
+            int[] pathStart = {row, col};
+            List<int[]> path = new ArrayList<>();
+            path.add(pathStart);
+
+            // Augment the matching along this path
+            augmentPath(path);
+            return true;
+        } else {
+            // There is a starred zero in this row
+            // Cover this row and uncover the column with the starred zero
+            rowCovered[row] = true;
+            colCovered[starCol] = false;
+
+            // Continue looking for augmenting path
+            return continueFindingPath();
+        }
+    }
+
+    private boolean continueFindingPath() {
+        List<int[]> path = new ArrayList<>();
+        int pathRow = -1, pathCol = -1;
+        boolean done = false;
+
+        // Start with the last primed zero that doesn't have a starred zero in its row
+        for (int j = 0; j < cols; j++) {
+            if (!colCovered[j]) {
+                for (int i = 0; i < rows; i++) {
+                    if (marked[i][j] == 2) {
+                        // Found a primed zero in uncovered column
+                        int starCol = findStarInRow(i);
+                        if (starCol == -1) {
+                            // No starred zero in this row - this is our starting point
+                            pathRow = i;
+                            pathCol = j;
+                            done = true;
+                            break;
+                        }
+                    }
+                }
+            }
+            if (done) break;
+        }
+
+        if (pathRow == -1) {
+            // No primed zeros with no starred zero in their row
+            // Try to find any uncovered zero
+            int[] zeroPos = findUncoveredZero();
+            int row = zeroPos[0];
+            int col = zeroPos[1];
+
+            if (row == -1) {
+                // No uncovered zeros left
+                return false;
+            }
+
+            // Mark this zero as primed
+            marked[row][col] = 2;
+
+            // Find starred zero in this row
+            int starCol = findStarInRow(row);
+
+            if (starCol == -1) {
+                // No starred zero - start augmenting path here
+                pathRow = row;
+                pathCol = col;
+            } else {
+                // Cover this row and uncover the column with the starred zero
+                rowCovered[row] = true;
+                colCovered[starCol] = false;
+                return continueFindingPath();
+            }
+        }
+
+        if (pathRow != -1) {
+            // Found a starting point for augmenting path
+            path.add(new int[]{pathRow, pathCol});
+            augmentPath(path);
+            return true;
+        }
+
+        return false;
+    }
+
+    private void augmentPath(List<int[]> initialPath) {
+        // The path starts with a primed zero
+        List<int[]> completePath = buildCompletePath(initialPath.get(0));
+
+        // Augment the path
+        for (int[] pos : completePath) {
+            int r = pos[0];
+            int c = pos[1];
+
+            if (marked[r][c] == 1) {
+                // Unstar starred zeros
+                marked[r][c] = 0;
+            } else if (marked[r][c] == 2) {
+                // Star primed zeros
+                marked[r][c] = 1;
+            }
+        }
+
+        // Clear all covers and primes
+        clearCovers();
+        clearPrimes();
+    }
+
+    private List<int[]> buildCompletePath(int[] start) {
+        List<int[]> path = new ArrayList<>();
+        path.add(start);
+
+        int row = start[0];
+        int col = start[1];
+
+        while (true) {
+            // Find a starred zero in the column
+            int starRow = findStarInColumn(col);
+
+            if (starRow == -1) {
+                // No starred zero in this column - path is complete
+                break;
+            }
+
+            // Add starred zero to path
+            path.add(new int[]{starRow, col});
+
+            // Find a primed zero in the row
+            col = findPrimeInRow(starRow);
+
+            // Add primed zero to path
+            path.add(new int[]{starRow, col});
+        }
+
+        return path;
+    }
+
+    private int findStarInRow(int row) {
+        for (int j = 0; j < cols; j++) {
+            if (marked[row][j] == 1) {
+                return j;
+            }
+        }
+        return -1;
+    }
+
+    private int findStarInColumn(int col) {
+        for (int i = 0; i < rows; i++) {
+            if (marked[i][col] == 1) {
+                return i;
+            }
+        }
+        return -1;
+    }
+
+    private int findPrimeInRow(int row) {
+        for (int j = 0; j < cols; j++) {
+            if (marked[row][j] == 2) {
+                return j;
+            }
+        }
+        return -1; // Should not happen in a valid path
+    }
+
+    private int[] findUncoveredZero() {
+        for (int i = 0; i < rows; i++) {
+            if (!rowCovered[i]) {
+                for (int j = 0; j < cols; j++) {
+                    if (!colCovered[j] && isZero(costMatrix[i][j])) {
+                        return new int[]{i, j};
+                    }
+                }
+            }
+        }
+        return new int[]{-1, -1}; // No uncovered zero found
+    }
+
+    private void clearCovers() {
+        Arrays.fill(rowCovered, false);
+        Arrays.fill(colCovered, false);
+    }
+
+    private void clearPrimes() {
+        for (int i = 0; i < rows; i++) {
+            for (int j = 0; j < cols; j++) {
+                if (marked[i][j] == 2) {
+                    marked[i][j] = 0;
+                }
+            }
+        }
+    }
+
+    private void adjustMatrix() {
+        // Find minimum uncovered value
+        double minValue = findMinUncovered();
+
+        // Add minValue to every covered row
+        for (int i = 0; i < rows; i++) {
+            if (rowCovered[i]) {
+                for (int j = 0; j < cols; j++) {
+                    if (costMatrix[i][j] < Double.MAX_VALUE / 2) {
+                        costMatrix[i][j] += minValue;
+                    }
+                }
+            }
+        }
+
+        // Subtract minValue from every uncovered column
+        for (int j = 0; j < cols; j++) {
+            if (!colCovered[j]) {
+                for (int i = 0; i < rows; i++) {
+                    if (costMatrix[i][j] < Double.MAX_VALUE / 2) {
+                        costMatrix[i][j] -= minValue;
+                    }
+                }
+            }
+        }
+    }
+
+    private double findMinUncovered() {
+        double min = Double.MAX_VALUE;
+
+        for (int i = 0; i < rows; i++) {
+            if (!rowCovered[i]) {
+                for (int j = 0; j < cols; j++) {
+                    if (!colCovered[j] && costMatrix[i][j] < min) {
+                        min = costMatrix[i][j];
+                    }
+                }
+            }
+        }
+
+        return min;
+    }
+
+    private boolean isZero(double value) {
+        return Math.abs(value) < EPSILON;
+    }
+
+    public double getOptimalCost() {
+        return minCost;
+    }
+
+    private void calculateMinCost(int[] result) {
+        minCost = 0;
+
+        for (int i = 0; i < originalRows; i++) {
+            int j = result[i];
+            if (j != -1 && j < originalCols) {
+                minCost += originalCostMatrix[i][j];
+            }
+        }
+    }
+}
+

scratch/Misc_301-400/323/ExpectedValueN.java

@@ -0,0 +1,94 @@
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.math.RoundingMode;
+
+/**
+ * Calculates the expected number of steps N until a 32-bit integer,
+ * initially zero, becomes all ones (2^32 - 1) by repeatedly performing
+ * a bitwise OR with a sequence of random 32-bit integers.
+ *
+ * The expected value is calculated using the formula:
+ * E[N] = sum_{j=0 to infinity} (1 - (1 - (1/2)^j)^M)
+ * where M = 32.
+ *
+ * This program approximates the sum by summing a fixed number of terms
+ * using high-precision BigDecimal arithmetic.
+ */
+public class ExpectedValueN {
+
+    public static void main(String[] args) {
+
+        // --- Configuration ---
+        final int M = 32; // Exponent (number of bits)
+        // Number of terms in the sum (j=0 to ITERATIONS-1).
+        // Chosen based on analysis to provide sufficient precision.
+        final int ITERATIONS = 70;
+        // Precision for BigDecimal calculations (number of significant digits)
+        // Chosen to be >> 10
+        final int PRECISION = 50;
+        // Final scale for the result (decimal places)
+        final int FINAL_SCALE = 10;
+        // Rounding mode for calculations and final result
+        final RoundingMode ROUNDING_MODE = RoundingMode.HALF_UP;
+
+        // --- Setup High Precision ---
+        // MathContext defines precision and rounding rules for operations
+        MathContext mc = new MathContext(PRECISION, ROUNDING_MODE);
+
+        // BigDecimal constants for efficiency and clarity
+        final BigDecimal ZERO = BigDecimal.ZERO;
+        final BigDecimal ONE = BigDecimal.ONE;
+        final BigDecimal TWO = new BigDecimal("2");
+
+        // --- Calculation ---
+        BigDecimal expectedValueSum = ZERO;
+        // Initialize halfPowerJ to (1/2)^0 = 1
+        BigDecimal halfPowerJ = ONE;
+
+        System.out.println("Calculating E[N] using " + ITERATIONS + " terms with precision " + PRECISION + "...");
+
+        for (int j = 0; j < ITERATIONS; j++) {
+            // Calculate the term Tj = 1 - (1 - (1/2)^j)^M
+
+            // 1. Calculate base: (1 - (1/2)^j)
+            BigDecimal termBase = ONE.subtract(halfPowerJ);
+
+            // 2. Calculate base^M: (1 - (1/2)^j)^M
+            // Use pow(int n, MathContext mc) for controlled precision during exponentiation
+            BigDecimal termBasePowM = termBase.pow(M, mc);
+
+            // 3. Calculate Tj = 1 - base^M
+            BigDecimal termTj = ONE.subtract(termBasePowM);
+
+            // 4. Add Tj to the running sum
+            expectedValueSum = expectedValueSum.add(termTj, mc);
+
+            // 5. Prepare (1/2)^j for the *next* iteration (j+1)
+            // (1/2)^(j+1) = (1/2)^j / 2
+            halfPowerJ = halfPowerJ.divide(TWO, mc);
+
+            // Optional: Print progress every 10 iterations
+            // if ((j + 1) % 10 == 0) {
+            //     System.out.println("j=" + j + ", Current Sum=" + expectedValueSum.round(new MathContext(FINAL_SCALE + 2)));
+            // }
+        }
+
+        // System.out.println("Calculation complete.");
+
+        // --- Round final result ---
+        BigDecimal finalResult = expectedValueSum.setScale(FINAL_SCALE, ROUNDING_MODE);
+
+        // --- Output ---
+        System.out.println("\nCalculated Expected Value E[N]:");
+        // Use toPlainString() to avoid potential scientific notation
+        System.out.println(finalResult.toPlainString());
+
+        // --- Sanity Check (Approximate value based on theoretical analysis) ---
+        double log2_M = Math.log(M) / Math.log(2.0);
+        double gamma = 0.57721566490153286060; // Euler-Mascheroni constant
+        double ln2 = Math.log(2.0);
+        double approxExpectedValue = log2_M + gamma / ln2 + 0.5;
+        System.out.println("\nApproximate theoretical value (for comparison): " + String.format("%.10f", approxExpectedValue));
+
+    }
+}

scratch/Misc_301-400/334/game.py

@@ -0,0 +1,164 @@
+import collections
+import time
+
+def simulate_bean_game(initial_bowls):
+    """
+    Simulates the bean game efficiently using a dictionary for state
+    and a set for eligible move tracking.
+
+    Args:
+        initial_bowls: A dictionary {index: count} representing the initial
+                       number of beans in each bowl. Bowls not in the dict
+                       are assumed to be empty.
+
+    Returns:
+        A dictionary {index: count} representing the final state where
+        all counts are 1. Keys with 0 count are excluded.
+        Returns it sorted by index for consistency.
+    """
+    # Use defaultdict for easier handling of empty bowls
+    state = collections.defaultdict(int)
+    state.update(initial_bowls)
+
+    # Keep track of bowls that can make a move (count >= 2)
+    # A set provides efficient add, remove, and arbitrary element retrieval (pop)
+    eligible_indices = set()
+    for i, count in initial_bowls.items():
+        if count >= 2:
+            eligible_indices.add(i)
+
+    # --- Simulation loop ---
+    while eligible_indices:
+        # Pick an index to process. Using pop() from set is arbitrary but fine.
+        # The uniqueness of the final state ensures the order doesn't matter.
+        i = eligible_indices.pop()
+
+        # Since we manage the set carefully, we can assume state[i] >= 2 here.
+        # If state[i] somehow dropped below 2 without i being removed from
+        # eligible_indices, that would indicate a bug in the update logic.
+
+        # --- Perform the move ---
+        # Store previous counts of neighbors *before* incrementing
+        prev_neighbor_counts = {
+            i - 1: state[i - 1], # Using defaultdict simplifies access
+            i + 1: state[i + 1]
+        }
+
+        state[i] -= 2
+        state[i - 1] += 1
+        state[i + 1] += 1
+        # --- Move performed ---
+
+
+        # --- Update eligibility ---
+
+        # Check bowl i (the one where the move happened)
+        if state[i] >= 2:
+            eligible_indices.add(i) # Still eligible, add it back
+
+        # Check bowl i-1 (left neighbor)
+        idx_minus_1 = i - 1
+        # If the count was < 2 before and is now >= 2, it becomes eligible
+        if prev_neighbor_counts[idx_minus_1] < 2 and state[idx_minus_1] >= 2:
+            eligible_indices.add(idx_minus_1)
+
+        # Check bowl i+1 (right neighbor)
+        idx_plus_1 = i + 1
+        # If the count was < 2 before and is now >= 2, it becomes eligible
+        if prev_neighbor_counts[idx_plus_1] < 2 and state[idx_plus_1] >= 2:
+            eligible_indices.add(idx_plus_1)
+
+        # Note: We don't explicitly remove bowls when their count drops below 2 here.
+        # They are implicitly removed by not being added back to eligible_indices
+        # after their turn (via pop) if their count drops below 2.
+
+    # --- Simulation complete ---
+
+    # Clean up the final state: remove zero-count bowls explicitly
+    # (although defaultdict wouldn't store them unless explicitly set to 0)
+    # Also ensures the result is a standard dict.
+    final_state = {idx: count for idx, count in state.items() if count > 0}
+
+    # Verify final state properties (all counts should be 1)
+    for idx, count in final_state.items():
+        if count != 1:
+             # This indicates a problem if the game theory holds true.
+             raise ValueError(f"Simulation ended unexpectedly: Bowl {idx} has count {count}")
+
+    # Return sorted dictionary for consistent output
+    return dict(sorted(final_state.items()))
+
+# --- Helper Function for Invariants (for testing) ---
+def calculate_invariants(bowls):
+    """Calculates the total number of beans (L) and the moment (M)."""
+    L = sum(bowls.values())
+    M = sum(i * count for i, count in bowls.items())
+    return L, M
+
+# --- Run the simulation for the given initial state ---
+def simulate(initial_state, debug=True):
+    if debug:
+        print("Initial State:", initial_state)
+    start_time = time.time()
+    final_state = simulate_bean_game(initial_state)
+    end_time = time.time()
+    if debug:
+        print("Final State:", sorted([int(k) for k in final_state.keys() if final_state[k]==1]))
+    print(f"(Took {end_time - start_time:.6f} seconds)")
+    
+    # Verify Invariants
+    initial_L, initial_M = calculate_invariants(initial_state)
+    final_L, final_M = calculate_invariants(final_state)
+    if debug:
+        print(f"Initial Invariants: L={initial_L}, M={initial_M}")
+        print(f"Final Invariants:   L={final_L}, M={final_M}")
+    assert initial_L == final_L and initial_M == final_M
+    if debug:
+        print("Invariants preserved: Yes")
+    print()
+
+
+print("--- Example from Problem ---")
+# Bowls at index 2 and 3 with 2 and 3 beans respectively
+initial_state_example = {2: 2, 3: 3}
+simulate(initial_state_example)
+
+print("\n--- Single Bowl Example ---")
+initial_state_complex = {7: 100}
+simulate(initial_state_complex)
+
+print("\n--- Complex Example ---")
+initial_state_complex = {0: 10, 80: 10}
+simulate(initial_state_complex)
+
+print("\n--- Spacing Example ---")
+initial_state_complex = {0: 10, 1: 10}
+simulate(initial_state_complex)
+
+initial_state_complex = {0: 10, 2: 10}
+simulate(initial_state_complex)
+
+initial_state_complex = {0: 10, 3: 10}
+simulate(initial_state_complex)
+
+initial_state_complex = {0: 10, 4: 10}
+simulate(initial_state_complex)
+
+initial_state_complex = {0: 10, 5: 10}
+simulate(initial_state_complex)
+
+### print("\n--- Large Number Examples ---")
+### # Put N beans in bowl 0. Expect final state to be spread out.
+### 
+### print("N = 10:")
+### initial_state_large = {0: 10}
+### simulate(initial_state_large, False)
+### 
+### print("N = 100:")
+### initial_state_large = {0: 100}
+### simulate(initial_state_large, False)
+### 
+### print("N = 1000:")
+### initial_state_large = {0: 1000}
+### simulate(initial_state_large, False)
+

scratch/Misc_301-400/334/invariants.py

@@ -0,0 +1,99 @@
+import time
+import math
+
+def calculate_invariants(bowls):
+    """Calculates the total number of beans (L) and the moment (M)."""
+    L = sum(bowls.values())
+    M = sum(i * count for i, count in bowls.items())
+    return L, M
+
+def calculate_final_state_directly(initial_bowls):
+    """
+    Calculates the final bean game state directly using invariants L and M,
+    based on minimizing the sum of squared indices.
+
+    Args:
+        initial_bowls: A dictionary {index: count} representing the initial state.
+
+    Returns:
+        A dictionary {index: 1} representing the final state, sorted by index.
+    """
+    # 1. Calculate Invariants
+    L, M = calculate_invariants(initial_bowls)
+
+    if L == 0:
+        return {}
+
+    # 2. Determine Base Moment (for interval 0 to L-1)
+    # Sum of 0 + 1 + ... + (L-1)
+    M_base = L * (L - 1) // 2
+
+    # 3. Calculate Moment Distribution
+    M_offset = M - M_base
+    q = M_offset // L  # Base shift for all indices
+    r = M_offset % L   # Number of indices getting an extra +1 shift
+
+    # Handle potential negative remainder in Python's % for negative M_offset
+    # Ensure 0 <= r < L
+    if r < 0:
+        r += L
+        q -= 1 # Adjust quotient accordingly if remainder was negative
+
+    # 4. Construct Final Indices
+    final_indices = set()
+    # Indices shifted by q
+    for i in range(L - r):
+        final_indices.add(i + q)
+    # Indices shifted by q+1 (the top r elements of the base sequence)
+    for i in range(L - r, L):
+        final_indices.add(i + q + 1)
+
+    # 5. Return Result dictionary
+    final_state = {idx: 1 for idx in final_indices}
+
+    # Return sorted dictionary for consistent output
+    return dict(sorted(final_state.items()))
+
+# --- Example Usage ---
+print("--- Example from Problem ---")
+initial_state_example = {2: 2, 3: 3}
+print("Initial State:", initial_state_example)
+L_ex, M_ex = calculate_invariants(initial_state_example)
+print(f"Invariants: L={L_ex}, M={M_ex}")
+start_time = time.time()
+final_state_example = calculate_final_state_directly(initial_state_example)
+end_time = time.time()
+print("Calculated Final State:", final_state_example)
+print(f"(Took {end_time - start_time:.6f} seconds)")
+# Verify matches previous simulation result {0: 1, 1: 1, 3: 1, 4: 1, 5: 1} - Yes!
+
+# --- Complex Example ---
+print("\n--- Complex Example ---")
+initial_state_complex = {0: 5, 5: 4} # L=9, M=0*5+5*4=20
+print("Initial State:", initial_state_complex)
+L_comp, M_comp = calculate_invariants(initial_state_complex)
+print(f"Invariants: L={L_comp}, M={M_comp}")
+start_time = time.time()
+final_state_complex = calculate_final_state_directly(initial_state_complex)
+end_time = time.time()
+print("Calculated Final State:", final_state_complex)
+print(f"(Took {end_time - start_time:.6f} seconds)")
+
+# --- Large Number Example ---
+print("\n--- Large Number Example ---")
+initial_state_large = {0: 1000000} # L=1,000,000, M=0
+print("Initial State:", initial_state_large)
+L_large, M_large = calculate_invariants(initial_state_large)
+print(f"Invariants: L={L_large}, M={M_large}")
+start_time = time.time()
+final_state_large = calculate_final_state_directly(initial_state_large)
+end_time = time.time()
+# Don't print the full state, just confirm calculation is fast
+print(f"Calculated Final State (Size: {len(final_state_large)})")
+# Optionally print first/last few elements
+keys = list(final_state_large.keys())
+if len(keys) > 10:
+    print(f"Example keys: {keys[:5]} ... {keys[-5:]}")
+else:
+     print(f"Final State Keys: {keys}")
+print(f"(Took {end_time - start_time:.6f} seconds)")

scratch/Misc_301-400/334/problem334.py

@@ -0,0 +1,143 @@
+# Final Python code (identical to user's provided code, validated as correct implementation of the derived algorithm)
+import functools
+import time
+import math
+
+t0 = 123456
+
+@functools.lru_cache(maxsize=None) # Use maxsize=None for potentially deep recursion
+def t(i):
+    if i < 0: # Added safety check
+        raise ValueError("Index i cannot be negative for t(i)")
+    if i == 0:
+        return t0
+
+    tim1 = t(i - 1)
+    if tim1 % 2 == 0:
+        # Use integer division
+        return tim1 // 2
+    else:
+        # Floor division is implicit with // for positive numbers
+        # For potentially negative t(i-1) (unlikely here), math.floor might differ
+        # Python's // behaves as floor division
+        # return math.floor(tim1 / 2) ^ 926252 # Original used floor + XOR -1
+
+        # Re-checking image: It looks like floor(t(i-1)/2) XOR 926252. Let's implement that.
+        # If t(i-1) is odd, t(i-1)//2 performs floor division.
+        return (tim1 // 2) ^ 926252
+
+@functools.lru_cache(maxsize=None)
+def b(i):
+    if i <= 0: # Changed condition to include 0 as invalid bowl index per problem
+         raise ValueError("Index i must be positive for b(i)")
+    # Use 2**11 directly instead of pow for clarity/efficiency
+    return t(i) % 2048 + 1
+
+def calculate_initial_stats(initial_bowls):
+    """Calculates L (beans), M (moment), and Q (quadratic sum) for the initial state."""
+    L = 0
+    M = 0
+    Q_initial = 0
+    for i, count in initial_bowls.items():
+        # Ensure keys and counts are integers
+        idx = int(i)
+        cnt = int(count)
+        L += cnt
+        M += idx * cnt
+        Q_initial += (idx**2) * cnt
+    return L, M, Q_initial
+
+def calculate_final_state_and_moves(initial_bowls):
+    """
+    Calculates the final bean game state directly using invariants L and M,
+    and also calculates the total number of moves required.
+    """
+    # 1. Calculate Initial Stats
+    L, M, Q_initial = calculate_initial_stats(initial_bowls)
+
+    if L == 0:
+        return {}, 0
+
+    # 2. Determine Final State Indices (Minimizes Q)
+    # Base moment (for interval 0 to L-1)
+    M_base = L * (L - 1) // 2
+    # Moment offset needed
+    M_offset = M - M_base
+    # Base shift (q) and number getting extra shift (r)
+    # Use explicit integer conversion just in case L or M_offset become non-int elsewhere (unlikely here)
+    L_int = int(L)
+    M_offset_int = int(M_offset)
+    q = M_offset_int // L_int
+    r = M_offset_int % L_int
+
+    # Ensure 0 <= r < L for Python's % behavior with negative numbers
+    if r < 0:
+        r += L_int
+        q -= 1
+
+    final_indices = set()
+    # Indices shifted by q
+    # Use L_int, r, q which are guaranteed integers
+    for i in range(L_int - r):
+        final_indices.add(i + q)
+    # Indices shifted by q+1 (the top r elements of the base sequence)
+    for i in range(L_int - r, L_int):
+        final_indices.add(i + q + 1)
+
+    # 3. Calculate Final Q
+    # Use explicit loop for clarity with large numbers, sum() is also fine
+    Q_final = 0
+    for idx in final_indices:
+        Q_final += idx**2 # Or idx*idx
+
+    # 4. Calculate Number of Moves
+    delta_Q = Q_final - Q_initial
+    if delta_Q < 0:
+         raise ValueError(f"Q_final ({Q_final}) should not be less than Q_initial ({Q_initial})")
+    if delta_Q % 2 != 0:
+         raise ValueError(f"Change in Q ({delta_Q}) must be even. Q_final={Q_final}, Q_initial={Q_initial}")
+
+    num_moves = delta_Q // 2
+
+    # 5. Return Result
+    final_state = {idx: 1 for idx in final_indices}
+    # Sorting the final dictionary is expensive for L=1.5M and not required by problem
+    # Return the set of indices and moves for efficiency, or unsorted dict
+    # Let's return unsorted dict as per original structure, but note the cost.
+    # return dict(sorted(final_state.items())), num_moves
+    return final_state, num_moves # Return unsorted dict for performance
+
+print("Generating initial state...")
+initial_state_large = {}
+# Generate b(i) for i from 1 to 1500
+# Increase recursion depth limit if needed, although caching should help
+# import sys
+# sys.setrecursionlimit(2000) # Example
+for i in range(1, 1500 + 1):
+    initial_state_large[i] = b(i)
+
+print("Initial state generated.")
+keys = list(initial_state_large.keys())
+if len(keys) > 10:
+    print(f"Example initial keys: {keys[:5]} ... {keys[-5:]}")
+
+print("Calculating invariants...")
+L_large, M_large, Q_init_large = calculate_initial_stats(initial_state_large)
+print(f"Invariants: L={L_large}, M={M_large}, Q_initial={Q_init_large}")
+
+print("Calculating final state and moves...")
+start_time = time.time()
+final_state_large_unsorted, moves_large = calculate_final_state_and_moves(initial_state_large)
+end_time = time.time()
+
+print(f"Calculation finished (Took {end_time - start_time:.6f} seconds)")
+print(f"Calculated Final State (Size: {len(final_state_large_unsorted)})")
+
+# Get keys for printing examples, sorting here just for display
+keys = sorted(final_state_large_unsorted.keys())
+if len(keys) > 10:
+    print(f"Example final keys: {keys[:5]} ... {keys[-5:]}")
+else:
+     print(f"Final State Keys: {keys}")
+
+print(f"Number of Moves: {moves_large}")

scratch/Misc_301-400/334/sequence.py

@@ -0,0 +1,25 @@
+import math
+import functools
+
+
+t0 = 123456
+
+@functools.lru_cache
+def t(i):
+    if i==0:
+        return t0
+
+    tim1 = t(i-1)
+    if tim1%2==0:
+        return tim1//2
+    else:
+        return math.floor(tim1/2)-1 ^ 926252
+
+@functools.lru_cache
+def b(i):
+    return t(i)%(pow(2, 11)) + 1
+
+s = 0
+for i in range(1500):
+    s += i
+print(math.log(s,2))