# Number Theory

## with SAGE

### Contents

[Prev] [Intro] [Next]

SAGE has basic commands and subroutines which implement all basic number theory (and many other things). A general introduction to SAGE use and how SAGE can be accessed at UMBC can be found in a separate document.

SAGE supports the use of integers of arbitrary length and a number of other basic data types such as polynomials. In addition to SAGE's native functions, the extensive functionality of PARI/GP can also be used from within SAGE (tutorial).

## Integer Operations

The basic arithmetic operations, `+, -, *, /`, are available from the keyboard. Note that the division operator returns a quotient (eg `7/2` returns 7/2) so the answer should be truncated with the "truncating division" `//` function (eg `7//2` returns 3). The modular reduction operator is represented by the `mod` command (eg `mod(7,2)` returns 1). Exponentiation is represented with the `^` symbol. The pow() function can be used for exponentiation if called with two parameters, or efficient modular exponentiation if called with three parameters. (pow(a,b,c) returns the same result as ((a^b) % c), but is much more efficient.)

Examples:

• `mod(8 + 7, 13)`
• `mod(5 * 4, 13)`
• `7 // 2`
• `3^2`
• `pow(2,6,11)`

In addition to the basic arithmetic operations for integers SAGE implements gcd and the extended gcd operations as basic functions:

Examples:
• `gcd(15,6)`
• `xgcd(23,7)`
(returns an answer of (1, 4, -13), as 23*4 + 7*(-13)=1)

Modular inversion is implemented with the `inverse_mod` command, so
`inverse_mod(3,19)`
returns the value 13.

Factoring and primality testing are implemented with the factor() and is_prime() commands. Examples of their use are
`factor(123456)`
which returns 2^6 * 3 * 643 and the command
`is_prime(643)`
which returns the value `True`.

Other useful functions include:

