Skip to content

Instantly share code, notes, and snippets.

@Hermann-SW
Created March 12, 2026 20:42
Show Gist options
  • Select an option

  • Save Hermann-SW/b538c991c922b1e8fc06abbab554b3ef to your computer and use it in GitHub Desktop.

Select an option

Save Hermann-SW/b538c991c922b1e8fc06abbab554b3ef to your computer and use it in GitHub Desktop.
SageMath diophantine S-Unit solve example: error free after many iterations of Google Gemini; changed to ℚ and added generator output by me
# 1. Setup
x = polygen(QQ, 'x')
# K.<i> = NumberField(x^2 + 1) # ℚ(i)
K.<i> = NumberField(x - 1)
# Using this, the root a is just 1. This forces Sage to wrap the rational
# numbers in a "NumberField object" which possesses the .S_unit_group() method.
#
S_list = K.primes_above(2) + K.primes_above(3)
S = tuple(S_list)
# 2. Get the actual generators of the S-unit group
U_S = K.S_unit_group(S=S)
gens = U_S.gens()
print("Generators:")
for g in U_S.gens_values():
print(g)
# 3. Solve
from sage.rings.number_field.S_unit_solver import solve_S_unit_equation
solutions = solve_S_unit_equation(K, S)
def vector_to_number(vec, generators):
"""Manually multiply generators by their exponents."""
res = K(1)
for i in range(len(vec)):
# We convert the generator to the field K to be safe
res *= K(generators[i])**vec[i]
return res
print(f"Found {len(solutions)} solutions:")
for entry in solutions:
# entry[0] and entry[1] are the exponent vectors
x_val = vector_to_number(entry[0], gens)
y_val = vector_to_number(entry[1], gens)
print(f"{x_val} + {y_val} == {x_val + y_val}")
print("-" * 20)
@Hermann-SW
Copy link
Author

Hermann-SW commented Mar 13, 2026

Instead S={2,3} using S={2,3,5,7} ...

hermann@7950x:~$ diff S-Unit.sol.sage S-Unit.sol2.sage 
9c9
< S_list = K.primes_above(2) + K.primes_above(3)
---
> S_list = K.primes_above(2) + K.primes_above(3) + K.primes_above(5) + K.primes_above(7)
hermann@7950x:~$ 

... leads to runtime explosion: aborted after 22.5h single core on AMD 7950X CPU with boost frequency of 5.4GHz:

hermann@7950x:~$ time Downloads/SageMath-10.8-x86_64.AppImage S-Unit.sol2.sage 
Generators:
-1
2
3
5
7
^CTraceback (most recent call last):
  File "/tmp/.mount_SageMaebpACL/app_root/local/bin/sage", line 5, in <module>
    from sage.cli.__main__ import main
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/cli/__main__.py", line 9, in <module>
    sys.exit(main())
             ~~~~^^
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/cli/__init__.py", line 56, in main
    return RunFileCmd(options).run()
           ~~~~~~~~~~~~~~~~~~~~~~~^^
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/cli/run_file_cmd.py", line 50, in run
    eval(compile(open(input_file, 'rb').read(), input_file, 'exec'), sage_globals())
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/tmp/sage_52ap2op7/S-Unit.sol271_34m73.sage.py", line 27, in <module>
    solutions = solve_S_unit_equation(K, S)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2835, in solve_S_unit_equation
    S_unit_solutions = sieve_below_bound(K, list(S), final_LLL_bound, verbose=verbose)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2716, in sieve_below_bound
    complement_exp_vec_dict = construct_complement_dictionaries(split_primes_list, SUK, verbose=verbose)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2343, in construct_complement_dictionaries
    if drop_vector(exp_vec, q, p, comp_exp_vec):
       ~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2130, in drop_vector
    if compatible_cv in complement_ev_dict[q][compatible_exp_vec]:
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "cysignals/signals.pyx", line 357, in cysignals.signals.python_check_interrupt
KeyboardInterrupt

real	1348m44.768s
user	1348m40.494s
sys	0m0.309s

hermann@7950x:~$ 

@Hermann-SW
Copy link
Author

Now tried with S={2,3,5}, hoping for better runtime than with {2,3,5,7} given that S={2,3} takes less than 5 minutes.

I aborted after 13.5h of 5.4GHz boost frequency run on AMD 7950X CPU:

hermann@7950x:~$ time Downloads/SageMath-10.8-x86_64.AppImage S-Unit.sol3.sage 
Generators:
-1
2
3
5
^CTraceback (most recent call last):
  File "/tmp/.mount_SageMaKhnGdE/app_root/local/bin/sage", line 5, in <module>
    from sage.cli.__main__ import main
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/cli/__main__.py", line 9, in <module>
    sys.exit(main())
             ~~~~^^
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/cli/__init__.py", line 56, in main
    return RunFileCmd(options).run()
           ~~~~~~~~~~~~~~~~~~~~~~~^^
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/cli/run_file_cmd.py", line 50, in run
    eval(compile(open(input_file, 'rb').read(), input_file, 'exec'), sage_globals())
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/tmp/sage_zm6ra6cv/S-Unit.sol3bjhh8xix.sage.py", line 27, in <module>
    solutions = solve_S_unit_equation(K, S)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2835, in solve_S_unit_equation
    S_unit_solutions = sieve_below_bound(K, list(S), final_LLL_bound, verbose=verbose)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2718, in sieve_below_bound
    cs_list = compatible_systems(split_primes_list, complement_exp_vec_dict)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2489, in compatible_systems
    old_systems = compatible_systems(S1, complement_exp_vec_dict)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2489, in compatible_systems
    old_systems = compatible_systems(S1, complement_exp_vec_dict)
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2496, in compatible_systems
    if all((compatible_vectors_check(exp_vec, exp_vec_qj, g, l) and
       ~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
            compatible_vectors_check(comp_vec, comp_vec_qj, g, l))
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
           for g, (exp_vec_qj, comp_vec_qj) in zip(gcds, old_system)):
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/var/tmp/sage-10.8/local/lib/python3.13/site-packages/sage/rings/number_field/S_unit_solver.py", line 2496, in <genexpr>
    if all((compatible_vectors_check(exp_vec, exp_vec_qj, g, l) and
    
  File "cysignals/signals.pyx", line 357, in cysignals.signals.python_check_interrupt
KeyboardInterrupt
^C
real	831m40.662s
user	831m35.096s
sys	0m2.705s

hermann@7950x:~$ 

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment