pruning_algorithm.js

/**
 * IVI (Integer Vector Inversion) Algorithm for Prime Factorization with Global Feasibility Pruning
 * 
 * This module implements the IVI algorithm with comprehensive global feasibility pruning.
 * 
 * Requires common.js to be loaded first.
 * 
 * Key properties:
 * - Linear complexity: O(n) where n is the number of digits
 * - Bounded state space: At most 50 admissible digit pairs per position (empirical)
 * - Deterministic: Logically decidable (LD) process
 * - Global pruning: Mathematically sound bounds eliminate impossible branches early
 * 
 * @module pruning_algorithm
 */

(function() {
    'use strict';

/**
 * Compares two partial numbers represented as LSD-first digit arrays.
 * Returns: -1 if p < q, 0 if p == q, 1 if p > q
 * 
 * @param {number[]} p_history - p digits (LSD-first)
 * @param {number[]} q_history - q digits (LSD-first)
 * @returns {number} Comparison result
 */
const M11 = 11n;

/** Mod-11 inverse for a in 1..10 (11 is prime: a^{-1} = a^{9} mod 11). */
function mod11Inverse(a) {
    const x = Number(a);
    if (x <= 0 || x >= 11) return 1n;
    let t = 1;
    for (let i = 0; i < 9; i++) t = (t * x) % 11;
    return BigInt(t);
}

/**
 * Build set of possible residues mod 11 for (P_value + powerK * A) where A in [0, M].
 * Step = powerK mod 11; set has at most 11 elements: { (p0 + step*j) mod 11 : j = 0..10 }.
 */
function possibleResiduesMod11(valueMod11, powerKMod11, M) {
    const set = new Set();
    const step = powerKMod11 === 0n ? 0n : powerKMod11;
    const cap = M < 10n ? M : 10n;
    for (let j = 0n; j <= cap; j++) {
        const r = (valueMod11 + step * j) % M11;
        set.add(r);
    }
    return set;
}

function comparePartialNumbers(p_history, q_history) {
    const len = Math.max(p_history.length, q_history.length);
    // Compare from most significant digit (end of array) to least significant
    for (let i = len - 1; i >= 0; i--) {
        const p_digit = i < p_history.length ? p_history[i] : 0;
        const q_digit = i < q_history.length ? q_history[i] : 0;
        if (p_digit < q_digit) return -1;
        if (p_digit > q_digit) return 1;
    }
    return 0; // Equal
}

/**
 * Global feasibility pruning: checks if a branch can possibly lead to a valid factorization.
 * 
 * Implements comprehensive tightening per specification:
 * 1. Exact termination rule - validates final state (P·Q = N; accept P > Q and swap on return)
 * 2. Immediate overshoot check - product already exceeds N
 * 3. Sqrt-based hard bound - P > sqrtN implies P*Q > N
 * 4. Explicit growth envelope (replaces rectangle bounds)
 * 5. Minimum contribution pruning
 * 6. Linear-term gap feasibility
 * 7. Upper tail tightening (division-based coupling)
 * 8. Length split feasibility (reserved)
 * (9. Mod-11 feasibility: tried and failed — does not help; left commented out.)
 *
 * All arithmetic uses BigInt. No floating point, no heuristics.
 * Checks ordered for fail-fast performance.
 * 
 * @param {Object} state - Current branch state
 * @param {bigint} P_value - Current partial p value as BigInt
 * @param {bigint} Q_value - Current partial q value as BigInt
 * @param {bigint} N - Target number to factorize as BigInt
 * @param {number} k - Number of digits processed (after adding current digits)
 * @param {number} totalDigits - Total number of digits in N
 * @param {bigint} sqrtN - Precomputed floor(sqrt(N)) as BigInt
 * @returns {number|boolean} Pruning step number (1-8) if pruned, 0 if feasible, true if exact match at termination
 */
function globalFeasible(state, P_value, Q_value, N, k, totalDigits, sqrtN) {
    const remainingDigits = totalDigits - k;

    // #1. Exact Termination Rule
    // If remaining === 0, accept only if Pk * Qk === N
    // (carry_in === 0 is checked separately in workFunction)
    // Symmetry: we accept P > Q at termination and swap when returning (see stepAlgorithm).
    if (remainingDigits === 0) {
        // Reject trivial factor 1
        if (P_value <= 1n || Q_value <= 1n) {
            return 1; // Pruned by step 1
        }
        
        // Reject actual zero values
        if (P_value === 0n || Q_value === 0n) {
            return 1; // Pruned by step 1
        }
        
        // Exact product match required
        const product = P_value * Q_value;
        return product === N ? true : 1; // true if match, 1 if pruned
    }
    
    // #2. Immediate Overshoot Check (fail fast)
    // Compute base product once and reuse
    const base = P_value * Q_value;
    if (base > N) {
        return 2; // Pruned by step 2
    }
    
    // #3. Sqrt-Based Hard Bound (Strengthened with ordering)
    // Since we return with P ≤ Q (swap), if P > sqrtN then P*Q > N
    if (P_value > sqrtN) {
        return 3; // Pruned by step 3
    }
    
    // Compute gap once
    const gap = N - base;
    
    // If gap < 0, we already pruned above, but check for safety
    if (gap < 0n) {
        return 2; // Pruned by step 2 (overshoot)
    }
    
    // Compute M = 10^remaining - 1 (maximum tail value)
    const M = powerOf10(remainingDigits) - 1n;
    const powerK = powerOf10(k);
    const power2K = powerOf10(2 * k);
    
    // #4. Explicit Growth Envelope (replaces Pmax * Qmax)
    // Maximum future contribution:
    // maxContribution = 10^k (Pk*M + Qk*M) + 10^(2k) M*M
    const maxLinearTerm = powerK * (P_value * M + Q_value * M);
    const maxQuadraticTerm = power2K * (M * M);
    const maxContribution = maxLinearTerm + maxQuadraticTerm;
    
    if (maxContribution < gap) {
        return 4; // Pruned by step 4
    }
    
    // #5. Minimum Contribution Pruning
    // Minimum occurs at A=0, B=0, but tighten for remaining === 1
    // However, if gap is already 0 or very small, we might be at the solution
    // Only apply minimum contribution check if gap is significant
    if (gap > 0n) {
        let minA = 0n;
        let minB = 0n;
        if (remainingDigits === 1) {
            // Highest digit must be ≥ 1, and with P ≤ Q, both need at least 1
            minA = 1n;
            minB = 1n;
        }
        
        const minContribution = powerK * (P_value * minB + Q_value * minA) + power2K * (minA * minB);
        if (minContribution > gap) {
            return 5; // Pruned by step 5
        }
    }
    
    // #6. Linear-Term Gap Feasibility (redundant with 4, for stats)
    if (maxLinearTerm + maxQuadraticTerm < gap) {
        return 4; // Pruned by step 4 (already checked above)
    }
    
    // #7. Upper Tail Tightening (Division-Based Coupling)
    // If Pmax ≤ sqrtN, then Q must satisfy Q ≥ ceil(N / Pmax)
    const Pmax = P_value + M * powerK;
    if (Pmax <= sqrtN) {
        // Compute minRequiredQ = ceil(N / Pmax)
        // ceil(a/b) = (a + b - 1) / b for BigInt
        const minRequiredQ = (N + Pmax - 1n) / Pmax;
        const Qmax = Q_value + M * powerK;
        if (Qmax < minRequiredQ) {
            return 7; // Pruned by step 7
        }
    }
    
    // #8. Length Split Feasibility (reserved)
    // Valid factors must satisfy: lenP + lenQ = lenN OR lenN + 1
    // More sophisticated length tracking could be added.

    // #9. Mod-11 feasibility (tried and failed — does not help; commented out)
    // const powerKMod11 = powerK % M11;
    // const setP = possibleResiduesMod11(P_value % M11, powerKMod11, M);
    // const setQ = possibleResiduesMod11(Q_value % M11, powerKMod11, M);
    // let mask = 0n;
    // for (const a of setP) {
    //     for (const b of setQ) {
    //         const idx = Number(a) * 11 + Number(b);
    //         mask |= 1n << BigInt(idx);
    //     }
    // }
    // const N_mod = N % M11;
    // const rowMask11 = (1n << 11n) - 1n;
    // let mod11Feasible = false;
    // for (let a = 0n; a < M11; a++) {
    //     const rowMask = (mask >> BigInt(Number(a) * 11)) & rowMask11;
    //     if (rowMask === 0n) continue;
    //     if (a === 0n) {
    //         if (N_mod === 0n) mod11Feasible = true;
    //     } else {
    //         const b = (N_mod * mod11Inverse(a)) % M11;
    //         if ((rowMask & (1n << b)) !== 0n) mod11Feasible = true;
    //     }
    //     if (mod11Feasible) break;
    // }
    // if (!mod11Feasible) return 9;

    return 0; // Feasible (not pruned)
}

function checkSolution(branch, N) {
    const p = branch.P_value;
    const q = branch.Q_value;
    const N_big = typeof N === 'bigint' ? N : BigInt(N);
    return p > 1n && q > 1n && p * q === N_big;
}

/**
 * Core IVI work function: finds valid digit pairs for position k.
 * 
 * This function implements the core recurrence relation of the IVI algorithm:
 * 
 *   Σ(i=1 to k) p_i · q_{k-i+1} + c_k = n_k + 10·c_{k+1}
 * 
 * For each possible digit pair (p_k, q_k), it:
 * 1. Computes the sum of all digit products contributing to position k
 * 2. Adds the incoming carry c_k
 * 3. Checks if the result satisfies the IVI constraint
 * 4. Applies global feasibility pruning
 * 5. Returns valid next states with computed carry_out
 * 
 * @param {Object} input - Input state for position k
 * @param {number} input.k - Current position (1-indexed, 1 = LSD)
 * @param {number[]} input.p_history - History of p digits determined so far (LSD-first)
 * @param {number[]} input.q_history - History of q digits determined so far (LSD-first)
 * @param {bigint} input.P_value - Current partial p value as BigInt
 * @param {bigint} input.Q_value - Current partial q value as BigInt
 * @param {number} input.carry_in - Incoming carry c_k (0 for k=1)
 * @param {number[]} input.N_digits - Target number N as digit array (LSD-first)
 * @param {bigint} input.N - Target number N as BigInt
 * @returns {Object[]} Array of valid next states, each containing:
 *   - k: next position (k+1)
 *   - p_history: extended p digit history
 *   - q_history: extended q digit history
 *   - P_value: updated partial p value as BigInt
 *   - Q_value: updated partial q value as BigInt
 *   - carry_in: outgoing carry c_{k+1}
 *   - pk: digit p_k chosen
 *   - qk: digit q_k chosen
 *   - lastTwoDigits: string representation of pk*10 + qk
 */
function workFunction(input) {
    const { k, p_history, q_history, P_value, Q_value, carry_in, N_digits, N, sqrtN } = input;
    
    if (!N_digits || k < 1 || k > N_digits.length) return { states: [], pruningStats: {} };
    if (P_value === undefined || Q_value === undefined || N === undefined || sqrtN === undefined) {
        throw new Error('P_value, Q_value, N (BigInt), and sqrtN (BigInt) are required');
    }
    
    const target_digit = N_digits[k - 1];
    const nextStates = [];
    const isLastDigit = k === N_digits.length;
    const totalDigits = N_digits.length;
    
    // Initialize pruning statistics for this workFunction call
    const pruningStats = {
        1: 0, // Exact Termination Rule
        2: 0, // Immediate Overshoot
        3: 0, // Sqrt-Based Hard Bound
        4: 0, // Explicit Growth Envelope
        5: 0, // Minimum Contribution
        6: 0, // Linear-Term Gap (redundant with 4, but track separately)
        7: 0, // Upper Tail Tightening
        8: 0  // Length Split (reserved)
        // 9: Mod-11 feasibility (tried and failed — does not help)
    };
    
    // #6. Carry Envelope Tightening (Mandatory)
    // Maximum convolution sum is bounded by: maxDigitContribution = 81 * k
    // maxSum = maxDigitContribution + carry_in
    // maxCarryOut = floor(maxSum / 10)
    const maxDigitContribution = 81 * k;
    const maxSum = maxDigitContribution + (carry_in || 0);
    const maxCarryOut = Math.floor(maxSum / 10);
    
    // Pre-compute base sum for terms i=2 to k-1 (terms that don't involve pk or qk)
    // These are: p_2*q_{k-1}, p_3*q_{k-2}, ..., p_{k-1}*q_2
    let baseSum = 0;
    for (let i = 2; i < k; i++) {
        const p_idx = i - 1;  // p_i is at index i-1
        const q_idx = k - i;  // q_{k-i+1} is at index k-i
        if (p_idx < p_history.length && q_idx >= 0 && q_idx < q_history.length) {
            baseSum += multiplyDigits(p_history[p_idx], q_history[q_idx]);
        }
    }
    
    // Pre-compute q_1 (first digit of q, at index 0) for reuse
    const q1 = q_history.length > 0 ? q_history[0] : 0;
    // Pre-compute p_1 (first digit of p, at index 0) for reuse
    const p1 = p_history.length > 0 ? p_history[0] : 0;
    
    // Explore all possible digit pairs (0-9 for each)
    for (let pk = 0; pk <= 9; pk++) {
        for (let qk = 0; qk <= 9; qk++) {
            // Compute sum_{i=1}^{k} p_i * q_{k-i+1}
            // For k=1: only p_1 * q_1 = pk * qk
            // For k>1: 
            //   - i=1: p_1 * q_k = p1 * qk
            //   - i=k: p_k * q_1 = pk * q1
            //   - i=2..k-1: already computed in baseSum
            const sumOfProducts = k === 1 
                ? multiplyDigits(pk, qk)
                : baseSum + multiplyDigits(p1, qk) + multiplyDigits(pk, q1);

            const total = sumOfProducts + carry_in;
            
            // IVI Constraint: total = n_k + 10*c_{k+1}
            // Early pruning: total must be >= target_digit
            if (total < target_digit) {
                continue;
            }
            
            const remainder = total - target_digit;
            // remainder = 10*c_{k+1}, meaning remainder must be >= 0 and divisible by 10
            if (remainder % 10 === 0) {
                const carry_out = remainder / 10;
                // #6. Carry Envelope Tightening: Reject if carry_out > maxCarryOut
                // At the last digit, final carry must be 0
                if (carry_out < 0 || carry_out > maxCarryOut || (isLastDigit && carry_out !== 0)) {
                    continue;
                }
                
                // Create arrays for new state
                const next_p_history = [...p_history, pk];
                const next_q_history = [...q_history, qk];
                
                // #9. Incremental Value Maintenance
                // Update P_value and Q_value incrementally (avoid rebuilding from scratch)
                // k is 1-indexed, so digit at position k goes into 10^(k-1) place
                const powerK = powerOf10(k - 1);
                const new_P_value = P_value + BigInt(pk) * powerK;
                const new_Q_value = Q_value + BigInt(qk) * powerK;
                
                // #1. Structural Rule: All pruning in globalFeasible
                // After adding digit at position k, we have k digits total
                const digitsProcessed = k;  // k is 1-indexed position, equals number of digits after adding
                const newState = {
                    k: digitsProcessed,
                    p_history: next_p_history,
                    q_history: next_q_history
                };
                
                const feasibleResult = globalFeasible(newState, new_P_value, new_Q_value, N, digitsProcessed, totalDigits, sqrtN);
                if (feasibleResult !== 0 && feasibleResult !== true) {
                    // Pruned by step feasibleResult (exclude step 1 final-state pruning from stats)
                    if (feasibleResult !== 1) {
                        pruningStats[feasibleResult] = (pruningStats[feasibleResult] || 0) + 1;
                    }
                    continue;
                }
                if (feasibleResult === false) {
                    // Should not happen with new return values, but handle for safety
                    continue;
                }
                
                const lastTwoDigits = `${pk}${qk}`.padStart(2, '0');
                const next_k = k + 1;
                
                nextStates.push({
                    k: next_k,
                    p_history: next_p_history,
                    q_history: next_q_history,
                    P_value: new_P_value,
                    Q_value: new_Q_value,
                    carry_in: carry_out,
                    pk: pk,
                    qk: qk,
                    lastTwoDigits: lastTwoDigits
                });
            }
        }
    }

    return { states: nextStates, pruningStats: pruningStats };
}

/**
 * Initializes the IVI algorithm state for factorizing N = p × q.
 * 
 * Creates the initial algorithm state with:
 * - Empty digit histories (starting at position k=1, the LSD)
 * - Zero initial carry (c_1 = 0)
 * - Target number N converted to LSD-first digit array
 * 
 * @param {number} p - First prime factor (for reference, not used in computation)
 * @param {number} q - Second prime factor (for reference, not used in computation)
 * @returns {Object} Initial algorithm state containing:
 *   - p, q, N: Original inputs and product
 *   - N_digits: N as digit array (LSD-first)
 *   - frontier: Initial frontier with one empty state at k=1
 *   - step: Current step counter (0 initially)
 *   - history: Empty history array
 */
function initializeAlgorithm(N) {
    const N_big = typeof N === 'bigint' ? N : BigInt(N);
    const N_digits = N_big.toString().split('').reverse().map(Number);
    
    // #5. Square Root Envelope Pruning: Precompute sqrtN once
    const sqrtN = integerSqrt(N_big);
    
    const N_display = N_big <= BigInt(Number.MAX_SAFE_INTEGER) ? Number(N_big) : N_big.toString();
    return {
        p: null,
        q: null,
        N: N_display,
        N_big: N_big,
        sqrtN: sqrtN,
        N_digits: N_digits,
        frontier: [{
            k: 1,
            p_history: [],
            q_history: [],
            P_value: 0n,
            Q_value: 0n,
            carry_in: 0,
            N_digits: N_digits
        }],
        step: 0,
        history: [],
        activeBranches: 1,
        maxActiveBranches: 1,
        // #11. Testing Requirements: Instrumentation
        nodesVisited: 0,
        nodesPruned: 0,
        maxFrontierWidth: 1
    };
}

/**
 * Executes one step of the IVI algorithm.
 * 
 * Processes the current position k by:
 * 1. Applying workFunction to all branches in the frontier
 * 2. Collecting valid next states with parent relationships
 * 3. Checking for solution completion (when k equals or exceeds N_digits.length)
 * 4. Updating state with new frontier and history
 * 
 * The algorithm terminates when:
 * - A valid solution is found (success: true)
 * - No valid branches exist (done: true)
 * - All digits processed without solution (done: true)
 * 
 * @param {Object} state - Current algorithm state
 * @param {number} state.step - Current step counter
 * @param {number} state.N - Target number to factorize
 * @param {number[]} state.N_digits - N as digit array (LSD-first)
 * @param {Object[]} state.frontier - Current frontier of active branches
 * @param {Object[]} state.history - History of all processed steps
 * @returns {Object} Updated state with one of:
 *   - success: true, foundP, foundQ, solutionPath (if solution found)
 *   - done: true (if no solution possible)
 *   - Updated step, frontier, and history (if continuing)
 */
function stepAlgorithm(state) {
    if (state.maxSteps != null && state.step >= state.maxSteps) {
        return { ...state, done: true };
    }
    const currentK = state.step + 1;
    
    // If we've already processed all digits, terminate
    if (currentK > state.N_digits.length) {
        return { ...state, done: true };
    }

    const target_digit = state.N_digits[currentK - 1];
    
    // Process all branches in frontier, tracking parent relationships
    const allResults = [];
    let nodesVisited = state.nodesVisited || 0;
    let nodesPruned = state.nodesPruned || 0;
    
    // Aggregate pruning statistics
    const stepPruningStats = {
        1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0
    };
    
    state.frontier.forEach((branch, parentIdx) => {
        const workResult = workFunction({ 
            ...branch, 
            k: currentK, 
            N_digits: state.N_digits,
            N: state.N_big || BigInt(state.N),
            sqrtN: state.sqrtN
        });
        
        const candidates = workResult.states || [];
        const pruningStats = workResult.pruningStats || {};
        
        // Aggregate pruning statistics
        for (let step = 1; step <= 8; step++) {
            stepPruningStats[step] = (stepPruningStats[step] || 0) + (pruningStats[step] || 0);
        }
        
        // #11. Instrumentation: Track visited and pruned nodes
        // Each digit pair (pk, qk) is a candidate node
        // We visit 100 candidates per branch (10*10), but only some pass pruning
        nodesVisited += 100; // 10*10 digit pairs per branch
        nodesPruned += (100 - candidates.length);
        
        candidates.forEach(result => allResults.push({ ...result, parentIdx }));
    });
    
    // Cap frontier size to avoid UI freeze when branching explodes (e.g. small N)
    if (state.maxFrontierSize != null && allResults.length > state.maxFrontierSize) {
        allResults = allResults.slice(0, state.maxFrontierSize);
    }
    
    // Update cumulative pruning statistics
    const updatedPruningStats = { ...state.pruningStats };
    for (let step = 1; step <= 8; step++) {
        updatedPruningStats[step] = (updatedPruningStats[step] || 0) + stepPruningStats[step];
    }

    // If no valid branches found, terminate
    if (allResults.length === 0) {
        return { 
            ...state, 
            done: true,
            activeBranches: 0,
            maxActiveBranches: state.maxActiveBranches || 0,
            nodesVisited: nodesVisited,
            nodesPruned: nodesPruned,
            maxFrontierWidth: state.maxFrontierWidth || 0,
            pruningStats: updatedPruningStats
        };
    }

    // Check for solution when processing the last digit
    if (currentK === state.N_digits.length) {
        for (let branchIdx = 0; branchIdx < allResults.length; branchIdx++) {
            const branch = allResults[branchIdx];
            if (branch.carry_in === 0 && checkSolution(branch, state.N)) {
                let foundP = branch.P_value;
                let foundQ = branch.Q_value;
                if (foundP > foundQ) [foundP, foundQ] = [foundQ, foundP];
                const activeBranches = allResults.length;
                const maxActiveBranches = Math.max(state.maxActiveBranches || 0, activeBranches);
                const maxFrontierWidth = Math.max(state.maxFrontierWidth || 0, activeBranches);
                
                return {
                    ...state,
                    step: currentK,
                    frontier: allResults,
                    history: [...state.history, {
                        k: currentK,
                        target_digit: target_digit,
                        branches: mapBranchesToHistory(allResults, branchIdx)
                    }],
                    success: true,
                    foundP: foundP.toString(),
                    foundQ: foundQ.toString(),
                    solutionPath: buildSolutionPath(branch),
                    activeBranches: activeBranches,
                    maxActiveBranches: maxActiveBranches,
                    nodesVisited: nodesVisited,
                    nodesPruned: nodesPruned,
                    maxFrontierWidth: maxFrontierWidth,
                    pruningStats: updatedPruningStats
                };
            }
        }
        // No solution found after processing last digit
        const activeBranches = allResults.length;
        const maxActiveBranches = Math.max(state.maxActiveBranches || 0, activeBranches);
        const maxFrontierWidth = Math.max(state.maxFrontierWidth || 0, activeBranches);
        return { 
            ...state, 
            done: true,
            activeBranches: activeBranches,
            maxActiveBranches: maxActiveBranches,
            nodesVisited: nodesVisited,
            nodesPruned: nodesPruned,
            maxFrontierWidth: maxFrontierWidth,
            pruningStats: updatedPruningStats
        };
    }

    // Store step history and continue
    const stepHistory = {
        k: currentK,
        target_digit: target_digit,
        branches: mapBranchesToHistory(allResults)
    };

    const activeBranches = allResults.length;
    const maxActiveBranches = Math.max(state.maxActiveBranches || 0, activeBranches);
    const maxFrontierWidth = Math.max(state.maxFrontierWidth || 0, activeBranches);

    return {
        ...state,
        step: currentK,
        frontier: allResults,
        history: [...state.history, stepHistory],
        activeBranches: activeBranches,
        maxActiveBranches: maxActiveBranches,
        nodesVisited: nodesVisited,
        nodesPruned: nodesPruned,
        maxFrontierWidth: maxFrontierWidth,
        pruningStats: updatedPruningStats
    };
}

    // Export functions to window for browser use
    if (typeof window !== 'undefined') {
        window.initializeAlgorithm = initializeAlgorithm;
        window.stepAlgorithm = stepAlgorithm;
        window.workFunction = workFunction; // Export for DCP workers
    }
})();