• `prime_factors(1234)`
• `divisors(1234)`
• `next_prime(1234)`
• `nth_prime(10)`
• `is_prime(1234)`
• `prime_pi(1234)` - number of primes less than n
• `euler_phi(1234)` - φ(n) = number of integers coprime to n
• `CRT(ra,rb,a,b)` - Chinese Remainder Theorem: find ra, rb such that ra*a+rb*b = gcd(a,b)
• `primitive_root(p)` - generator of Zp* for prime p
• `continued_fraction(x)`
• `5.jacobi(1237)` - Jacobi/Legendre/Kronecker symbol (also `kronecker(5,1237)`
• `mod(123,1001).sqrt()` - Modular square root (Shanks/Tonelli/RESSOL algorithm)

## Polynomial Operations

Polynomials are constructed explicitly in SAGE, with the coefficients specified. This is similar to MAGMA, but differs from the much more casual approach in Maple and Mathematica, where polynomials are viewed formally as symbols with at most a generic treatment of coefficients. Some examples are:

• Z[x] (Integer coefficients) ``` sage: ZP.<x> = ZZ[] sage: (x^5 + 3*x^2 - 2*x + 7) // (x + 1) sage: (x^5 + 3*x^2 - 2*x + 7).quo_rem(x + 1) sage: gcd(3*x^2+6*x-9,5*x^3-2*x+2) sage: xgcd(3*x^2+6*x-9,5*x^3-2*x+2) sage: factor(3*x^5+5*x-8) sage: (3*x^5+5*x-8).factor_mod(3)``` (Factor with coeffs in Z3)
• Z5[y] (Coefficients mod 5) ``` sage: Z5P.<y> = GF(5)[] sage: Z5P([1,2,3]) sage: (y^5 + 3*y^2 - 2*y + 7) // (y + 1) sage: xgcd(3*y^2+6*y-9,2*y^3-2*y+3) sage: factor(2*y^3-2*y+3) sage: (2*y^3-2*y+3).is_irreducible()``` (is it irreducible?)

Polynomials are manipulated and combined with the same basic operations as are used on integers, `+, -, *, ^, //, %`.

Polynomials with coefficients over a finite field can be enumerated with the following trick: sage: Z2P.<z> = GF(2)[] sage: for num in range(32): # Enumerate polys of degree 5 ....: print z^5 + sum([coeff*z^i for i,coeff in enumerate([(num//2**i)%2 for i in range(5)])])

and a search for low density primitive polynomials can be implemented with a python generator for low-density lists (ref). Here we list all density 5 primitive polynomials over GF2 of degree 32. You can also check and confirm that there are no lower density primitive polynomials of that degree. sage: def fixeddensity(len, density): ....: if len == density: yield *len ....: elif density == 0: yield *len ....: else: ....: for leftlist in fixeddensity(len-1, density): ....: yield leftlist +  ....: for leftlist in fixeddensity(len-1, density-1): ....: yield leftlist +  sage: for it in fixeddensity(31,4): ....: thepoly = z^32 + sum([coeff*z^i for i,coeff in enumerate(it)]) ....: if thepoly.is_primitive(): print thepoly

## Gaussian Integer Operations

SAGE currently has no convenient and complete way to work with Gaussian integers.

A kludge which allows simple operations with the Gaussian integers is to represent them as a formal polynomial ring quotient. This allows you to add and multiply elements, but there is no concept of division, modular reduction or factoring:

``````sage: ZZP = PolynomialRing(ZZ,'x'); x = ZZP.gen(); GI = ZZP.quotient(x^2 + 1, 'i'); i=GI.gen()
sage: (i+2)*(3*i-5)
sage: (i+2)^10
``````

A more exotic representation (but more appropriate for more general quadratic number fields) is to represent the Gaussian Integers as the maximal order in the number field Q(√-1).

``````sage: QFi.<sr1> = NumberField(x^2 + 1)
sage: GI = QFi.order(sr1)     # Create the Gaussian Integers
sage: i1 = QFi(sr1)
sage: (2+i1)^10
sage: 120 // (2-sr1)
sage: factor(GI.ideal(25+sr1))  # Factor as ideals
# So 25 + sqrt(-1) = (1 + 4 sqrt(-1))(1 + sqrt(-1))(1 + 2 sqrt(-1))
(Fractional ideal (4*sr1 + 1)) * (Fractional ideal (sr1 + 1)) * (Fractional ideal (2*sr1 + 1))
``````

The ideal factorization is the best we can do for all of the fields without unique factorization, and it is not too hard to translate back to usual factorization by eye in the case of the Gaussian integers.

Some examples are the Eisenstein Integers, the ring of integers in Q(√-3), and the integers in Q(√-5), where unique factorization fails, so many of the ideals are not principal (cannot be generated by a single number in the ring).

``````sage: QF3.<sr3> = NumberField(x^2 + 3)
sage: ZQ3 = QF3.maximal_order()
sage: ZQ3.basis()
[1/2*sr3 + 1/2, sr3]
sage: factor(ZQ3.ideal(7))   # 7 splits as 7 = (2 - sqrt(-3))(2 + sqrt(-3))
(Fractional ideal (sr3 - 2)) * (Fractional ideal (-sr3 - 2))
``````
```sage: QF5.<sr5> = NumberField(x^2 + 5) sage: ZQ5 = QF5.maximal_order() sage: ZQ5.basis() [1, sr5] sage: factor(ZQ5.ideal(12)) (Fractional ideal (2, sr5 + 1))^4 * (Fractional ideal (3, sr5 + 1)) * (Fractional ideal (3, sr5 - 1)) ```

## Finite Field Operations

Finite fields are cleanly represented in SAGE by the `FiniteField` (equivalently `GF`) class. A finite field can be created with SAGE choosing the representation, or a specific defining polynomial can be specified. This is well documented in the SAGE reference manual and in API documentation.

``````sage: GF125a.<b> = GF(5^3)  # Let SAGE choose the representation
sage: b^15
sage: GF5.<x5>=GF(5)[]; GF125.<a> = GF(5^3,'a',modulus=x5^3+x5+1)
sage: GF125.multiplicative_generator()
sage: for i in range(124):  # Compute order(a) the hard way
print i,a^i
sage: a.multiplicative_order()  # The easy way``````

You can also specify the defining polynomial for the finite field explicitly. The following code illustrates how to search for a primitive polynomial (i.e. so that the root w has full multiplicative order in the resulting Galois field) of a particular form.

``````sage: Z2P.<z> = GF(2)[]
sage: for num in range(2**10):  # Find primitive degree 128 polys
thepoly = z^128 + sum([coeff*z^i for i,coeff in enumerate([(num//2**i)%2 for i in range(10)])])
if thepoly.is_primitive(): print thepoly
sage: GF2p128.<w> = GF(2^128,'w',modulus=z^128+z^7+z^2+z+1)``````

And now you can confirm that the root w has full order, by checking all possible factors of the group order.

``````sage: factor(2**128-1)
3 * 5 * 17 * 257 * 641 * 65537 * 274177 * 6700417 * 67280421310721
sage: w^(2^128-1)
1
sage: w^((2^128-1)/3)
w^125 + w^123 + w^120 + w^118 + ... (certainly not 1)
sage: w^((2^128-1)/5)
w^126 + w^125 + w^124 + w^123 + ... (again not 1)
sage: for p in prime_factors(2**128-1): # Check all possible cofactors
....:     if w^((2^128-1)/p)==1:
....:         print "Not primitive, w^((2**128-1)/",p," is 1"``````

(Confirm the operation of the code by running it against (w+1) instead of w, observing that (w+1) does not have full order.)

## Continued Fraction Operations

The `continued_fraction` function produces continued fraction expansions and the `value` function can convert continued fraction expansions to numbers. SAGE currently does not correctly deal with periodic continued fractions produced by quadratic surds. ``` sage: a = continued_fraction(74/42); a [1, 1, 3, 5] sage: a.value() 37/21 sage: a=continued_fraction(sqrt(6)); a # a truncated value [2, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 4, 2, 2, 1] sage: a.value() 159018721/64919121 sage: a=continued_fraction(pi,bits=30); a # 30-bit accuracy (about 9 digits) [3, 7, 15, 1, 294] sage: float(pi-a.value()) -1.2349787859022854e-09 sage: convergents(a) [3, 22/7, 333/106, 355/113, 104703/33328] ```

[Prev] [Intro] [Next]
Robert Campbell
10 March, 2009