Prime-Square-Sum

📓 View the Mathematica Notebook (Requires Mathematica or Wolfram Player)

The python program is a computational platform for validating mathematical relationships. It checks LHS expressions against RHS targets and stores sequences for future analysis. The main focus of prime-square-sum is on summing triangular row values together (where stf() below is an acronym for "sum triangular factors") to determine if:

is equal to the series of squared/cubed primes.

b = triangular number (also the number base);              //equal to: tri(r)=(r^2+r)/2
r = qg(b) = size of the base row of the triangular number; //qg(b) = 1/2(-1+sqrt(1+8b)
z = row in the triangular number;  //ex. tf(10,4)=0123; tf(10,3)=456; tf(10,2)=78, etc.)

Where tf(), the individual "triangular row factors", are defined to be:

There is an interesting relationship when {b=10, r=4} where the sum of the triangular rows in base-10, 0123 + 456 + 78 + 9, happens to work out to be the sum of the first seven squared primes.

stf(10) = 2² + 3² + 5² + 7² + 11² + 13² + 17² = 666

What I find fascinating about this relationship is the resultant value 666 is a triangular number itself. So the question then is if we were able to sum the rows of a 666 element triangle with 36 rows in base-666 would the result also be the sum of squared or cubed primes?

This program attempts to provide an answer. The base-10 number from stf(666) unfortunately though is massively large at 98 digits:

37005443752611483714216385166550857181329086284892731078593232926279977894581784762614450464857290

While stf(b) is a closed-form function, primesum(n,p) is computationally irreducible. This means there is no readily available formula to directly find which n satisfies primesum(n,p) = stf(b). So, in the absence of a better method, we enumerate prime sums and check for matches. See docs/rationale.md for additional details as to why this tool needs to be as flexible as it is.

Due to the huge size of the 98-digits I've further adapted prime-square-sum to work with multiprocessing and CUDA (via CuPy) to speed up the computations. See the Mathematica notebook in paper and notes/ for more details.

Usage

# Find n where sum of first n squared primes equals 666 (this is the known result: stf(10))
python prime-square-sum.py --expr "does_exist primesum(n,2) == 666"
# Output: Found: n=7  (because 2² + 3² + 5² + 7² + 11² + 13² + 17² = 666)

# Shorthand using default expression
python prime-square-sum.py --target 666

The --target flag searches against the default expression primesum(n,2) (sum of squared primes). This default is defined in equations.json and can be customized via config.json. See docs/equations.md for details.

# Find matches between prime sums and triangular numbers
python prime-square-sum.py --expr "for_any primesum(n,2) == tri(m)" --max n:100 --max m:50

# List available functions
python prime-square-sum.py --list functions

The Open Question

The core investigation: does the sum of triangular row factors as stf(666), the 98-digit number, equal a sum of squared, cubed, or some p-th power primes?

# The main search — does stf(666) = primesum(n, 2) for some n?
python prime-square-sum.py --expr "does_exist primesum(n,2) == trisum(666)" \
    --algorithm sieve:primesieve --max n:1000000000

# The pattern may continue at a higher power — search p = 3, 4, 5
python prime-square-sum.py --expr "does_exist primesum(n,3) == trisum(666)" \
    --algorithm sieve:primesieve --max n:50000000
python prime-square-sum.py --expr "does_exist primesum(n,4) == trisum(666)" \
    --algorithm sieve:primesieve --max n:10000000
python prime-square-sum.py --expr "does_exist primesum(n,5) == trisum(666)" \
    --algorithm sieve:primesieve --max n:5000000

# Skip early values that are too small to match a 98-digit target
python prime-square-sum.py --expr "does_exist primesum(n,p) == trisum(666)" \
    --algorithm sieve:primesieve -min n:50000000 -max n:1000000000 --max p:5

This search requires billions of prime sums and benefits from --algorithm sieve:primesieve (C++ primesieve library via conda). See Installation for setup. The answer to this question would establish whether the chain primesum(3,1)=10 → stf(primesum(3,1))=666 → primesum(7,2)=666 → stf(stf(primesum(3,1))) → primesum(?,?), or perhaps easier to see primesum(3,1) → stf(10) → primesum(7,2) → stf(666) → primesum(?,?), continues — and whether it may continue indefinitely. The formal proofs provide the algebraic framework (closed-form stf(b), verified in Lean 4), while the numerical search here would supply the concrete witness needed for induction.

Arithmetic Expressions (v0.7.12+)

Use standard arithmetic operators directly in expressions:

# Arithmetic in expressions
python prime-square-sum.py --expr "does_exist n**2 == 25"                     # n=5
python prime-square-sum.py --expr "does_exist tri(n) + 1 == 11" --max n:10    # n=4
python prime-square-sum.py --expr "verify (2 + 3) * 4 == 20"                  # true

# Operators: +  -  *  /  //  %  **  ^  (unary: -x, +x)
# Python-compatible precedence and associativity

