0
# Number Theory
1

2
Comprehensive number theory functions including prime generation and testing, integer factorization, modular arithmetic, and mathematical utility functions. This module provides the foundation for cryptographic applications and mathematical computation.
3

4
## Capabilities
5

6
### Prime Number Generation
7

8
Functions for generating prime numbers, including special classes of primes and efficient algorithms for finding primes within specified ranges.
9

10
```python { .api }
11
def primes(n):
12
    """
13
    Generate all prime numbers ≤ n using sieve methods.
14
    
15
    Parameters:
16
    - n (int): Upper bound for prime generation
17
    
18
    Returns:
19
    np.ndarray: Array of all primes ≤ n
20
    """
21

22
def kth_prime(k):
23
    """
24
    Return the k-th prime number (1-indexed).
25
    
26
    Parameters:
27
    - k (int): Prime index (k=1 returns 2, k=2 returns 3, etc.)
28
    
29
    Returns:
30
    int: The k-th prime number
31
    """
32

33
def prev_prime(n):
34
    """
35
    Find the largest prime ≤ n.
36
    
37
    Parameters:
38
    - n (int): Upper bound
39
    
40
    Returns:
41
    int: Largest prime ≤ n, or None if no such prime exists
42
    """
43

44
def next_prime(n):
45
    """
46
    Find the smallest prime > n.
47
    
48
    Parameters:
49
    - n (int): Lower bound
50
    
51
    Returns:
52
    int: Smallest prime > n
53
    """
54

55
def random_prime(bits, seed=None):
56
    """
57
    Generate random prime with specified bit length.
58
    
59
    Parameters:
60
    - bits (int): Bit length of desired prime
61
    - seed (int, optional): Random seed for reproducibility
62
    
63
    Returns:
64
    int: Random prime with specified bit length
65
    """
66

67
def mersenne_exponents(n=None):
68
    """
69
    Return known Mersenne prime exponents.
70
    
71
    Parameters:
72
    - n (int, optional): Maximum exponent, returns all known if None
73
    
74
    Returns:
75
    list: Known Mersenne prime exponents p where 2^p - 1 is prime
76
    """
77

78
def mersenne_primes(n=None):
79
    """
80
    Return known Mersenne primes.
81
    
82
    Parameters:
83
    - n (int, optional): Maximum exponent, returns all known if None
84
    
85
    Returns:
86
    list: Known Mersenne primes 2^p - 1
87
    """
88
```
89

90
### Primality Testing
91

92
Probabilistic and deterministic algorithms for testing whether numbers are prime, composite, or have special properties.
93

94
```python { .api }
95
def is_prime(n, rounds=None):
96
    """
97
    Test if integer is prime using multiple algorithms.
98
    
99
    Combines deterministic tests for small numbers with probabilistic
100
    Miller-Rabin test for larger numbers.
101
    
102
    Parameters:
103
    - n (int): Integer to test
104
    - rounds (int, optional): Number of Miller-Rabin rounds for probabilistic test
105
    
106
    Returns:
107
    bool: True if n is prime, False otherwise
108
    """
109

110
def is_composite(n):
111
    """
112
    Test if integer is composite (not prime and > 1).
113
    
114
    Parameters:
115
    - n (int): Integer to test
116
    
117
    Returns:
118
    bool: True if n is composite, False otherwise
119
    """
120

121
def fermat_primality_test(n, a, rounds=1):
122
    """
123
    Fermat primality test with base a.
124
    
125
    Tests if a^(n-1) ≡ 1 (mod n). If False, n is definitely composite.
126
    If True, n is probably prime.
127
    
128
    Parameters:
129
    - n (int): Integer to test
130
    - a (int): Test base, must be coprime to n
131
    - rounds (int): Number of test rounds, default 1
132
    
133
    Returns:
134
    bool: True if n passes Fermat test, False if definitely composite
135
    """
136

137
def miller_rabin_primality_test(n, a=None, rounds=None):
138
    """
139
    Miller-Rabin probabilistic primality test.
140
    
141
    More reliable than Fermat test, can detect Carmichael numbers.
142
    
143
    Parameters:
144
    - n (int): Integer to test  
145
    - a (int, optional): Test base, random if None
146
    - rounds (int, optional): Number of test rounds, auto-selected if None
147
    
148
    Returns:
149
    bool: True if n is probably prime, False if definitely composite
150
    """
151

152
def is_prime_power(n):
153
    """
154
    Test if n = p^k for some prime p and integer k ≥ 1.
155
    
156
    Parameters:
157
    - n (int): Integer to test
158
    
159
    Returns:
160
    bool or tuple: False if not prime power, (p, k) if n = p^k
161
    """
162

163
def is_perfect_power(n):
164
    """
165
    Test if n = c^e for integers c, e with e > 1.
166
    
167
    Parameters:
168
    - n (int): Integer to test
169
    
170
    Returns:
171
    bool or tuple: False if not perfect power, (c, e) if n = c^e
172
    """
173

174
def is_smooth(n, B):
175
    """
176
    Test if n is B-smooth (all prime factors ≤ B).
177
    
178
    Parameters:
179
    - n (int): Integer to test
180
    - B (int): Smoothness bound
181
    
182
    Returns:
183
    bool: True if all prime factors of n are ≤ B
184
    """
185

186
def is_powersmooth(n, B):
187
    """
188
    Test if n is B-powersmooth (all prime powers in factorization ≤ B).
189
    
190
    Parameters:
191
    - n (int): Integer to test  
192
    - B (int): Powersmooth bound
193
    
194
    Returns:
195
    bool: True if all prime powers dividing n are ≤ B
196
    """
197
```
198

199
### Modular Arithmetic
200

201
Functions for modular arithmetic operations, including totient functions, primitive roots, and multiplicative group analysis.
202

203
```python { .api }
204
def totatives(n):
205
    """
206
    Find all totatives of n (integers 1 ≤ k < n where gcd(k, n) = 1).
207
    
208
    Parameters:
209
    - n (int): Modulus
210
    
211
    Returns:
212
    np.ndarray: Array of totatives of n
213
    """
214

215
def euler_phi(n):
216
    """
217
    Euler's totient function φ(n) - count of totatives of n.
218
    
219
    Parameters:
220
    - n (int): Input integer
221
    
222
    Returns:
223
    int: Number of integers 1 ≤ k ≤ n where gcd(k, n) = 1
224
    """
225

226
def carmichael_lambda(n):
227
    """
228
    Carmichael lambda function λ(n) - exponent of multiplicative group Z*_n.
229
    
230
    Parameters:
231
    - n (int): Input integer
232
    
233
    Returns:
234
    int: λ(n), the exponent of the multiplicative group modulo n
235
    """
236

237
def is_cyclic(n):
238
    """
239
    Test if multiplicative group Z*_n is cyclic.
240
    
241
    Parameters:
242
    - n (int): Modulus
243
    
244
    Returns:
245
    bool: True if Z*_n is cyclic (has primitive roots)
246
    """
247

248
def primitive_root(n, start=1, stop=None, method="min"):
249
    """
250
    Find primitive root modulo n (generator of Z*_n).
251
    
252
    Parameters:
253
    - n (int): Modulus (must have primitive roots)
254
    - start (int): Start search from this value, default 1
255
    - stop (int, optional): Stop search at this value, defaults to n
256
    - method (str): Search method - "min", "max", or "random", default "min"
257
    
258
    Returns:
259
    int: A primitive root modulo n
260
    
261
    Raises:
262
    ValueError: If n has no primitive roots
263
    """
264

265
def primitive_roots(n, start=1, stop=None, reverse=False):
266
    """
267
    Iterator through all primitive roots modulo n.
268
    
269
    Parameters:
270
    - n (int): Modulus (must have primitive roots)
271
    - start (int): Start iteration from this value, default 1
272
    - stop (int, optional): Stop iteration at this value, defaults to n
273
    - reverse (bool): Iterate in reverse order, default False
274
    
275
    Yields:
276
    int: Primitive roots modulo n
277
    """
278

279
def is_primitive_root(g, n):
280
    """
281
    Test if g is a primitive root modulo n.
282
    
283
    Parameters:
284
    - g (int): Candidate primitive root
285
    - n (int): Modulus
286
    
287
    Returns:
288
    bool: True if g is a primitive root modulo n
289
    """
290
```
291

292
### Number-Theoretic Symbols
293

294
Functions for computing Legendre, Jacobi, and Kronecker symbols used in quadratic residue theory and cryptography.
295

296
```python { .api }
297
def legendre_symbol(a, p):
298
    """
299
    Compute Legendre symbol (a/p).
300
    
301
    The Legendre symbol indicates whether a is a quadratic residue modulo p.
302
    
303
    Parameters:
304
    - a (int): Numerator
305
    - p (int): Denominator (must be odd prime)
306
    
307
    Returns:
308
    int: -1, 0, or 1
309
      - 1 if a is quadratic residue mod p
310
      - -1 if a is quadratic non-residue mod p  
311
      - 0 if a ≡ 0 (mod p)
312
    """
313

314
def jacobi_symbol(a, n):
315
    """
316
    Compute Jacobi symbol (a/n).
317
    
318
    Generalization of Legendre symbol to composite moduli.
319
    
320
    Parameters:
321
    - a (int): Numerator
322
    - n (int): Denominator (must be positive odd)
323
    
324
    Returns:
325
    int: -1, 0, or 1
326
    """
327

328
def kronecker_symbol(a, n):
329
    """
330
    Compute Kronecker symbol (a/n).
331
    
332
    Further generalization of Jacobi symbol to all integer moduli.
333
    
334
    Parameters:
335
    - a (int): Numerator
336
    - n (int): Denominator
337
    
338
    Returns:
339
    int: -1, 0, or 1
340
    """
341
```
342

343
### Integer Factorization
344

345
Algorithms for factoring integers and analyzing their divisor structure.
346

347
```python { .api }
348
def perfect_power(n):
349
    """
350
    Find perfect power decomposition n = c^e with maximal e > 1.
351
    
352
    Parameters:
353
    - n (int): Integer to decompose
354
    
355
    Returns:
356
    tuple or None: (c, e) if n = c^e with e > 1, None otherwise
357
    """
358

359
def trial_division(n, B=None):
360
    """
361
    Factor n using trial division up to bound B.
362
    
363
    Parameters:
364
    - n (int): Integer to factor
365
    - B (int, optional): Trial division bound, defaults to sqrt(n)
366
    
367
    Returns:
368
    tuple: (factors, remainder) where factors are found prime factors
369
    """
370

371
def pollard_p1(n, B=None, B2=None):
372
    """
373
    Pollard p-1 factorization algorithm.
374
    
375
    Effective when n has prime factor p where p-1 is smooth.
376
    
377
    Parameters:
378
    - n (int): Integer to factor
379
    - B (int, optional): Stage 1 bound
380
    - B2 (int, optional): Stage 2 bound
381
    
382
    Returns:
383
    int or None: Non-trivial factor of n, or None if unsuccessful
384
    """
385

386
def pollard_rho(n, c=1):
387
    """
388
    Pollard rho factorization algorithm.
389
    
390
    Probabilistic algorithm effective for finding small to medium factors.
391
    
392
    Parameters:
393
    - n (int): Integer to factor
394
    - c (int): Parameter for iteration function, default 1
395
    
396
    Returns:
397
    int or None: Non-trivial factor of n, or None if unsuccessful
398
    """
399

400
def divisors(n):
401
    """
402
    Find all positive divisors of n.
403
    
404
    Parameters:
405
    - n (int): Integer to find divisors of
406
    
407
    Returns:
408
    list: All positive divisors of n in ascending order
409
    """
410

411
def divisor_sigma(n, k=1):
412
    """
413
    Compute sum of k-th powers of divisors: σ_k(n) = Σ d^k.
414
    
415
    Parameters:
416
    - n (int): Integer
417
    - k (int): Power, default 1
418
    
419
    Returns:
420
    int: σ_k(n) = sum of d^k over all divisors d of n
421
    """
422
```
423

424
### Polymorphic Mathematical Functions
425

426
Generic mathematical functions that work with both integers and polynomials over finite fields.
427

