prime sieve formatted by symmetric zero-padding

prime-zeroes9.py

#!/usr/bin/env python3
import math
import time
from itertools import product
import multiprocessing as mp
from collections import defaultdict, Counter
import random

def is_probable_prime(n, rounds=20):
    """
    Miller-Rabin primality test for large integers.
    Returns True if n is probably prime, False otherwise.
    
    This implementation is optimized for deterministic exploration.
    """
    if n < 2:
        return False
    
    # Quick checks for small primes and their multiples
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    for sp in small_primes:
        if n == sp:
            return True
        if n % sp == 0:
            return False
    
    # Factor out powers of 2 from n-1
    d = n - 1
    r = 0
    while d % 2 == 0:
        d //= 2
        r += 1
    
    # Miller-Rabin test
    for _ in range(rounds):
        a = random.randrange(2, n - 1)
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for __ in range(r - 1):
            x = (x * x) % n
            if x == n - 1:
                break
        else:
            return False
    
    return True

def build_double_membrane(outer_bound, k_outer, inner_bound, k_inner, middle=""):
    """
    Builds a double-membrane structured number of the form:
    outer_bound + (k_outer zeros) + inner_bound + (k_inner zeros) + 
    middle + (k_inner zeros) + inner_bound + (k_outer zeros) + outer_bound
    
    Returns it as a string.
    """
    outer_zeros = "0" * k_outer
    inner_zeros = "0" * k_inner
    
    return f"{outer_bound}{outer_zeros}{inner_bound}{inner_zeros}{middle}{inner_zeros}{inner_bound}{outer_zeros}{outer_bound}"

def generate_all_middle_values(max_length):
    """
    Generates all possible middle values up to the specified maximum length.
    Returns a dictionary mapping length to a list of all possible middle values of that length.
    """
    all_middles = {0: [""]}  # Empty middle
    
    for length in range(1, max_length + 1):
        middles = []
        # First digit non-zero (1-9)
        first_digits = range(1, 10)
        if length > 1:
            # Remaining digits (0-9)
            remaining_combinations = product(range(10), repeat=length-1)
            for first in first_digits:
                for combo in remaining_combinations:
                    middle = str(first) + ''.join(str(d) for d in combo)
                    middles.append(middle)
        else:
            # Single-digit case
            middles = [str(d) for d in first_digits]
        all_middles[length] = middles
    
    return all_middles

def test_configuration(params):
    """
    Tests a specific double-membrane configuration for primality.
    Uses deterministic exhaustive search over all possible middle values.
    
    params: (outer_bound, k_outer, inner_bound, k_inner, max_middle_length)
    
    Returns a list of tuples (number, middle) for discovered primes.
    """
    outer_bound, k_outer, inner_bound, k_inner, max_middle_length = params
    
    # Generate all possible middle values
    all_middles = generate_all_middle_values(max_middle_length)
    
    # Test each middle value
    results = []
    
    for mid_len in range(max_middle_length + 1):
        for mid_str in all_middles[mid_len]:
            # Build the double-membrane number
            num_str = build_double_membrane(outer_bound, k_outer, inner_bound, k_inner, mid_str)
            num = int(num_str)
            
            # Test for primality
            if is_probable_prime(num):
                results.append((num, mid_str))
    
    return results, len(all_middles[0]) + sum(len(all_middles[i]) for i in range(1, max_middle_length + 1))

def worker_task(params):
    """Worker task for parallel processing."""
    try:
        o_bd, o_k, i_bd, i_k, max_mid_len = params
        results, tested = test_configuration(params)
        return results, tested, params
    except Exception as e:
        return [], 0, params, str(e)

def format_number(num):
    """Format large numbers with commas for readability."""
    return f"{num:,}"

def calculate_total_configurations(outer_bounds, outer_ks, inner_bounds, inner_ks, max_middle_length):
    """Calculate the total number of configurations that will be tested."""
    total_configs = len(outer_bounds) * len(outer_ks) * len(inner_bounds) * len(inner_ks)
    
    # Calculate total middle values
    total_tests = 0
    for mid_len in range(max_middle_length + 1):
        if mid_len == 0:
            # Empty middle
            middle_count = 1
        else:
            # First digit can be 1-9, rest can be 0-9
            middle_count = 9 * (10 ** (mid_len - 1))
        total_tests += middle_count
    
    return total_configs, total_tests * total_configs