Boolean & Bitwise Operators (v0.7.14+)

# Boolean logic with short-circuit evaluation
python prime-square-sum.py --expr "does_exist n > 0 and n < 10 and tri(n) == 28" --max n:20

# Bitwise operations via keyword operators
python prime-square-sum.py --expr "verify 5 xor 3 == 6"
python prime-square-sum.py --expr "verify 5 band 3 == 1"

# Chained comparisons
python prime-square-sum.py --expr "does_exist 1 < n < 10 and tri(n) == 28" --max n:20

# Context blocks for operator disambiguation
python prime-square-sum.py --expr "verify bit[2^3] == 1"    # ^ is XOR in bit context
python prime-square-sum.py --expr "verify 2^3 == 8"          # ^ is power by default

See docs/expressions.md for the precedence table and behavioral notes.

Solve & Enumerate (v0.8.0+)

Use solve as a calculator or to enumerate sequence values:

# Calculator mode — evaluate an expression directly
python prime-square-sum.py --expr "solve tri(36)"                    # 666
python prime-square-sum.py --expr "tri(36)"                          # 666 (implicit)

# Enumerate sequence values
python prime-square-sum.py --expr "solve tri(n)" --max n:10          # n=1: 1, n=2: 3, ...
python prime-square-sum.py --expr "for_any primesum(n,2)" --max n:7  # tabulate squared prime sums

Verify Mode (v0.7.6+)

Verify known results without iteration using the verify directive:

# Explicit verify - returns true/false
python prime-square-sum.py --expr "verify primesum(7,2) == 666"  # Returns: true

# Implicit verify - auto-detected when no free variables
python prime-square-sum.py --expr "primesum(7,2) == trisum(10)"  # Returns: true

# Using saved equation
python prime-square-sum.py --equation verify-stf10

# JSON output for programmatic use
python prime-square-sum.py --expr "verify primesum(7,2) == trisum(10)" --format json
# Returns: {"verified": true}

The verify directive evaluates closed formulas (no free variables) directly, returning true or false without any iteration.

Algorithm Selection (v0.7.5+)

Control which sieve algorithm is used for prime generation:

# List available algorithms
python prime-square-sum.py --list algorithms

# Force segmented sieve (bounded memory)
python prime-square-sum.py --target 666 --algorithm sieve:segmented

# Use minimal memory mode
python prime-square-sum.py --target trisum(10) --prefer minimal

# Set memory limit
python prime-square-sum.py --target trisum(10) --max-memory 512

Available algorithms: auto, primesieve, basic, segmented, individual

Settings can also be saved in config.json (see config.json.example).

Installation

Recommended: Conda

Use the provided environment.yml for easiest setup:

# Create environment from file
conda env create -f environment.yml
conda activate prime-square-sum

# Verify
python prime-square-sum.py --rhs 666

Or install manually:

conda create -n prime-square-sum python=3.10
conda activate prime-square-sum
conda install -c conda-forge python-primesieve numpy pytest cupy lark

Note: Windows requires Python 3.10 for primesieve (conda-forge builds only available up to 3.10).

This installs:

• python-primesieve - Fast prime generation (~125M primes/sec)
• numpy - Array operations
• cupy - GPU acceleration (requires NVIDIA GPU + CUDA)
• lark - Expression parser for --expr queries (v0.7.1+)

Other Installation Methods

For venv, pip, Docker, GPU setup, and detailed troubleshooting, see docs/install.md.

Documentation

Verify GPU Setup

python -c "from utils.gpu import init_gpu, print_gpu_status; init_gpu(); print_gpu_status()"

Precomputed Primes

Large prime datasets available via GitHub Releases:

Project Structure

Prime-Square-Sum/
├── prime-square-sum.py      # Main program
├── run_tests.py             # Test runner
├── utils/                   # Sieve, I/O, GPU utilities
├── tests/                   # Unit tests
├── proofs/                  # Lean 4 + Mathlib formal proofs
├── verification/            # Verification scripts
├── paper and notes/         # Mathematica notebook
└── data/                    # Prime files (download from Releases)
    ├── npy_files/           # NumPy format (recommended)
    └── dat_files/           # Pickle format (legacy)

Roadmap

See ROADMAP.md for version history and planned features.

Current: v0.8.x — sequence enumeration, solve directive, iterator improvements, CLI cleanup.

Related

See Zero_AG to The Scarcity Framework for related theoretical work.

Contributing

Pull requests welcome. Start a conversation in Discussions.

License

prime-square-sum, Copyright (C) 2010 to 2026 Dustin Darcy. This work is licensed under CC BY-NC-ND 4.0 (Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International).

• Attribution: You must give appropriate credit when referencing this work
• NonCommercial: You may not use this for commercial purposes
• NoDerivatives: You may not distribute modified versions without permission

See LICENSE for full details.