428
```python { .api }
429
def gcd(*args):
430
    """
431
    Greatest common divisor of integers or polynomials.
432
    
433
    Parameters:
434
    - *args: Integers or polynomials (must be same type)
435
    
436
    Returns:
437
    int or Poly: GCD of all arguments
438
    """
439

440
def egcd(a, b):
441
    """
442
    Extended Euclidean algorithm: find gcd(a,b) and Bézout coefficients.
443
    
444
    Parameters:
445
    - a, b: Integers or polynomials (same type)
446
    
447
    Returns:
448
    tuple: (gcd, x, y) where gcd = x*a + y*b
449
    """
450

451
def lcm(*args):
452
    """
453
    Least common multiple of integers or polynomials.
454
    
455
    Parameters:
456
    - *args: Integers or polynomials (must be same type)
457
    
458
    Returns:
459
    int or Poly: LCM of all arguments
460
    """
461

462
def prod(args):
463
    """
464
    Product of sequence of integers or polynomials.
465
    
466
    Parameters:
467
    - args: Iterable of integers or polynomials
468
    
469
    Returns:
470
    int or Poly: Product of all elements
471
    """
472

473
def are_coprime(*args):
474
    """
475
    Test if arguments are pairwise coprime.
476
    
477
    Parameters:
478
    - *args: Integers or polynomials (must be same type)
479
    
480
    Returns:
481
    bool: True if gcd of any pair equals 1
482
    """
483

484
def crt(remainders, moduli):
485
    """
486
    Chinese Remainder Theorem solver.
487
    
488
    Find x such that x ≡ remainders[i] (mod moduli[i]) for all i.
489
    
490
    Parameters:
491
    - remainders: Sequence of remainders
492
    - moduli: Sequence of pairwise coprime moduli
493
    
494
    Returns:
495
    int: Solution x in range [0, product(moduli))
496
    """
497

498
def factors(n):
499
    """
500
    Prime factorization of integer or irreducible factorization of polynomial.
501
    
502
    Parameters:
503
    - n: Integer or polynomial
504
    
505
    Returns:
506
    dict: Factorization as {factor: multiplicity, ...}
507
    """
508

509
def is_square_free(n):
510
    """
511
    Test if integer/polynomial is square-free (no repeated prime factors).
512
    
513
    Parameters:
514
    - n: Integer or polynomial
515
    
516
    Returns:
517
    bool: True if square-free
518
    """
519
```
520

521
### Integer Mathematics
522

523
Basic integer mathematical functions for roots, logarithms, and special computations.
524

525
```python { .api }
526
def isqrt(n):
527
    """
528
    Integer square root: largest integer k such that k² ≤ n.
529
    
530
    Parameters:
531
    - n (int): Non-negative integer
532
    
533
    Returns:
534
    int: floor(sqrt(n))
535
    """
536

537
def iroot(n, k):
538
    """
539
    Integer k-th root: largest integer r such that r^k ≤ n.
540
    
541
    Parameters:
542
    - n (int): Non-negative integer
543
    - k (int): Root degree (positive integer)
544
    
545
    Returns:
546
    int: floor(n^(1/k))
547
    """
548

549
def ilog(n, base):
550
    """
551
    Integer logarithm: largest integer k such that base^k ≤ n.
552
    
553
    Parameters:
554
    - n (int): Positive integer
555
    - base (int): Logarithm base (> 1)
556
    
557
    Returns:
558
    int: floor(log_base(n))
559
    """
560
```
561

562
## Usage Examples
563

564
### Prime Generation and Testing
565

566
```python
567
import galois
568

569
# Generate all primes up to 100
570
prime_list = galois.primes(100)
571
print(f"Primes ≤ 100: {prime_list}")
572

573
# Get specific primes
574
print(f"10th prime: {galois.kth_prime(10)}")
575
print(f"Previous prime ≤ 50: {galois.prev_prime(50)}")
576
print(f"Next prime > 50: {galois.next_prime(50)}")
577

578
# Generate random prime
579
random_prime = galois.random_prime(bits=256)
580
print(f"256-bit random prime: {random_prime}")
581

582
# Test primality
583
candidates = [97, 98, 99, 101]
584
for n in candidates:
585
    print(f"{n} is prime: {galois.is_prime(n)}")
586

587
# Mersenne primes
588
mersenne_exp = galois.mersenne_exponents(100)
589
print(f"Mersenne exponents ≤ 100: {mersenne_exp}")
590
```
591

592
### Modular Arithmetic
593

594
```python
595
import galois
596

597
# Euler's totient function
598
n = 12
599
phi_n = galois.euler_phi(n)
600
totatives_n = galois.totatives(n)
601
print(f"φ({n}) = {phi_n}")
602
print(f"Totatives of {n}: {totatives_n}")
603

604
# Primitive roots
605
p = 17  # Prime with primitive roots
606
print(f"Is Z*_{p} cyclic: {galois.is_cyclic(p)}")
607

608
prim_root = galois.primitive_root(p)
609
print(f"Primitive root mod {p}: {prim_root}")
610

611
all_prim_roots = list(galois.primitive_roots(p))
612
print(f"All primitive roots mod {p}: {all_prim_roots}")
613

614
# Test specific element
615
g = 3
616
is_prim = galois.is_primitive_root(g, p)
617
print(f"{g} is primitive root mod {p}: {is_prim}")
618
```
619

620
### Factorization
621