def deterministic_search(outer_bounds, outer_ks, inner_bounds, inner_ks, max_middle_length, max_processes=None):
    """
    Performs a deterministic search for double-membrane primes using specified parameters.
    Uses parallel processing for efficiency.
    
    Returns a comprehensive analysis of discovered primes.
    """
    # Calculate total configurations for progress reporting
    total_configs, total_tests = calculate_total_configurations(
        outer_bounds, outer_ks, inner_bounds, inner_ks, max_middle_length
    )
    
    print(f"Deterministic Double-Membrane Prime Search")
    print(f"------------------------------------------")
    print(f"Parameters:")
    print(f"  Outer bounding digits: {outer_bounds}")
    print(f"  Outer k values: {outer_ks}")
    print(f"  Inner bounding digits: {inner_bounds}")
    print(f"  Inner k values: {inner_ks}")
    print(f"  Maximum middle length: {max_middle_length}")
    print(f"Search space:")
    print(f"  Total configurations: {format_number(total_configs)}")
    print(f"  Total tests: {format_number(total_tests)}")
    
    # Generate all parameter combinations
    all_params = [
        (o_bd, o_k, i_bd, i_k, max_middle_length)
        for o_bd in outer_bounds
        for o_k in outer_ks
        for i_bd in inner_bounds
        for i_k in inner_ks
    ]
    
    # Set up multiprocessing
    if max_processes is None:
        max_processes = mp.cpu_count()
    print(f"Using {max_processes} CPU cores for parallel processing")
    
    # Store results in a nested dictionary
    results = {
        o_bd: {
            o_k: {
                i_bd: {
                    i_k: []
                    for i_k in inner_ks
                } for i_bd in inner_bounds
            } for o_k in outer_ks
        } for o_bd in outer_bounds
    }
    
    # Statistics tracking
    total_primes = 0
    total_tested = 0
    config_tested = 0
    start_time = time.time()
    
    # Process configurations in parallel
    with mp.Pool(processes=max_processes) as pool:
        for batch_results in pool.imap_unordered(worker_task, all_params, chunksize=1):
            if len(batch_results) == 3:
                primes, tested, params = batch_results
                error = None
            else:
                primes, tested, params, error = batch_results
            
            o_bd, o_k, i_bd, i_k, _ = params
            
            # Store results
            results[o_bd][o_k][i_bd][i_k] = primes
            
            # Update statistics
            config_tested += 1
            total_tested += tested
            total_primes += len(primes)
            
            # Progress reporting
            elapsed = time.time() - start_time
            completed_pct = (config_tested / total_configs) * 100
            rate = total_tested / elapsed if elapsed > 0 else 0
            eta = (total_tests - total_tested) / rate if rate > 0 else 0
            
            if error:
                print(f"Error in configuration {params}: {error}")
            
            if config_tested % 10 == 0 or config_tested == total_configs:
                print(f"Progress: {config_tested}/{total_configs} configurations "
                      f"({completed_pct:.2f}%) | "
                      f"Found: {total_primes} primes | "
                      f"Rate: {format_number(int(rate))}/sec | "
                      f"ETA: {int(eta/60)} minutes")
    
    # Final statistics
    end_time = time.time()
    total_time = end_time - start_time
    
    print(f"\nSearch completed in {total_time:.2f} seconds")
    print(f"Tested {format_number(total_tested)} numbers, found {total_primes} primes")
    print(f"Prime density: {total_primes / total_tested * 100:.6f}%")
    
    return results, {
        'total_primes': total_primes,
        'total_tested': total_tested,
        'total_time': total_time,
        'prime_density': total_primes / total_tested
    }

def analyze_results(results, stats):
    """
    Performs a comprehensive analysis of the discovered primes.
    
    Returns analysis metrics and insights.
    """
    # Extract all primes into a flat list for analysis
    all_primes = []
    for o_bd in results:
        for o_k in results[o_bd]:
            for i_bd in results[o_bd][o_k]:
                for i_k in results[o_bd][o_k][i_bd]:
                    for prime, middle in results[o_bd][o_k][i_bd][i_k]:
                        all_primes.append({
                            'prime': prime,
                            'outer_bd': o_bd,
                            'outer_k': o_k,
                            'inner_bd': i_bd,
                            'inner_k': i_k,
                            'middle': middle,
                            'digits': len(str(prime))
                        })
    
    # Count by parameter
    outer_bd_count = Counter(p['outer_bd'] for p in all_primes)
    outer_k_count = Counter(p['outer_k'] for p in all_primes)
    inner_bd_count = Counter(p['inner_bd'] for p in all_primes)
    inner_k_count = Counter(p['inner_k'] for p in all_primes)
    middle_count = Counter(p['middle'] for p in all_primes)
    middle_len_count = Counter(len(p['middle']) for p in all_primes)
    digit_count = Counter(p['digits'] for p in all_primes)
    
    # Top middle values
    top_middles = middle_count.most_common(20)
    
    # Parameter combinations analysis
    param_combo_count = Counter((p['outer_bd'], p['outer_k'], p['inner_bd'], p['inner_k']) for p in all_primes)
    top_combos = param_combo_count.most_common(20)
    
    # Density by parameter (requires total tests by parameter)
    # For simplicity, we'll just use counts for comparative analysis
    
    return {
        'total_unique_primes': len(all_primes),
        'by_outer_bd': dict(outer_bd_count),
        'by_outer_k': dict(outer_k_count),
        'by_inner_bd': dict(inner_bd_count),
        'by_inner_k': dict(inner_k_count),
        'by_middle_len': dict(middle_len_count),
        'by_digits': dict(digit_count),
        'top_middles': top_middles,
        'top_combos': top_combos,
        'all_primes': all_primes
    }

def print_histogram(data, title, sort_by_key=True):
    """Prints a simple text histogram of frequency data."""
    print(f"\n{title}:")
    
    # Sort data for presentation
    if sort_by_key:
        sorted_data = sorted(data.items())
    else:
        sorted_data = sorted(data.items(), key=lambda x: x[1], reverse=True)
    
    # Find maximum value for scaling
    max_value = max(data.values())
    scale_factor = 40.0 / max_value if max_value > 0 else 1
    
    # Print histogram
    for key, value in sorted_data:
        bar_length = int(value * scale_factor)
        bar = '█' * bar_length
        print(f"  {key}: {value:5d} {bar}")

def print_top_items(items, title, format_key=None):
    """Prints top items with formatting."""
    print(f"\n{title}:")
    for i, (key, count) in enumerate(items, 1):
        key_str = format_key(key) if format_key else key
        print(f"  {i:2d}. {key_str}: {count}")

def format_combo(combo):
    """Format parameter combination for display."""
    o_bd, o_k, i_bd, i_k = combo
    return f"outer={o_bd}, outer_k={o_k}, inner={i_bd}, inner_k={i_k}"

def print_results(results):
    """Prints the results in a structured, readable format."""
    # ANSI color codes for visual clarity
    bold_yellow = "\033[1;33m"
    bold_green = "\033[1;32m"
    reset = "\033[0m"
    
    # Column widths for alignment
    col_widths = [35, 10, 8, 8, 10, 8, 11]
    
    # Header labels
    headers = ["Prime", "Middle", "Digits", "Outer BD", "Outer K", "Inner BD", "Inner K"]
    
    # Print results organized by outer bounding digit and outer k
    for o_bd in sorted(results.keys()):
        print(f"\n{bold_yellow}=== Outer Bounding Digit = '{o_bd}' ==={reset}")
        
        for o_k in sorted(results[o_bd].keys()):
            print(f"  -- Outer k = {o_k} --")
            
            # Print header row
            header_row = "  "
            for i, header in enumerate(headers):
                header_row += f"{header:<{col_widths[i]}}"
            print(header_row)
            print("  " + "-" * (sum(col_widths)))
            
            # Track if we found any primes for this outer_bd/outer_k combination
            found_primes = False
            total_primes = 0
            
            # Print all primes for this outer bounding digit and k value
            for i_bd in sorted(results[o_bd][o_k].keys()):
                for i_k in sorted(results[o_bd][o_k][i_bd].keys()):
                    primes_list = results[o_bd][o_k][i_bd][i_k]
                    
                    if not primes_list:
                        continue
                    
                    found_primes = True
                    total_primes += len(primes_list)
                    
                    # Sort primes in ascending order
                    for prime_val, mid_str in sorted(primes_list, key=lambda x: x[0]):
                        dcount = len(str(prime_val))
                        
                        prime_colored = f"{bold_green}{prime_val}{reset}"
                        
                        print(
                            "  "
                            f"{prime_colored:<{col_widths[0]}}"
                            f"{mid_str:<{col_widths[1]}}"
                            f"{str(dcount):<{col_widths[2]}}"
                            f"{o_bd:<{col_widths[3]}}"
                            f"{str(o_k):<{col_widths[4]}}"
                            f"{i_bd:<{col_widths[5]}}"
                            f"{str(i_k):<{col_widths[6]}}"
                        )
            
            if not found_primes:
                print("  No primes found for this configuration.")
            else:
                print(f"  Total primes for this configuration: {total_primes}")
            
            print("")

def print_analysis(analysis):
    """Prints a comprehensive analysis of the search results."""
    print("\n==========================================")
    print("Double-Membrane Prime Number Analysis")
    print("==========================================")
    
    print(f"\nTotal unique primes found: {analysis['total_unique_primes']}")
    
    # Print histograms for parameter frequencies
    print_histogram(analysis['by_outer_bd'], "Primes by Outer Bounding Digit")
    print_histogram(analysis['by_outer_k'], "Primes by Outer K Value")
    print_histogram(analysis['by_inner_bd'], "Primes by Inner Bounding Digit")
    print_histogram(analysis['by_inner_k'], "Primes by Inner K Value")
    print_histogram(analysis['by_middle_len'], "Primes by Middle Length")
    print_histogram(analysis['by_digits'], "Primes by Total Digits")
    
    # Print top items
    print_top_items(analysis['top_middles'], "Top 20 Most Common Middle Values")
    print_top_items(analysis['top_combos'], "Top 20 Most Productive Parameter Combinations", format_combo)
    
    # Print some example primes
    print("\nExample Primes from Different Configurations:")
    samples = {}
    for prime_info in analysis['all_primes']:
        key = (prime_info['outer_bd'], prime_info['outer_k'], prime_info['inner_bd'], prime_info['inner_k'])
        if key not in samples and len(samples) < 15:
            samples[key] = prime_info
    
    for i, (_, prime_info) in enumerate(sorted(samples.items()), 1):
        print(f"  {i:2d}. {prime_info['prime']} with outer={prime_info['outer_bd']}, "
              f"outer_k={prime_info['outer_k']}, inner={prime_info['inner_bd']}, "
              f"inner_k={prime_info['inner_k']}, middle='{prime_info['middle']}'")

def main():
    """Main function to execute the deterministic double-membrane prime search."""
    print("Deterministic Double-Membrane Prime Number Search")
    print("------------------------------------------------")
    
    # Define search parameters
    outer_bounds = ['1', '3', '7', '9']  # Prime or coprime to 10
    outer_ks = [2, 3, 5, 7]  # Small prime values
    inner_bounds = ['1', '3', '5', '7', '9']  # Odd digits
    inner_ks = [2, 3, 5, 7]  # Small prime values
    max_middle_length = 2  # Maximum length of middle values
    
    # Use all available CPU cores
    max_processes = mp.cpu_count()
    
    # Run the deterministic search
    results, stats = deterministic_search(
        outer_bounds, 
        outer_ks, 
        inner_bounds, 
        inner_ks, 
        max_middle_length,
        max_processes
    )
    
    # Analyze results
    analysis = analyze_results(results, stats)
    
    # Print analysis
    print_analysis(analysis)
    
    # Print detailed results
    print("\nDetailed Prime Listings:")
    print_results(results)

if __name__ == "__main__":
    # Set random seed for reproducible primality testing
    random.seed(42)
    main()