622
```python
623
import galois
624

625
# Factor integers
626
n = 1024
627
factors_dict = galois.factors(n)
628
print(f"Prime factorization of {n}: {factors_dict}")
629

630
# Trial division
631
n = 12345
632
factors, remainder = galois.trial_division(n, B=100)
633
print(f"Trial division of {n}: factors={factors}, remainder={remainder}")
634

635
# Pollard methods for larger numbers
636
n = 8051  # 83 * 97
637
factor = galois.pollard_rho(n)
638
print(f"Pollard rho found factor of {n}: {factor}")
639

640
# Divisors
641
n = 60
642
all_divisors = galois.divisors(n)
643
sigma_1 = galois.divisor_sigma(n, 1)  # Sum of divisors
644
print(f"Divisors of {n}: {all_divisors}")
645
print(f"Sum of divisors σ₁({n}): {sigma_1}")
646
```
647

648
### Number-Theoretic Symbols
649

650
```python
651
import galois
652

653
# Legendre symbol
654
a, p = 7, 11
655
legendre = galois.legendre_symbol(a, p)
656
print(f"Legendre symbol ({a}/{p}): {legendre}")
657

658
# Jacobi symbol  
659
a, n = 7, 15
660
jacobi = galois.jacobi_symbol(a, n)
661
print(f"Jacobi symbol ({a}/{n}): {jacobi}")
662

663
# Test quadratic residues
664
p = 13
665
quadratic_residues = []
666
for a in range(1, p):
667
    if galois.legendre_symbol(a, p) == 1:
668
        quadratic_residues.append(a)
669
print(f"Quadratic residues mod {p}: {quadratic_residues}")
670
```
671

672
### Chinese Remainder Theorem
673

674
```python
675
import galois
676

677
# Solve system of congruences
678
# x ≡ 2 (mod 3)
679
# x ≡ 3 (mod 5)  
680
# x ≡ 2 (mod 7)
681
remainders = [2, 3, 2]
682
moduli = [3, 5, 7]
683

684
solution = galois.crt(remainders, moduli)
685
print(f"CRT solution: x ≡ {solution} (mod {galois.prod(moduli)})")
686

687
# Verify solution
688
for r, m in zip(remainders, moduli):
689
    print(f"{solution} ≡ {solution % m} (mod {m}) [expected {r}]")
690
```
691

692
### GCD and Extended Euclidean Algorithm
693

694
```python
695
import galois
696

697
# Basic GCD
698
a, b, c = 48, 18, 24
699
result = galois.gcd(a, b, c)
700
print(f"gcd({a}, {b}, {c}) = {result}")
701

702
# Extended GCD with Bézout coefficients
703
a, b = 48, 18
704
gcd_val, x, y = galois.egcd(a, b)
705
print(f"gcd({a}, {b}) = {gcd_val}")
706
print(f"Bézout coefficients: {gcd_val} = {x}×{a} + {y}×{b}")
707
print(f"Verification: {x*a + y*b} = {gcd_val}")
708

709
# LCM
710
lcm_val = galois.lcm(a, b, c)
711
print(f"lcm({a}, {b}, {c}) = {lcm_val}")
712

713
# Test coprimality
714
numbers = [9, 16, 25]
715
coprime = galois.are_coprime(*numbers)
716
print(f"Are {numbers} pairwise coprime: {coprime}")
717
```
718

719
### Integer Roots and Logarithms
720

721
```python
722
import galois
723

724
# Integer square root
725
n = 15
726
sqrt_n = galois.isqrt(n)
727
print(f"isqrt({n}) = {sqrt_n} (since {sqrt_n}² = {sqrt_n**2} ≤ {n} < {(sqrt_n+1)**2})")
728

729
# Integer k-th roots
730
n, k = 100, 3
731
root_n = galois.iroot(n, k)
732
print(f"iroot({n}, {k}) = {root_n} (since {root_n}³ = {root_n**k} ≤ {n})")
733

734
# Integer logarithms
735
n, base = 1000, 10
736
log_n = galois.ilog(n, base)
737
print(f"ilog({n}, {base}) = {log_n} (since {base}^{log_n} = {base**log_n} ≤ {n})")
738
```
739

740
## Performance and Implementation Notes
741

742
- **Sieve Methods**: Efficient prime generation using optimized sieve algorithms
743
- **Probabilistic Tests**: Miller-Rabin with appropriate round counts for security levels
744
- **Optimized Algorithms**: Fast modular exponentiation and extended GCD implementations
745
- **Memory Efficiency**: Lazy evaluation for large prime iterations
746
- **Cryptographic Quality**: Suitable randomness sources for cryptographic prime generation
747
- **Field Integration**: Seamless operation with finite field elements where applicable