Returns a solution to a Chinese Remainder Theorem problem.
INPUT:
OUTPUT:
If m, n are not None, returns a solution \(x\) to the simultaneous congruences \(x\equiv a \bmod m\) and \(x\equiv b \bmod n\), if one exists. By the Chinese Remainder Theorem, a solution to the simultaneous congruences exists if and only if \(a\equiv b\pmod{\gcd(m,n)}\). The solution \(x\) is only well-defined modulo \(\text{lcm}(m,n)\).
If a and b are lists, returns a simultaneous solution to the congruences \(x\equiv a_i\pmod{b_i}\), if one exists.
See also
EXAMPLES:
Using crt by giving it pairs of residues and moduli:
sage: crt(2, 1, 3, 5)
11
sage: crt(13, 20, 100, 301)
28013
sage: crt([2, 1], [3, 5])
11
sage: crt([13, 20], [100, 301])
28013
You can also use upper case:
sage: c = CRT(2,3, 3, 5); c
8
sage: c % 3 == 2
True
sage: c % 5 == 3
True
Note that this also works for polynomial rings:
sage: K.<a> = NumberField(x^3 - 7)
sage: R.<y> = K[]
sage: f = y^2 + 3
sage: g = y^3 - 5
sage: CRT(1,3,f,g)
-3/26*y^4 + 5/26*y^3 + 15/26*y + 53/26
sage: CRT(1,a,f,g)
(-3/52*a + 3/52)*y^4 + (5/52*a - 5/52)*y^3 + (15/52*a - 15/52)*y + 27/52*a + 25/52
You can also do this for any number of moduli:
sage: K.<a> = NumberField(x^3 - 7)
sage: R.<x> = K[]
sage: CRT([], [])
0
sage: CRT([a], [x])
a
sage: f = x^2 + 3
sage: g = x^3 - 5
sage: h = x^5 + x^2 - 9
sage: k = CRT([1, a, 3], [f, g, h]); k
(127/26988*a - 5807/386828)*x^9 + (45/8996*a - 33677/1160484)*x^8 + (2/173*a - 6/173)*x^7 + (133/6747*a - 5373/96707)*x^6 + (-6/2249*a + 18584/290121)*x^5 + (-277/8996*a + 38847/386828)*x^4 + (-135/4498*a + 42673/193414)*x^3 + (-1005/8996*a + 470245/1160484)*x^2 + (-1215/8996*a + 141165/386828)*x + 621/8996*a + 836445/386828
sage: k.mod(f)
1
sage: k.mod(g)
a
sage: k.mod(h)
3
If the moduli are not coprime, a solution may not exist:
sage: crt(4,8,8,12)
20
sage: crt(4,6,8,12)
Traceback (most recent call last):
...
ValueError: No solution to crt problem since gcd(8,12) does not divide 4-6
sage: x = polygen(QQ)
sage: crt(2,3,x-1,x+1)
-1/2*x + 5/2
sage: crt(2,x,x^2-1,x^2+1)
-1/2*x^3 + x^2 + 1/2*x + 1
sage: crt(2,x,x^2-1,x^3-1)
Traceback (most recent call last):
...
ValueError: No solution to crt problem since gcd(x^2 - 1,x^3 - 1) does not divide 2-x
sage: crt(int(2), int(3), int(7), int(11))
58
Returns a CRT basis for the given moduli.
INPUT:
extended Euclidean algorithm
OUTPUT:
Note
The pairwise coprimality of the input is not checked.
EXAMPLES:
sage: a1 = ZZ(mod(42,5))
sage: a2 = ZZ(mod(42,13))
sage: c1,c2 = CRT_basis([5,13])
sage: mod(a1*c1+a2*c2,5*13)
42
A polynomial example:
sage: x=polygen(QQ)
sage: mods = [x,x^2+1,2*x-3]
sage: b = CRT_basis(mods)
sage: b
[-2/3*x^3 + x^2 - 2/3*x + 1, 6/13*x^3 - x^2 + 6/13*x, 8/39*x^3 + 8/39*x]
sage: [[bi % mj for mj in mods] for bi in b]
[[1, 0, 0], [0, 1, 0], [0, 0, 1]]
Given a list v of elements and a list of corresponding moduli, find a single element that reduces to each element of v modulo the corresponding moduli.
See also
EXAMPLES:
sage: CRT_list([2,3,2], [3,5,7])
23
sage: x = polygen(QQ)
sage: c = CRT_list([3], [x]); c
3
sage: c.parent()
Univariate Polynomial Ring in x over Rational Field
It also works if the moduli are not coprime:
sage: CRT_list([32,2,2],[60,90,150])
452
But with non coprime moduli there is not always a solution:
sage: CRT_list([32,2,1],[60,90,150])
Traceback (most recent call last):
...
ValueError: No solution to crt problem since gcd(180,150) does not divide 92-1
The arguments must be lists:
sage: CRT_list([1,2,3],"not a list")
Traceback (most recent call last):
...
ValueError: Arguments to CRT_list should be lists
sage: CRT_list("not a list",[2,3])
Traceback (most recent call last):
...
ValueError: Arguments to CRT_list should be lists
The list of moduli must have the same length as the list of elements:
sage: CRT_list([1,2,3],[2,3,5])
23
sage: CRT_list([1,2,3],[2,3])
Traceback (most recent call last):
...
ValueError: Arguments to CRT_list should be lists of the same length
sage: CRT_list([1,2,3],[2,3,5,7])
Traceback (most recent call last):
...
ValueError: Arguments to CRT_list should be lists of the same length
TESTS:
sage: CRT([32r,2r,2r],[60r,90r,150r])
452
Vector form of the Chinese Remainder Theorem: given a list of integer vectors \(v_i\) and a list of coprime moduli \(m_i\), find a vector \(w\) such that \(w = v_i \pmod m_i\) for all \(i\). This is more efficient than applying CRT() to each entry.
INPUT:
OUTPUT:
EXAMPLES:
sage: CRT_vectors([[3,5,7],[3,5,11]], [2,3])
[3, 5, 5]
sage: CRT_vectors([vector(ZZ, [2,3,1]), Sequence([1,7,8],ZZ)], [8,9])
[10, 43, 17]
Return the value of the Euler phi function on the integer n. We defined this to be the number of positive integers <= n that are relatively prime to n. Thus if n<=0 then euler_phi(n) is defined and equals 0.
INPUT:
EXAMPLES:
sage: euler_phi(1)
1
sage: euler_phi(2)
1
sage: euler_phi(3)
2
sage: euler_phi(12)
4
sage: euler_phi(37)
36
Notice that euler_phi is defined to be 0 on negative numbers and 0.
sage: euler_phi(-1)
0
sage: euler_phi(0)
0
sage: type(euler_phi(0))
<type 'sage.rings.integer.Integer'>
We verify directly that the phi function is correct for 21.
sage: euler_phi(21)
12
sage: [i for i in range(21) if gcd(21,i) == 1]
[1, 2, 4, 5, 8, 10, 11, 13, 16, 17, 19, 20]
The length of the list of integers ‘i’ in range(n) such that the gcd(i,n) == 1 equals euler_phi(n).
sage: len([i for i in range(21) if gcd(21,i) == 1]) == euler_phi(21)
True
The phi function also has a special plotting method.
sage: P = plot(euler_phi, -3, 71)
AUTHORS:
Plot the Euler phi function.
INPUT:
EXAMPLES:
sage: p = Euler_Phi().plot()
sage: p.ymax()
46.0
The greatest common divisor of a and b, or if a is a list and b is omitted the greatest common divisor of all elements of a.
INPUT:
Additional keyword arguments are passed to the respectively called methods.
OUTPUT:
The given elements are first coerced into a common parent. Then, their greatest common divisor in that common parent is returned.
EXAMPLES:
sage: GCD(97,100)
1
sage: GCD(97*10^15, 19^20*97^2)
97
sage: GCD(2/3, 4/5)
2/15
sage: GCD([2,4,6,8])
2
sage: GCD(srange(0,10000,10)) # fast !!
10
Note that to take the gcd of \(n\) elements for \(n \not= 2\) you must put the elements into a list by enclosing them in [..]. Before #4988 the following wrongly returned 3 since the third parameter was just ignored:
sage: gcd(3,6,2)
Traceback (most recent call last):
...
TypeError: gcd() takes at most 2 arguments (3 given)
sage: gcd([3,6,2])
1
Similarly, giving just one element (which is not a list) gives an error:
sage: gcd(3)
Traceback (most recent call last):
...
TypeError: 'sage.rings.integer.Integer' object is not iterable
By convention, the gcd of the empty list is (the integer) 0:
sage: gcd([])
0
sage: type(gcd([]))
<type 'sage.rings.integer.Integer'>
TESTS:
The following shows that indeed coercion takes place before computing the gcd. This behaviour was introduced in trac ticket #10771:
sage: R.<x>=QQ[]
sage: S.<x>=ZZ[]
sage: p = S.random_element(degree=(0,10))
sage: q = R.random_element(degree=(0,10))
sage: parent(gcd(1/p,q))
Fraction Field of Univariate Polynomial Ring in x over Rational Field
sage: parent(gcd([1/p,q]))
Fraction Field of Univariate Polynomial Ring in x over Rational Field
Make sure we try QQ and not merely ZZ (trac ticket #13014):
sage: bool(gcd(2/5, 3/7) == gcd(SR(2/5), SR(3/7)))
True
Make sure that the gcd of Expressions stays symbolic:
sage: parent(gcd(2, 4))
Integer Ring
sage: parent(gcd(SR(2), 4))
Symbolic Ring
sage: parent(gcd(2, SR(4)))
Symbolic Ring
sage: parent(gcd(SR(2), SR(4)))
Symbolic Ring
Verify that objects without gcd methods but which can’t be coerced to ZZ or QQ raise an error:
sage: F.<a,b> = FreeMonoid(2)
sage: gcd(a,b)
Traceback (most recent call last):
...
TypeError: unable to find gcd
The least common multiple of a and b, or if a is a list and b is omitted the least common multiple of all elements of a.
Note that LCM is an alias for lcm.
INPUT:
OUTPUT:
First, the given elements are coerced into a common parent. Then, their least common multiple in that parent is returned.
EXAMPLES:
sage: lcm(97,100)
9700
sage: LCM(97,100)
9700
sage: LCM(0,2)
0
sage: LCM(-3,-5)
15
sage: LCM([1,2,3,4,5])
60
sage: v = LCM(range(1,10000)) # *very* fast!
sage: len(str(v))
4349
TESTS:
The following tests against a bug that was fixed in trac ticket #10771:
sage: lcm(4/1,2)
4
The following shows that indeed coercion takes place before computing the least common multiple:
sage: R.<x>=QQ[]
sage: S.<x>=ZZ[]
sage: p = S.random_element(degree=(0,5))
sage: q = R.random_element(degree=(0,5))
sage: parent(lcm([1/p,q]))
Fraction Field of Univariate Polynomial Ring in x over Rational Field
Make sure we try QQ and not merely ZZ (trac ticket #13014):
sage: bool(lcm(2/5, 3/7) == lcm(SR(2/5), SR(3/7)))
True
Make sure that the lcm of Expressions stays symbolic:
sage: parent(lcm(2, 4))
Integer Ring
sage: parent(lcm(SR(2), 4))
Symbolic Ring
sage: parent(lcm(2, SR(4)))
Symbolic Ring
sage: parent(lcm(SR(2), SR(4)))
Symbolic Ring
Verify that objects without lcm methods but which can’t be coerced to ZZ or QQ raise an error:
sage: F.<a,b> = FreeMonoid(2)
sage: lcm(a,b)
Traceback (most recent call last):
...
TypeError: unable to find lcm
Check rational and integers (trac ticket #17852):
sage: lcm(1/2, 4)
4
sage: lcm(4, 1/2)
4
Returns the value of the Moebius function of abs(n), where n is an integer.
DEFINITION: \(\mu(n)\) is 0 if \(n\) is not square free, and otherwise equals \((-1)^r\), where \(n\) has \(r\) distinct prime factors.
For simplicity, if \(n=0\) we define \(\mu(n) = 0\).
IMPLEMENTATION: Factors or - for integers - uses the PARI C library.
INPUT:
OUTPUT: 0, 1, or -1
EXAMPLES:
sage: moebius(-5)
-1
sage: moebius(9)
0
sage: moebius(12)
0
sage: moebius(-35)
1
sage: moebius(-1)
1
sage: moebius(7)
-1
sage: moebius(0) # potentially nonstandard!
0
The moebius function even makes sense for non-integer inputs.
sage: x = GF(7)['x'].0
sage: moebius(x+2)
-1
Plot the Moebius function.
INPUT:
EXAMPLES:
sage: p = Moebius().plot()
sage: p.ymax()
1.0
Return the Moebius function evaluated at the given range of values, i.e., the image of the list range(start, stop, step) under the Mobius function.
This is much faster than directly computing all these values with a list comprehension.
EXAMPLES:
sage: v = moebius.range(-10,10); v
[1, 0, 0, -1, 1, -1, 0, -1, -1, 1, 0, 1, -1, -1, 0, -1, 1, -1, 0, 0]
sage: v == [moebius(n) for n in range(-10,10)]
True
sage: v = moebius.range(-1000, 2000, 4)
sage: v == [moebius(n) for n in range(-1000,2000, 4)]
True
Return the sum of the k-th powers of the divisors of n.
INPUT:
OUTPUT: integer
EXAMPLES:
sage: sigma(5)
6
sage: sigma(5,2)
26
The sigma function also has a special plotting method.
sage: P = plot(sigma, 1, 100)
This method also works with k-th powers.
sage: P = plot(sigma, 1, 100, k=2)
AUTHORS:
TESTS:
sage: sigma(100,4)
106811523
sage: sigma(factorial(100),3).mod(144169)
3672
sage: sigma(factorial(150),12).mod(691)
176
sage: RR(sigma(factorial(133),20))
2.80414775675747e4523
sage: sigma(factorial(100),0)
39001250856960000
sage: sigma(factorial(41),1)
229199532273029988767733858700732906511758707916800
Plot the sigma (sum of k-th powers of divisors) function.
INPUT:
EXAMPLES:
sage: p = Sigma().plot()
sage: p.ymax()
124.0
Return a triple (g,s,t) such that \(g = s\cdot a+t\cdot b = \gcd(a,b)\).
Note
One exception is if \(a\) and \(b\) are not in a principal ideal domain (see Wikipedia article Principal_ideal_domain), e.g., they are both polynomials over the integers. Then this function can’t in general return (g,s,t) as above, since they need not exist. Instead, over the integers, we first multiply \(g\) by a divisor of the resultant of \(a/g\) and \(b/g\), up to sign.
INPUT:
OUTPUT:
Note
There is no guarantee that the returned cofactors (s and t) are minimal.
EXAMPLES:
sage: xgcd(56, 44)
(4, 4, -5)
sage: 4*56 + (-5)*44
4
sage: g, a, b = xgcd(5/1, 7/1); g, a, b
(1, 3, -2)
sage: a*(5/1) + b*(7/1) == g
True
sage: x = polygen(QQ)
sage: xgcd(x^3 - 1, x^2 - 1)
(x - 1, 1, -x)
sage: K.<g> = NumberField(x^2-3)
sage: g.xgcd(g+2)
(1, 1/3*g, 0)
sage: R.<a,b> = K[]
sage: S.<y> = R.fraction_field()[]
sage: xgcd(y^2, a*y+b)
(1, a^2/b^2, ((-a)/b^2)*y + 1/b)
sage: xgcd((b+g)*y^2, (a+g)*y+b)
(1, (a^2 + (2*g)*a + 3)/(b^3 + (g)*b^2), ((-a + (-g))/b^2)*y + 1/b)
Here is an example of a xgcd for two polynomials over the integers, where the linear combination is not the gcd but the gcd multiplied by the resultant:
sage: R.<x> = ZZ[]
sage: gcd(2*x*(x-1), x^2)
x
sage: xgcd(2*x*(x-1), x^2)
(2*x, -1, 2)
sage: (2*(x-1)).resultant(x)
2
Returns a polynomial of degree at most \(degree\) which is approximately satisfied by the number \(z\). Note that the returned polynomial need not be irreducible, and indeed usually won’t be if \(z\) is a good approximation to an algebraic number of degree less than \(degree\).
You can specify the number of known bits or digits of \(z\) with known_bits=k or known_digits=k. PARI is then told to compute the result using \(0.8k\) of these bits/digits. Or, you can specify the precision to use directly with use_bits=k or use_digits=k. If none of these are specified, then the precision is taken from the input value.
A height bound may be specified to indicate the maximum coefficient size of the returned polynomial; if a sufficiently small polynomial is not found, then None will be returned. If proof=True then the result is returned only if it can be proved correct (i.e. the only possible minimal polynomial satisfying the height bound, or no such polynomial exists). Otherwise a ValueError is raised indicating that higher precision is required.
ALGORITHM: Uses LLL for real/complex inputs, PARI C-library algdep command otherwise.
Note that algebraic_dependency is a synonym for algdep.
INPUT:
z - real, complex, or \(p\)-adic number
degree - an integer
coefficient size for the returned polynomial
proof - a boolean (default: False), requires height_bound to be set
EXAMPLES:
sage: algdep(1.888888888888888, 1)
9*x - 17
sage: algdep(0.12121212121212,1)
33*x - 4
sage: algdep(sqrt(2),2)
x^2 - 2
This example involves a complex number:
sage: z = (1/2)*(1 + RDF(sqrt(3)) *CC.0); z
0.500000000000000 + 0.866025403784439*I
sage: p = algdep(z, 6); p
x^3 + 1
sage: p.factor()
(x + 1) * (x^2 - x + 1)
sage: z^2 - z + 1 # abs tol 2e-16
0.000000000000000
This example involves a \(p\)-adic number:
sage: K = Qp(3, print_mode = 'series')
sage: a = K(7/19); a
1 + 2*3 + 3^2 + 3^3 + 2*3^4 + 2*3^5 + 3^8 + 2*3^9 + 3^11 + 3^12 + 2*3^15 + 2*3^16 + 3^17 + 2*3^19 + O(3^20)
sage: algdep(a, 1)
19*x - 7
These examples show the importance of proper precision control. We compute a 200-bit approximation to \(sqrt(2)\) which is wrong in the 33’rd bit:
sage: z = sqrt(RealField(200)(2)) + (1/2)^33
sage: p = algdep(z, 4); p
227004321085*x^4 - 216947902586*x^3 - 99411220986*x^2 + 82234881648*x - 211871195088
sage: factor(p)
227004321085*x^4 - 216947902586*x^3 - 99411220986*x^2 + 82234881648*x - 211871195088
sage: algdep(z, 4, known_bits=32)
x^2 - 2
sage: algdep(z, 4, known_digits=10)
x^2 - 2
sage: algdep(z, 4, use_bits=25)
x^2 - 2
sage: algdep(z, 4, use_digits=8)
x^2 - 2
Using the height_bound and proof parameters, we can see that \(pi\) is not the root of an integer polynomial of degree at most 5 and coefficients bounded above by 10:
sage: algdep(pi.n(), 5, height_bound=10, proof=True) is None
True
For stronger results, we need more precicion:
sage: algdep(pi.n(), 5, height_bound=100, proof=True) is None
Traceback (most recent call last):
...
ValueError: insufficient precision for non-existence proof
sage: algdep(pi.n(200), 5, height_bound=100, proof=True) is None
True
sage: algdep(pi.n(), 10, height_bound=10, proof=True) is None
Traceback (most recent call last):
...
ValueError: insufficient precision for non-existence proof
sage: algdep(pi.n(200), 10, height_bound=10, proof=True) is None
True
We can also use proof=True to get positive results:
sage: a = sqrt(2) + sqrt(3) + sqrt(5)
sage: algdep(a.n(), 8, height_bound=1000, proof=True)
Traceback (most recent call last):
...
ValueError: insufficient precision for uniqueness proof
sage: f = algdep(a.n(1000), 8, height_bound=1000, proof=True); f
x^8 - 40*x^6 + 352*x^4 - 960*x^2 + 576
sage: f(a).expand()
0
TESTS:
sage: algdep(complex("1+2j"), 4)
x^2 - 2*x + 5
Returns a polynomial of degree at most \(degree\) which is approximately satisfied by the number \(z\). Note that the returned polynomial need not be irreducible, and indeed usually won’t be if \(z\) is a good approximation to an algebraic number of degree less than \(degree\).
You can specify the number of known bits or digits of \(z\) with known_bits=k or known_digits=k. PARI is then told to compute the result using \(0.8k\) of these bits/digits. Or, you can specify the precision to use directly with use_bits=k or use_digits=k. If none of these are specified, then the precision is taken from the input value.
A height bound may be specified to indicate the maximum coefficient size of the returned polynomial; if a sufficiently small polynomial is not found, then None will be returned. If proof=True then the result is returned only if it can be proved correct (i.e. the only possible minimal polynomial satisfying the height bound, or no such polynomial exists). Otherwise a ValueError is raised indicating that higher precision is required.
ALGORITHM: Uses LLL for real/complex inputs, PARI C-library algdep command otherwise.
Note that algebraic_dependency is a synonym for algdep.
INPUT:
z - real, complex, or \(p\)-adic number
degree - an integer
coefficient size for the returned polynomial
proof - a boolean (default: False), requires height_bound to be set
EXAMPLES:
sage: algdep(1.888888888888888, 1)
9*x - 17
sage: algdep(0.12121212121212,1)
33*x - 4
sage: algdep(sqrt(2),2)
x^2 - 2
This example involves a complex number:
sage: z = (1/2)*(1 + RDF(sqrt(3)) *CC.0); z
0.500000000000000 + 0.866025403784439*I
sage: p = algdep(z, 6); p
x^3 + 1
sage: p.factor()
(x + 1) * (x^2 - x + 1)
sage: z^2 - z + 1 # abs tol 2e-16
0.000000000000000
This example involves a \(p\)-adic number:
sage: K = Qp(3, print_mode = 'series')
sage: a = K(7/19); a
1 + 2*3 + 3^2 + 3^3 + 2*3^4 + 2*3^5 + 3^8 + 2*3^9 + 3^11 + 3^12 + 2*3^15 + 2*3^16 + 3^17 + 2*3^19 + O(3^20)
sage: algdep(a, 1)
19*x - 7
These examples show the importance of proper precision control. We compute a 200-bit approximation to \(sqrt(2)\) which is wrong in the 33’rd bit:
sage: z = sqrt(RealField(200)(2)) + (1/2)^33
sage: p = algdep(z, 4); p
227004321085*x^4 - 216947902586*x^3 - 99411220986*x^2 + 82234881648*x - 211871195088
sage: factor(p)
227004321085*x^4 - 216947902586*x^3 - 99411220986*x^2 + 82234881648*x - 211871195088
sage: algdep(z, 4, known_bits=32)
x^2 - 2
sage: algdep(z, 4, known_digits=10)
x^2 - 2
sage: algdep(z, 4, use_bits=25)
x^2 - 2
sage: algdep(z, 4, use_digits=8)
x^2 - 2
Using the height_bound and proof parameters, we can see that \(pi\) is not the root of an integer polynomial of degree at most 5 and coefficients bounded above by 10:
sage: algdep(pi.n(), 5, height_bound=10, proof=True) is None
True
For stronger results, we need more precicion:
sage: algdep(pi.n(), 5, height_bound=100, proof=True) is None
Traceback (most recent call last):
...
ValueError: insufficient precision for non-existence proof
sage: algdep(pi.n(200), 5, height_bound=100, proof=True) is None
True
sage: algdep(pi.n(), 10, height_bound=10, proof=True) is None
Traceback (most recent call last):
...
ValueError: insufficient precision for non-existence proof
sage: algdep(pi.n(200), 10, height_bound=10, proof=True) is None
True
We can also use proof=True to get positive results:
sage: a = sqrt(2) + sqrt(3) + sqrt(5)
sage: algdep(a.n(), 8, height_bound=1000, proof=True)
Traceback (most recent call last):
...
ValueError: insufficient precision for uniqueness proof
sage: f = algdep(a.n(1000), 8, height_bound=1000, proof=True); f
x^8 - 40*x^6 + 352*x^4 - 960*x^2 + 576
sage: f(a).expand()
0
TESTS:
sage: algdep(complex("1+2j"), 4)
x^2 - 2*x + 5
Return the n-th Bernoulli number, as a rational number.
INPUT:
EXAMPLES:
sage: bernoulli(12)
-691/2730
sage: bernoulli(50)
495057205241079648212477525/66
We demonstrate each of the alternative algorithms:
sage: bernoulli(12, algorithm='flint')
-691/2730
sage: bernoulli(12, algorithm='gap')
-691/2730
sage: bernoulli(12, algorithm='gp')
-691/2730
sage: bernoulli(12, algorithm='magma') # optional - magma
-691/2730
sage: bernoulli(12, algorithm='pari')
-691/2730
sage: bernoulli(12, algorithm='bernmm')
-691/2730
sage: bernoulli(12, algorithm='bernmm', num_threads=4)
-691/2730
TESTS:
sage: algs = ['gap','gp','pari','bernmm','flint']
sage: test_list = [ZZ.random_element(2, 2255) for _ in range(500)]
sage: vals = [[bernoulli(i,algorithm = j) for j in algs] for i in test_list] # long time (up to 21s on sage.math, 2011)
sage: union([len(union(x))==1 for x in vals]) # long time (depends on previous line)
[True]
sage: algs = ['gp','pari','bernmm']
sage: test_list = [ZZ.random_element(2256, 5000) for _ in range(500)]
sage: vals = [[bernoulli(i,algorithm = j) for j in algs] for i in test_list] # long time (up to 30s on sage.math, 2011)
sage: union([len(union(x))==1 for x in vals]) # long time (depends on previous line)
[True]
AUTHOR:
Return the binomial coefficient
which is defined for \(m \in \ZZ\) and any \(x\). We extend this definition to include cases when \(x-m\) is an integer but \(m\) is not by
If \(m < 0\), return \(0\).
INPUT:
OUTPUT: number or symbolic expression (if input is symbolic)
EXAMPLES:
sage: from sage.rings.arith import binomial
sage: binomial(5,2)
10
sage: binomial(2,0)
1
sage: binomial(1/2, 0)
1
sage: binomial(3,-1)
0
sage: binomial(20,10)
184756
sage: binomial(-2, 5)
-6
sage: binomial(-5, -2)
0
sage: binomial(RealField()('2.5'), 2)
1.87500000000000
sage: n=var('n'); binomial(n,2)
1/2*(n - 1)*n
sage: n=var('n'); binomial(n,n)
1
sage: n=var('n'); binomial(n,n-1)
n
sage: binomial(2^100, 2^100)
1
sage: x = polygen(ZZ)
sage: binomial(x, 3)
1/6*x^3 - 1/2*x^2 + 1/3*x
sage: binomial(x, x-3)
1/6*x^3 - 1/2*x^2 + 1/3*x
If \(x \in \ZZ\), there is an optional ‘algorithm’ parameter, which can be ‘mpir’ (faster for small values) or ‘pari’ (faster for large values):
sage: a = binomial(100, 45, algorithm='mpir')
sage: b = binomial(100, 45, algorithm='pari')
sage: a == b
True
TESTS:
We test that certain binomials are very fast (this should be instant) – see trac ticket #3309:
sage: a = binomial(RR(1140000.78), 23310000)
We test conversion of arguments to Integers – see trac ticket #6870:
sage: binomial(1/2,1/1)
1/2
sage: binomial(10^20+1/1,10^20)
100000000000000000001
sage: binomial(SR(10**7),10**7)
1
sage: binomial(3/2,SR(1/1))
3/2
Some floating point cases – see trac ticket #7562, trac ticket #9633, and trac ticket #12448:
sage: binomial(1.,3)
0.000000000000000
sage: binomial(-2.,3)
-4.00000000000000
sage: binomial(0.5r, 5)
0.02734375
sage: a = binomial(float(1001), float(1)); a
1001.0
sage: type(a)
<type 'float'>
sage: binomial(float(1000), 1001)
0.0
Test more output types:
sage: type(binomial(5r, 2))
<type 'int'>
sage: type(binomial(5, 2r))
<type 'sage.rings.integer.Integer'>
sage: type(binomial(5.0r, 2))
<type 'float'>
sage: type(binomial(5/1, 2))
<type 'sage.rings.rational.Rational'>
sage: R = Integers(11)
sage: b = binomial(R(7), R(3))
sage: b
2
sage: b.parent()
Ring of integers modulo 11
Test symbolic and uni/multivariate polynomials:
sage: x = polygen(ZZ)
sage: binomial(x, 3)
1/6*x^3 - 1/2*x^2 + 1/3*x
sage: binomial(x, 3).parent()
Univariate Polynomial Ring in x over Rational Field
sage: K.<x,y> = Integers(7)[]
sage: binomial(y,3)
-y^3 + 3*y^2 - 2*y
sage: binomial(y,3).parent()
Multivariate Polynomial Ring in x, y over Ring of integers modulo 7
sage: n = var('n')
sage: binomial(n,2)
1/2*(n - 1)*n
Invalid inputs:
sage: x = polygen(ZZ)
sage: binomial(x, x^2)
Traceback (most recent call last):
...
TypeError: either m or x-m must be an integer
sage: k, i = var('k,i')
sage: binomial(k,i)
Traceback (most recent call last):
...
TypeError: either m or x-m must be an integer
sage: R6 = Zmod(6)
sage: binomial(R6(5), 2)
Traceback (most recent call last):
...
ZeroDivisionError: factorial(2) not invertible in Ring of integers modulo 6
sage: R7 = Zmod(7)
sage: binomial(R7(10), 7)
Traceback (most recent call last):
...
ZeroDivisionError: factorial(7) not invertible in Ring of integers modulo 7
The last two examples failed to execute since \(2!\) and \(7!\) are respectively not invertible in \(\ZZ/6\ZZ\) and \(\ZZ/7\ZZ\). One can check that there is no well defined value for that binomial coefficient in the quotient:
sage: R6(binomial(5,2))
4
sage: R6(binomial(5+6,2))
1
sage: R7(binomial(3, 7))
0
sage: R7(binomial(10, 7))
1
sage: R7(binomial(17, 7))
2
For symbolic manipulation, you should use the function binomial() from the module sage.functions.other:
sage: from sage.functions.other import binomial
sage: binomial(k, i)
binomial(k, i)
Return a dictionary containing pairs \(\{(k_1,k_2) : C_{k,n}\}\) where \(C_{k_n}\) are binomial coefficients and \(n = k_1 + k_2\).
INPUT:
OUTPUT: dict
EXAMPLES:
sage: sorted(binomial_coefficients(3).items())
[((0, 3), 1), ((1, 2), 3), ((2, 1), 3), ((3, 0), 1)]
Notice the coefficients above are the same as below:
sage: R.<x,y> = QQ[]
sage: (x+y)^3
x^3 + 3*x^2*y + 3*x*y^2 + y^3
AUTHORS:
Function returns the continuant of the sequence \(v\) (list or tuple).
Definition: see Graham, Knuth and Patashnik, Concrete Mathematics, section 6.7: Continuants. The continuant is defined by
If n = None or n > len(v) the default n = len(v) is used.
INPUT:
OUTPUT: element of ring (integer, polynomial, etcetera).
EXAMPLES:
sage: continuant([1,2,3])
10
sage: p = continuant([2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10])
sage: q = continuant([1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10])
sage: p/q
517656/190435
sage: continued_fraction([2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10]).convergent(14)
517656/190435
sage: x = PolynomialRing(RationalField(),'x',5).gens()
sage: continuant(x)
x0*x1*x2*x3*x4 + x0*x1*x2 + x0*x1*x4 + x0*x3*x4 + x2*x3*x4 + x0 + x2 + x4
sage: continuant(x, 3)
x0*x1*x2 + x0 + x2
sage: continuant(x,2)
x0*x1 + 1
We verify the identity
for \(n = 6\) using polynomial arithmetic:
sage: z = QQ['z'].0
sage: continuant((z,z,z,z,z,z,z,z,z,z,z,z,z,z,z),6)
z^6 + 5*z^4 + 6*z^2 + 1
sage: continuant(9)
Traceback (most recent call last):
...
TypeError: object of type 'sage.rings.integer.Integer' has no len()
AUTHORS:
Returns a solution to a Chinese Remainder Theorem problem.
INPUT:
OUTPUT:
If m, n are not None, returns a solution \(x\) to the simultaneous congruences \(x\equiv a \bmod m\) and \(x\equiv b \bmod n\), if one exists. By the Chinese Remainder Theorem, a solution to the simultaneous congruences exists if and only if \(a\equiv b\pmod{\gcd(m,n)}\). The solution \(x\) is only well-defined modulo \(\text{lcm}(m,n)\).
If a and b are lists, returns a simultaneous solution to the congruences \(x\equiv a_i\pmod{b_i}\), if one exists.
See also
EXAMPLES:
Using crt by giving it pairs of residues and moduli:
sage: crt(2, 1, 3, 5)
11
sage: crt(13, 20, 100, 301)
28013
sage: crt([2, 1], [3, 5])
11
sage: crt([13, 20], [100, 301])
28013
You can also use upper case:
sage: c = CRT(2,3, 3, 5); c
8
sage: c % 3 == 2
True
sage: c % 5 == 3
True
Note that this also works for polynomial rings:
sage: K.<a> = NumberField(x^3 - 7)
sage: R.<y> = K[]
sage: f = y^2 + 3
sage: g = y^3 - 5
sage: CRT(1,3,f,g)
-3/26*y^4 + 5/26*y^3 + 15/26*y + 53/26
sage: CRT(1,a,f,g)
(-3/52*a + 3/52)*y^4 + (5/52*a - 5/52)*y^3 + (15/52*a - 15/52)*y + 27/52*a + 25/52
You can also do this for any number of moduli:
sage: K.<a> = NumberField(x^3 - 7)
sage: R.<x> = K[]
sage: CRT([], [])
0
sage: CRT([a], [x])
a
sage: f = x^2 + 3
sage: g = x^3 - 5
sage: h = x^5 + x^2 - 9
sage: k = CRT([1, a, 3], [f, g, h]); k
(127/26988*a - 5807/386828)*x^9 + (45/8996*a - 33677/1160484)*x^8 + (2/173*a - 6/173)*x^7 + (133/6747*a - 5373/96707)*x^6 + (-6/2249*a + 18584/290121)*x^5 + (-277/8996*a + 38847/386828)*x^4 + (-135/4498*a + 42673/193414)*x^3 + (-1005/8996*a + 470245/1160484)*x^2 + (-1215/8996*a + 141165/386828)*x + 621/8996*a + 836445/386828
sage: k.mod(f)
1
sage: k.mod(g)
a
sage: k.mod(h)
3
If the moduli are not coprime, a solution may not exist:
sage: crt(4,8,8,12)
20
sage: crt(4,6,8,12)
Traceback (most recent call last):
...
ValueError: No solution to crt problem since gcd(8,12) does not divide 4-6
sage: x = polygen(QQ)
sage: crt(2,3,x-1,x+1)
-1/2*x + 5/2
sage: crt(2,x,x^2-1,x^2+1)
-1/2*x^3 + x^2 + 1/2*x + 1
sage: crt(2,x,x^2-1,x^3-1)
Traceback (most recent call last):
...
ValueError: No solution to crt problem since gcd(x^2 - 1,x^3 - 1) does not divide 2-x
sage: crt(int(2), int(3), int(7), int(11))
58
Return the Dedekind sum \(s(p,q)\) defined for integers \(p\), \(q\) as
where
Warning
Caution is required as the Dedekind sum sometimes depends on the algorithm or is left undefined when \(p\) and \(q\) are not coprime.
INPUT:
OUTPUT: a rational number
EXAMPLES:
Several small values:
sage: for q in range(10): print [dedekind_sum(p,q) for p in range(q+1)]
[0]
[0, 0]
[0, 0, 0]
[0, 1/18, -1/18, 0]
[0, 1/8, 0, -1/8, 0]
[0, 1/5, 0, 0, -1/5, 0]
[0, 5/18, 1/18, 0, -1/18, -5/18, 0]
[0, 5/14, 1/14, -1/14, 1/14, -1/14, -5/14, 0]
[0, 7/16, 1/8, 1/16, 0, -1/16, -1/8, -7/16, 0]
[0, 14/27, 4/27, 1/18, -4/27, 4/27, -1/18, -4/27, -14/27, 0]
Check relations for restricted arguments:
sage: q = 23; dedekind_sum(1, q); (q-1)*(q-2)/(12*q)
77/46
77/46
sage: p, q = 100, 723 # must be coprime
sage: dedekind_sum(p, q) + dedekind_sum(q, p)
31583/86760
sage: -1/4 + (p/q + q/p + 1/(p*q))/12
31583/86760
We check that evaluation works with large input:
sage: dedekind_sum(3^54 - 1, 2^93 + 1)
459340694971839990630374299870/29710560942849126597578981379
sage: dedekind_sum(3^54 - 1, 2^93 + 1, algorithm='pari')
459340694971839990630374299870/29710560942849126597578981379
We check consistency of the results:
sage: dedekind_sum(5, 7, algorithm='default')
-1/14
sage: dedekind_sum(5, 7, algorithm='flint')
-1/14
sage: dedekind_sum(5, 7, algorithm='pari')
-1/14
sage: dedekind_sum(6, 8, algorithm='default')
-1/8
sage: dedekind_sum(6, 8, algorithm='flint')
-1/8
sage: dedekind_sum(6, 8, algorithm='pari')
-1/8
REFERENCES:
| [Apostol] | T. Apostol, Modular functions and Dirichlet series in number theory, Springer, 1997 (2nd ed), section 3.7–3.9. |
Returns the \(n\) successive differences of the elements in \(lis\).
EXAMPLES:
sage: differences(prime_range(50))
[1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4]
sage: differences([i^2 for i in range(1,11)])
[3, 5, 7, 9, 11, 13, 15, 17, 19]
sage: differences([i^3 + 3*i for i in range(1,21)])
[10, 22, 40, 64, 94, 130, 172, 220, 274, 334, 400, 472, 550, 634, 724, 820, 922, 1030, 1144]
sage: differences([i^3 - i^2 for i in range(1,21)], 2)
[10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76, 82, 88, 94, 100, 106, 112]
sage: differences([p - i^2 for i, p in enumerate(prime_range(50))], 3)
[-1, 2, -4, 4, -4, 4, 0, -6, 8, -6, 0, 4]
AUTHORS:
Returns a list of all positive integer divisors of the nonzero integer n.
INPUT:
EXAMPLES:
sage: divisors(-3)
[1, 3]
sage: divisors(6)
[1, 2, 3, 6]
sage: divisors(28)
[1, 2, 4, 7, 14, 28]
sage: divisors(2^5)
[1, 2, 4, 8, 16, 32]
sage: divisors(100)
[1, 2, 4, 5, 10, 20, 25, 50, 100]
sage: divisors(1)
[1]
sage: divisors(0)
Traceback (most recent call last):
...
ValueError: n must be nonzero
sage: divisors(2^3 * 3^2 * 17)
[1, 2, 3, 4, 6, 8, 9, 12, 17, 18, 24, 34, 36, 51, 68, 72, 102, 136, 153, 204, 306, 408, 612, 1224]
This function works whenever one has unique factorization:
sage: K.<a> = QuadraticField(7)
sage: divisors(K.ideal(7))
[Fractional ideal (1), Fractional ideal (a), Fractional ideal (7)]
sage: divisors(K.ideal(3))
[Fractional ideal (1), Fractional ideal (3), Fractional ideal (-a + 2), Fractional ideal (-a - 2)]
sage: divisors(K.ideal(35))
[Fractional ideal (1), Fractional ideal (5), Fractional ideal (a), Fractional ideal (7), Fractional ideal (5*a), Fractional ideal (35)]
TESTS:
sage: divisors(int(300))
[1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 25, 30, 50, 60, 75, 100, 150, 300]
Return a list of the primes \(\leq n\).
This is extremely slow and is for educational purposes only.
INPUT:
OUTPUT:
EXAMPLES:
sage: eratosthenes(3)
[2, 3]
sage: eratosthenes(50)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
sage: len(eratosthenes(100))
25
sage: eratosthenes(213) == prime_range(213)
True
Returns the factorization of n. The result depends on the type of n.
If n is an integer, returns the factorization as an object of type Factorization.
If n is not an integer, n.factor(proof=proof, **kwds) gets called. See n.factor?? for more documentation in this case.
Warning
This means that applying factor to an integer result of a symbolic computation will not factor the integer, because it is considered as an element of a larger symbolic ring.
EXAMPLE:
sage: f(n)=n^2
sage: is_prime(f(3))
False
sage: factor(f(3))
9
INPUT:
OUTPUT:
The qsieve and ecm commands give access to highly optimized implementations of algorithms for doing certain integer factorization problems. These implementations are not used by the generic factor command, which currently just calls PARI (note that PARI also implements sieve and ecm algorithms, but they aren’t as optimized). Thus you might consider using them instead for certain numbers.
The factorization returned is an element of the class Factorization; see Factorization?? for more details, and examples below for usage. A Factorization contains both the unit factor (+1 or -1) and a sorted list of (prime, exponent) pairs.
The factorization displays in pretty-print format but it is easy to obtain access to the (prime,exponent) pairs and the unit, to recover the number from its factorization, and even to multiply two factorizations. See examples below.
EXAMPLES:
sage: factor(500)
2^2 * 5^3
sage: factor(-20)
-1 * 2^2 * 5
sage: f=factor(-20)
sage: list(f)
[(2, 2), (5, 1)]
sage: f.unit()
-1
sage: f.value()
-20
sage: factor( -next_prime(10^2) * next_prime(10^7) )
-1 * 101 * 10000019
sage: factor(-500, algorithm='kash') # optional - kash
-1 * 2^2 * 5^3
sage: factor(-500, algorithm='magma') # optional - magma
-1 * 2^2 * 5^3
sage: factor(0)
Traceback (most recent call last):
...
ArithmeticError: Prime factorization of 0 not defined.
sage: factor(1)
1
sage: factor(-1)
-1
sage: factor(2^(2^7)+1)
59649589127497217 * 5704689200685129054721
Sage calls PARI’s factor, which has proof False by default. Sage has a global proof flag, set to True by default (see sage.structure.proof.proof, or proof.[tab]). To override the default, call this function with proof=False.
sage: factor(3^89-1, proof=False)
2 * 179 * 1611479891519807 * 5042939439565996049162197
sage: factor(2^197 + 1) # long time (2s)
3 * 197002597249 * 1348959352853811313 * 251951573867253012259144010843
Any object which has a factor method can be factored like this:
sage: K.<i> = QuadraticField(-1)
sage: factor(122 - 454*i)
(-3*i - 2) * (-i - 2)^3 * (i + 1)^3 * (i + 4)
To access the data in a factorization:
sage: f = factor(420); f
2^2 * 3 * 5 * 7
sage: [x for x in f]
[(2, 2), (3, 1), (5, 1), (7, 1)]
sage: [p for p,e in f]
[2, 3, 5, 7]
sage: [e for p,e in f]
[2, 1, 1, 1]
sage: [p^e for p,e in f]
[4, 3, 5, 7]
Compute the factorial of \(n\), which is the product \(1\cdot 2\cdot 3 \cdots (n-1)\cdot n\).
INPUT:
OUTPUT: an integer
EXAMPLES:
sage: from sage.rings.arith import factorial
sage: factorial(0)
1
sage: factorial(4)
24
sage: factorial(10)
3628800
sage: factorial(1) == factorial(0)
True
sage: factorial(6) == 6*5*4*3*2
True
sage: factorial(1) == factorial(0)
True
sage: factorial(71) == 71* factorial(70)
True
sage: factorial(-32)
Traceback (most recent call last):
...
ValueError: factorial -- must be nonnegative
PERFORMANCE: This discussion is valid as of April 2006. All timings below are on a Pentium Core Duo 2Ghz MacBook Pro running Linux with a 2.6.16.1 kernel.
Returns the falling factorial \((x)_a\).
The notation in the literature is a mess: often \((x)_a\), but there are many other notations: GKP: Concrete Mathematics uses \(x^{\underline{a}}\).
Definition: for integer \(a \ge 0\) we have \(x(x-1) \cdots (x-a+1)\). In all other cases we use the GAMMA-function: \(\frac {\Gamma(x+1)} {\Gamma(x-a+1)}\).
INPUT:
OR
OUTPUT: the falling factorial
EXAMPLES:
sage: falling_factorial(10, 3)
720
sage: falling_factorial(10, RR('3.0'))
720.000000000000
sage: falling_factorial(10, RR('3.3'))
1310.11633396601
sage: falling_factorial(10, 10)
3628800
sage: factorial(10)
3628800
sage: a = falling_factorial(1+I, I); a
gamma(I + 2)
sage: CC(a)
0.652965496420167 + 0.343065839816545*I
sage: falling_factorial(1+I, 4)
4*I + 2
sage: falling_factorial(I, 4)
-10
sage: M = MatrixSpace(ZZ, 4, 4)
sage: A = M([1,0,1,0,1,0,1,0,1,0,10,10,1,0,1,1])
sage: falling_factorial(A, 2) # A(A - I)
[ 1 0 10 10]
[ 1 0 10 10]
[ 20 0 101 100]
[ 2 0 11 10]
sage: x = ZZ['x'].0
sage: falling_factorial(x, 4)
x^4 - 6*x^3 + 11*x^2 - 6*x
TESTS:
Check that trac ticket #14858 is fixed:
sage: falling_factorial(-4, SR(2))
20
Check that trac ticket #16770 is fixed:
sage: d = var('d')
sage: type(falling_factorial(d, 0))
<type 'sage.symbolic.expression.Expression'>
AUTHORS:
Write the integer \(n\) as a sum of four integer squares.
INPUT:
OUTPUT: a tuple \((a,b,c,d)\) of non-negative integers such that \(n = a^2 + b^2 + c^2 + d^2\) with \(a <= b <= c <= d\).
EXAMPLES:
sage: four_squares(3)
(0, 1, 1, 1)
sage: four_squares(13)
(0, 0, 2, 3)
sage: four_squares(130)
(0, 0, 3, 11)
sage: four_squares(1101011011004)
(90, 102, 1220, 1049290)
sage: four_squares(10^100-1)
(155024616290, 2612183768627, 14142135623730950488016887, 99999999999999999999999999999999999999999999999999)
sage: for i in range(2^129, 2^129+10000): # long time
....: S = four_squares(i)
....: assert sum(x^2 for x in S) == i
TESTS:
sage: for _ in xrange(100):
....: n = ZZ.random_element(2**32,2**34)
....: aa,bb,cc,dd = four_squares(n)
....: assert aa**2 + bb**2 + cc**2 + dd**2 == n
Return the discriminant of the quadratic extension \(K=Q(\sqrt{D})\), i.e. an integer d congruent to either 0 or 1, mod 4, and such that, at most, the only square dividing it is 4.
INPUT:
OUTPUT:
EXAMPLES:
sage: fundamental_discriminant(102)
408
sage: fundamental_discriminant(720)
5
sage: fundamental_discriminant(2)
8
The greatest common divisor of a and b, or if a is a list and b is omitted the greatest common divisor of all elements of a.
INPUT:
Additional keyword arguments are passed to the respectively called methods.
OUTPUT:
The given elements are first coerced into a common parent. Then, their greatest common divisor in that common parent is returned.
EXAMPLES:
sage: GCD(97,100)
1
sage: GCD(97*10^15, 19^20*97^2)
97
sage: GCD(2/3, 4/5)
2/15
sage: GCD([2,4,6,8])
2
sage: GCD(srange(0,10000,10)) # fast !!
10
Note that to take the gcd of \(n\) elements for \(n \not= 2\) you must put the elements into a list by enclosing them in [..]. Before #4988 the following wrongly returned 3 since the third parameter was just ignored:
sage: gcd(3,6,2)
Traceback (most recent call last):
...
TypeError: gcd() takes at most 2 arguments (3 given)
sage: gcd([3,6,2])
1
Similarly, giving just one element (which is not a list) gives an error:
sage: gcd(3)
Traceback (most recent call last):
...
TypeError: 'sage.rings.integer.Integer' object is not iterable
By convention, the gcd of the empty list is (the integer) 0:
sage: gcd([])
0
sage: type(gcd([]))
<type 'sage.rings.integer.Integer'>
TESTS:
The following shows that indeed coercion takes place before computing the gcd. This behaviour was introduced in trac ticket #10771:
sage: R.<x>=QQ[]
sage: S.<x>=ZZ[]
sage: p = S.random_element(degree=(0,10))
sage: q = R.random_element(degree=(0,10))
sage: parent(gcd(1/p,q))
Fraction Field of Univariate Polynomial Ring in x over Rational Field
sage: parent(gcd([1/p,q]))
Fraction Field of Univariate Polynomial Ring in x over Rational Field
Make sure we try QQ and not merely ZZ (trac ticket #13014):
sage: bool(gcd(2/5, 3/7) == gcd(SR(2/5), SR(3/7)))
True
Make sure that the gcd of Expressions stays symbolic:
sage: parent(gcd(2, 4))
Integer Ring
sage: parent(gcd(SR(2), 4))
Symbolic Ring
sage: parent(gcd(2, SR(4)))
Symbolic Ring
sage: parent(gcd(SR(2), SR(4)))
Symbolic Ring
Verify that objects without gcd methods but which can’t be coerced to ZZ or QQ raise an error:
sage: F.<a,b> = FreeMonoid(2)
sage: gcd(a,b)
Traceback (most recent call last):
...
TypeError: unable to find gcd
Return the fastest gcd function for integers of size no larger than order.
EXAMPLES:
sage: sage.rings.arith.get_gcd(4000)
<built-in method gcd_int of sage.rings.fast_arith.arith_int object at ...>
sage: sage.rings.arith.get_gcd(400000)
<built-in method gcd_longlong of sage.rings.fast_arith.arith_llong object at ...>
sage: sage.rings.arith.get_gcd(4000000000)
<function gcd at ...>
Return the fastest inverse_mod function for integers of size no larger than order.
EXAMPLES:
sage: sage.rings.arith.get_inverse_mod(6000)
<built-in method inverse_mod_int of sage.rings.fast_arith.arith_int object at ...>
sage: sage.rings.arith.get_inverse_mod(600000)
<built-in method inverse_mod_longlong of sage.rings.fast_arith.arith_llong object at ...>
sage: sage.rings.arith.get_inverse_mod(6000000000)
<function inverse_mod at ...>
This is the product of all (finite) primes where the Hilbert symbol is -1. What is the same, this is the (reduced) discriminant of the quaternion algebra \((a,b)\) over \(\QQ\).
INPUT:
OUTPUT:
EXAMPLES:
sage: hilbert_conductor(-1, -1)
2
sage: hilbert_conductor(-1, -11)
11
sage: hilbert_conductor(-2, -5)
5
sage: hilbert_conductor(-3, -17)
17
AUTHOR:
Finds a pair of integers \((a,b)\) such that hilbert_conductor(a,b) == d. The quaternion algebra \((a,b)\) over \(\QQ\) will then have (reduced) discriminant \(d\).
INPUT:
OUTPUT: pair of integers
EXAMPLES:
sage: hilbert_conductor_inverse(2)
(-1, -1)
sage: hilbert_conductor_inverse(3)
(-1, -3)
sage: hilbert_conductor_inverse(6)
(-1, 3)
sage: hilbert_conductor_inverse(30)
(-3, -10)
sage: hilbert_conductor_inverse(4)
Traceback (most recent call last):
...
ValueError: d needs to be squarefree
sage: hilbert_conductor_inverse(-1)
Traceback (most recent call last):
...
ValueError: d needs to be positive
AUTHOR:
TESTS:
sage: for i in xrange(100):
... d = ZZ.random_element(2**32).squarefree_part()
... if hilbert_conductor(*hilbert_conductor_inverse(d)) != d:
... print "hilbert_conductor_inverse failed for d =", d
Returns 1 if \(ax^2 + by^2\) \(p\)-adically represents a nonzero square, otherwise returns \(-1\). If either a or b is 0, returns 0.
INPUT:
OUTPUT: integer (0, -1, or 1)
EXAMPLES:
sage: hilbert_symbol (-1, -1, -1, algorithm='all')
-1
sage: hilbert_symbol (2,3, 5, algorithm='all')
1
sage: hilbert_symbol (4, 3, 5, algorithm='all')
1
sage: hilbert_symbol (0, 3, 5, algorithm='all')
0
sage: hilbert_symbol (-1, -1, 2, algorithm='all')
-1
sage: hilbert_symbol (1, -1, 2, algorithm='all')
1
sage: hilbert_symbol (3, -1, 2, algorithm='all')
-1
sage: hilbert_symbol(QQ(-1)/QQ(4), -1, 2) == -1
True
sage: hilbert_symbol(QQ(-1)/QQ(4), -1, 3) == 1
True
AUTHORS:
Return the ceiling of x.
EXAMPLES:
sage: integer_ceil(5.4)
6
sage: integer_ceil(x)
Traceback (most recent call last):
...
NotImplementedError: computation of ceil of x not implemented
Return the largest integer \(\leq x\).
INPUT:
OUTPUT: an Integer
EXAMPLES:
sage: integer_floor(5.4)
5
sage: integer_floor(float(5.4))
5
sage: integer_floor(-5/2)
-3
sage: integer_floor(RDF(-5/2))
-3
sage: integer_floor(x)
Traceback (most recent call last):
...
NotImplementedError: computation of floor of x not implemented
The inverse of the ring element a modulo m.
If no special inverse_mod is defined for the elements, it tries to coerce them into integers and perform the inversion there
sage: inverse_mod(7,1)
0
sage: inverse_mod(5,14)
3
sage: inverse_mod(3,-5)
2
This function returns True if and only if n is a power of 2
INPUT:
OUTPUT:
EXAMPLES:
sage: is_power_of_two(1024)
True
sage: is_power_of_two(1)
True
sage: is_power_of_two(24)
False
sage: is_power_of_two(0)
False
sage: is_power_of_two(-4)
False
Return True if \(n\) is a prime number, and False otherwise.
Use a provable primality test or a strong pseudo-primality test depending on the global arithmetic proof flag.
INPUT:
AUTHORS:
EXAMPLES:
sage: is_prime(389)
True
sage: is_prime(2000)
False
sage: is_prime(2)
True
sage: is_prime(-1)
False
sage: is_prime(1)
False
sage: is_prime(-2)
False
sage: a = 2**2048 + 981
sage: is_prime(a) # not tested - takes ~ 1min
sage: proof.arithmetic(False)
sage: is_prime(a) # instantaneous!
True
sage: proof.arithmetic(True)
Test whether n is a positive power of a prime number
This function simply calls the method Integer.is_prime_power() of Integers.
INPUT:
EXAMPLES:
sage: is_prime_power(389)
True
sage: is_prime_power(2000)
False
sage: is_prime_power(2)
True
sage: is_prime_power(1024)
True
sage: is_prime_power(1024, get_data=True)
(2, 10)
The same results can be obtained with:
sage: 389.is_prime_power()
True
sage: 2000.is_prime_power()
False
sage: 2.is_prime_power()
True
sage: 1024.is_prime_power()
True
sage: 1024.is_prime_power(get_data=True)
(2, 10)
TESTS:
sage: is_prime_power(-1)
False
sage: is_prime_power(1)
False
sage: is_prime_power(QQ(997^100))
True
sage: is_prime_power(1/2197)
Traceback (most recent call last):
...
TypeError: no conversion of this rational to integer
sage: is_prime_power("foo")
Traceback (most recent call last):
...
TypeError: unable to convert 'foo' to an integer
Test whether n is a pseudo-prime
The result is NOT proven correct - this is a pseudo-primality test!.
INPUT:
Note
We do not consider negatives of prime numbers as prime.
EXAMPLES:
sage: is_pseudoprime(389)
True
sage: is_pseudoprime(2000)
False
sage: is_pseudoprime(2)
True
sage: is_pseudoprime(-1)
False
sage: factor(-6)
-1 * 2 * 3
sage: is_pseudoprime(1)
False
sage: is_pseudoprime(-2)
False
TESTS:
Deprecation warning from trac ticket #16878:
sage: is_pseudoprime(127, flag=0)
doctest:...: DeprecationWarning: the keyword 'flag' is deprecated and no longer used
See http://trac.sagemath.org/16878 for details.
True
Test if n is a power of a pseudoprime.
The result is NOT proven correct - this IS a pseudo-primality test!. Note that a prime power is a positive power of a prime number so that 1 is not a prime power.
INPUT:
EXAMPLES:
sage: is_pseudoprime_power(389)
True
sage: is_pseudoprime_power(2000)
False
sage: is_pseudoprime_power(2)
True
sage: is_pseudoprime_power(1024)
True
sage: is_pseudoprime_power(-1)
False
sage: is_pseudoprime_power(1)
False
sage: is_pseudoprime_power(997^100)
True
Use of the get_data keyword:
sage: is_pseudoprime_power(3^1024, get_data=True)
(3, 1024)
sage: is_pseudoprime_power(2^256, get_data=True)
(2, 256)
sage: is_pseudoprime_power(31, get_data=True)
(31, 1)
sage: is_pseudoprime_power(15, get_data=True)
(15, 0)
Deprecated version of is_pseudoprime_power.
EXAMPLES:
sage: is_pseudoprime_small_power(1234)
doctest:...: DeprecationWarning: the function is_pseudoprime_small_power() is deprecated, use is_pseudoprime_power() instead.
See http://trac.sagemath.org/16878 for details.
False
sage: is_pseudoprime_small_power(3^1024, get_data=True)
[(3, 1024)]
Returns whether or not n is square, and if n is a square also returns the square root. If n is not square, also returns None.
INPUT:
OUTPUT:
EXAMPLES:
sage: is_square(2)
False
sage: is_square(4)
True
sage: is_square(2.2)
True
sage: is_square(-2.2)
False
sage: is_square(CDF(-2.2))
True
sage: is_square((x-1)^2)
True
sage: is_square(4, True)
(True, 2)
Test whether n is square free.
EXAMPLES:
sage: is_squarefree(100)
False
sage: is_squarefree(101)
True
sage: R = ZZ['x']
sage: x = R.gen()
sage: is_squarefree((x^2+x+1) * (x-2))
True
sage: is_squarefree((x-1)**2 * (x-3))
False
sage: O = ZZ[sqrt(-1)]
sage: I = O.gen(1)
sage: is_squarefree(I+1)
True
sage: is_squarefree(O(2))
False
sage: O(2).factor()
(-I) * (I + 1)^2
This method fails on domains which are not Unique Factorization Domains:
sage: O = ZZ[sqrt(-5)]
sage: a = O.gen(1)
sage: is_squarefree(a - 3)
Traceback (most recent call last):
...
ArithmeticError: non-principal ideal in factorization
The Jacobi symbol of integers a and b, where b is odd.
Note
The kronecker_symbol() command extends the Jacobi symbol to all integers b.
If
\(b = p_1^{e_1} * ... * p_r^{e_r}\)
then
\((a|b) = (a|p_1)^{e_1} ... (a|p_r)^{e_r}\)
where \((a|p_j)\) are Legendre Symbols.
INPUT:
EXAMPLES:
sage: jacobi_symbol(10,777)
-1
sage: jacobi_symbol(10,5)
0
sage: jacobi_symbol(10,2)
Traceback (most recent call last):
...
ValueError: second input must be odd, 2 is not odd
Synonym for kronecker_symbol().
The Kronecker symbol \((x|y)\).
INPUT:
OUTPUT:
EXAMPLES:
sage: kronecker(3,5)
-1
sage: kronecker(3,15)
0
sage: kronecker(2,15)
1
sage: kronecker(-2,15)
-1
sage: kronecker(2/3,5)
1
The Kronecker symbol \((x|y)\).
INPUT:
EXAMPLES:
sage: kronecker_symbol(13,21)
-1
sage: kronecker_symbol(101,4)
1
IMPLEMENTATION: Using GMP.
The least common multiple of a and b, or if a is a list and b is omitted the least common multiple of all elements of a.
Note that LCM is an alias for lcm.
INPUT:
OUTPUT:
First, the given elements are coerced into a common parent. Then, their least common multiple in that parent is returned.
EXAMPLES:
sage: lcm(97,100)
9700
sage: LCM(97,100)
9700
sage: LCM(0,2)
0
sage: LCM(-3,-5)
15
sage: LCM([1,2,3,4,5])
60
sage: v = LCM(range(1,10000)) # *very* fast!
sage: len(str(v))
4349
TESTS:
The following tests against a bug that was fixed in trac ticket #10771:
sage: lcm(4/1,2)
4
The following shows that indeed coercion takes place before computing the least common multiple:
sage: R.<x>=QQ[]
sage: S.<x>=ZZ[]
sage: p = S.random_element(degree=(0,5))
sage: q = R.random_element(degree=(0,5))
sage: parent(lcm([1/p,q]))
Fraction Field of Univariate Polynomial Ring in x over Rational Field
Make sure we try QQ and not merely ZZ (trac ticket #13014):
sage: bool(lcm(2/5, 3/7) == lcm(SR(2/5), SR(3/7)))
True
Make sure that the lcm of Expressions stays symbolic:
sage: parent(lcm(2, 4))
Integer Ring
sage: parent(lcm(SR(2), 4))
Symbolic Ring
sage: parent(lcm(2, SR(4)))
Symbolic Ring
sage: parent(lcm(SR(2), SR(4)))
Symbolic Ring
Verify that objects without lcm methods but which can’t be coerced to ZZ or QQ raise an error:
sage: F.<a,b> = FreeMonoid(2)
sage: lcm(a,b)
Traceback (most recent call last):
...
TypeError: unable to find lcm
Check rational and integers (trac ticket #17852):
sage: lcm(1/2, 4)
4
sage: lcm(4, 1/2)
4
The Legendre symbol \((x|p)\), for \(p\) prime.
Note
The kronecker_symbol() command extends the Legendre symbol to composite moduli and \(p=2\).
INPUT:
EXAMPLES:
sage: legendre_symbol(2,3)
-1
sage: legendre_symbol(1,3)
1
sage: legendre_symbol(1,2)
Traceback (most recent call last):
...
ValueError: p must be odd
sage: legendre_symbol(2,15)
Traceback (most recent call last):
...
ValueError: p must be a prime
sage: kronecker_symbol(2,15)
1
sage: legendre_symbol(2/3,7)
-1
Maximal Quotient Rational Reconstruction.
For research purposes only - this is pure Python, so slow.
INPUT:
OUTPUT:
Either integers \(n,d\) such that \(d>0\), \(\mathop{\mathrm{gcd}}(n,d)=1\), \(n/d=u \bmod m\), and \(T \cdot d \cdot |n| < m\), or None.
Reference: Monagan, Maximal Quotient Rational Reconstruction: An Almost Optimal Algorithm for Rational Reconstruction (page 11)
This algorithm is probabilistic.
EXAMPLES:
sage: mqrr_rational_reconstruction(21,3100,13)
(21, 1)
Return the multinomial coefficient
INPUT:
OUTPUT:
Returns the integer:
EXAMPLES:
sage: multinomial(0, 0, 2, 1, 0, 0)
3
sage: multinomial([0, 0, 2, 1, 0, 0])
3
sage: multinomial(3, 2)
10
sage: multinomial(2^30, 2, 1)
618970023101454657175683075
sage: multinomial([2^30, 2, 1])
618970023101454657175683075
AUTHORS:
Return a dictionary containing pairs \(\{(k_1, k_2, ..., k_m) : C_{k, n}\}\) where \(C_{k, n}\) are multinomial coefficients such that \(n = k_1 + k_2 + ...+ k_m\).
INPUT:
OUTPUT: dict
EXAMPLES:
sage: sorted(multinomial_coefficients(2, 5).items())
[((0, 5), 1), ((1, 4), 5), ((2, 3), 10), ((3, 2), 10), ((4, 1), 5), ((5, 0), 1)]
Notice that these are the coefficients of \((x+y)^5\):
sage: R.<x,y> = QQ[]
sage: (x+y)^5
x^5 + 5*x^4*y + 10*x^3*y^2 + 10*x^2*y^3 + 5*x*y^4 + y^5
sage: sorted(multinomial_coefficients(3, 2).items())
[((0, 0, 2), 1), ((0, 1, 1), 2), ((0, 2, 0), 1), ((1, 0, 1), 2), ((1, 1, 0), 2), ((2, 0, 0), 1)]
ALGORITHM: The algorithm we implement for computing the multinomial coefficients is based on the following result:
..math:
\binom{n}{k_1, \cdots, k_m} =
\frac{k_1+1}{n-k_1}\sum_{i=2}^m \binom{n}{k_1+1, \cdots, k_i-1, \cdots}
e.g.:
sage: k = (2, 4, 1, 0, 2, 6, 0, 0, 3, 5, 7, 1) # random value
sage: n = sum(k)
sage: s = 0
sage: for i in range(1, len(k)):
... ki = list(k)
... ki[0] += 1
... ki[i] -= 1
... s += multinomial(n, *ki)
sage: multinomial(n, *k) == (k[0] + 1) / (n - k[0]) * s
True
TESTS:
sage: multinomial_coefficients(0, 0)
{(): 1}
sage: multinomial_coefficients(0, 3)
{}
The next prime greater than the integer n. If n is prime, then this function does not return n, but the next prime after n. If the optional argument proof is False, this function only returns a pseudo-prime, as defined by the PARI nextprime function. If it is None, uses the global default (see sage.structure.proof.proof)
INPUT:
EXAMPLES:
sage: next_prime(-100)
2
sage: next_prime(1)
2
sage: next_prime(2)
3
sage: next_prime(3)
5
sage: next_prime(4)
5
Notice that the next_prime(5) is not 5 but 7.
sage: next_prime(5)
7
sage: next_prime(2004)
2011
Return the smallest prime power greater than n.
Note that if n is a prime power, then this function does not return n, but the next prime power after n.
This function just calls the method Integer.next_prime_power() of Integers.
See also
EXAMPLES:
sage: next_prime_power(1)
2
sage: next_prime_power(2)
3
sage: next_prime_power(10)
11
sage: next_prime_power(7)
8
sage: next_prime_power(99)
101
The same results can be obtained with:
sage: 1.next_prime_power()
2
sage: 2.next_prime_power()
3
sage: 10.next_prime_power()
11
Note that \(2\) is the smallest prime power:
sage: next_prime_power(-10)
2
sage: next_prime_power(0)
2
Returns the next probable prime after self, as determined by PARI.
INPUT:
EXAMPLES:
sage: next_probable_prime(-100)
2
sage: next_probable_prime(19)
23
sage: next_probable_prime(int(999999999))
1000000007
sage: next_probable_prime(2^768)
1552518092300708935148979488462502555256886017116696611139052038026050952686376886330878408828646477950487730697131073206171580044114814391444287275041181139204454976020849905550265285631598444825262999193716468750892846853816058039
Return the n-th prime number (1-indexed, so that 2 is the 1st prime.)
INPUT:
OUTPUT:
EXAMPLES:
sage: nth_prime(3)
5
sage: nth_prime(10)
29
sage: nth_prime(0)
Traceback (most recent call last):
...
ValueError: nth prime meaningless for non-positive n (=0)
TESTS:
sage: all(prime_pi(nth_prime(j)) == j for j in range(1, 1000, 10))
True
Return the number of divisors of the integer n.
INPUT:
OUTPUT:
EXAMPLES:
sage: number_of_divisors(100)
9
sage: number_of_divisors(-720)
30
The odd part of the integer \(n\). This is \(n / 2^v\), where \(v = \mathrm{valuation}(n,2)\).
EXAMPLES:
sage: odd_part(5)
5
sage: odd_part(4)
1
sage: odd_part(factorial(31))
122529844256906551386796875
The n-th power of a modulo the integer m.
EXAMPLES:
sage: power_mod(0,0,5)
Traceback (most recent call last):
...
ArithmeticError: 0^0 is undefined.
sage: power_mod(2,390,391)
285
sage: power_mod(2,-1,7)
4
sage: power_mod(11,1,7)
4
sage: R.<x> = ZZ[]
sage: power_mod(3*x, 10, 7)
4*x^10
sage: power_mod(11,1,0)
Traceback (most recent call last):
...
ZeroDivisionError: modulus must be nonzero.
The largest prime < n. The result is provably correct. If n <= 1, this function raises a ValueError.
EXAMPLES:
sage: previous_prime(10)
7
sage: previous_prime(7)
5
sage: previous_prime(8)
7
sage: previous_prime(7)
5
sage: previous_prime(5)
3
sage: previous_prime(3)
2
sage: previous_prime(2)
Traceback (most recent call last):
...
ValueError: no previous prime
sage: previous_prime(1)
Traceback (most recent call last):
...
ValueError: no previous prime
sage: previous_prime(-20)
Traceback (most recent call last):
...
ValueError: no previous prime
Return the largest prime power smaller than n.
The result is provably correct. If n is smaller or equal than 2 this function raises an error.
This function simply call the method Integer.previous_prime_power() of Integers.
See also
EXAMPLES:
sage: previous_prime_power(3)
2
sage: previous_prime_power(10)
9
sage: previous_prime_power(7)
5
sage: previous_prime_power(127)
125
The same results can be obtained with:
sage: 3.previous_prime_power()
2
sage: 10.previous_prime_power()
9
sage: 7.previous_prime_power()
5
sage: 127.previous_prime_power()
125
Input less than or equal to \(2\) raises errors:
sage: previous_prime_power(2)
Traceback (most recent call last):
...
ValueError: no prime power less than 2
sage: previous_prime_power(-10)
Traceback (most recent call last):
...
ValueError: no prime power less than 2
sage: n = previous_prime_power(2^16 - 1)
sage: while is_prime(n):
....: n = previous_prime_power(n)
sage: factor(n)
251^2
The prime divisors of n.
INPUT:
OUTPUT:
A list of prime factors of n. For integers, this list is sorted in increasing order.
EXAMPLES:
sage: prime_divisors(1)
[]
sage: prime_divisors(100)
[2, 5]
sage: prime_divisors(2004)
[2, 3, 167]
If n is negative, we do not include -1 among the prime divisors, since -1 is not a prime number:
sage: prime_divisors(-100)
[2, 5]
For polynomials we get all irreducible factors:
sage: R.<x> = PolynomialRing(QQ)
sage: prime_divisors(x^12 - 1)
[x - 1, x + 1, x^2 - x + 1, x^2 + 1, x^2 + x + 1, x^4 - x^2 + 1]
The prime divisors of n.
INPUT:
OUTPUT:
A list of prime factors of n. For integers, this list is sorted in increasing order.
EXAMPLES:
sage: prime_divisors(1)
[]
sage: prime_divisors(100)
[2, 5]
sage: prime_divisors(2004)
[2, 3, 167]
If n is negative, we do not include -1 among the prime divisors, since -1 is not a prime number:
sage: prime_divisors(-100)
[2, 5]
For polynomials we get all irreducible factors:
sage: R.<x> = PolynomialRing(QQ)
sage: prime_divisors(x^12 - 1)
[x - 1, x + 1, x^2 - x + 1, x^2 + 1, x^2 + x + 1, x^4 - x^2 + 1]
List of all positive primes powers between start and stop-1, inclusive. If the second argument is omitted, returns the prime powers up to the first argument.
INPUT:
OUTPUT:
The set of all prime powers between start and stop or, if only one argument is passed, the set of all prime powers between 1 and start. The number \(n\) is a prime power if \(n=p^k\), where \(p\) is a prime number and \(k\) is a positive integer. Thus, \(1\) is not a prime power.
EXAMPLES:
sage: prime_powers(20)
[2, 3, 4, 5, 7, 8, 9, 11, 13, 16, 17, 19]
sage: len(prime_powers(1000))
193
sage: len(prime_range(1000))
168
sage: a = [z for z in range(95,1234) if is_prime_power(z)]
sage: b = prime_powers(95,1234)
sage: len(b)
194
sage: len(a)
194
sage: a[:10]
[97, 101, 103, 107, 109, 113, 121, 125, 127, 128]
sage: b[:10]
[97, 101, 103, 107, 109, 113, 121, 125, 127, 128]
sage: a == b
True
sage: prime_powers(100) == [i for i in range(100) if is_prime_power(i)]
True
sage: prime_powers(10,7)
[]
sage: prime_powers(-5)
[]
sage: prime_powers(-1,3)
[2]
TESTS:
Check that output are always Sage integers (trac ticket #922):
sage: v = prime_powers(10)
sage: type(v[0])
<type 'sage.rings.integer.Integer'>
sage: prime_powers(0,1)
[]
sage: prime_powers(2)
[]
sage: prime_powers(3)
[2]
sage: prime_powers("foo")
Traceback (most recent call last):
...
TypeError: unable to convert 'foo' to an integer
sage: prime_powers(6, "bar")
Traceback (most recent call last):
...
TypeError: unable to convert 'bar' to an integer
Check that long input are accepted (trac ticket #17852):
sage: prime_powers(6l)
[2, 3, 4, 5]
sage: prime_powers(6l,10l)
[7, 8, 9]
Returns the prime-to-m part of n, i.e., the largest divisor of n that is coprime to m.
INPUT:
OUTPUT: Integer
EXAMPLES:
sage: 240.prime_to_m_part(2)
15
sage: 240.prime_to_m_part(3)
80
sage: 240.prime_to_m_part(5)
48
sage: 43434.prime_to_m_part(20)
21717
Returns an iterator over all primes between start and stop-1, inclusive. This is much slower than prime_range, but potentially uses less memory. As with next_prime(), the optional argument proof controls whether the numbers returned are guaranteed to be prime or not.
This command is like the xrange command, except it only iterates over primes. In some cases it is better to use primes than prime_range, because primes does not build a list of all primes in the range in memory all at once. However, it is potentially much slower since it simply calls the next_prime() function repeatedly, and next_prime() is slow.
INPUT:
OUTPUT:
EXAMPLES:
sage: for p in primes(5,10):
... print p
...
5
7
sage: list(primes(13))
[2, 3, 5, 7, 11]
sage: list(primes(10000000000, 10000000100))
[10000000019, 10000000033, 10000000061, 10000000069, 10000000097]
sage: max(primes(10^100, 10^100+10^4, proof=False))
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000009631
sage: next(p for p in primes(10^20, infinity) if is_prime(2*p+1))
100000000000000001243
TESTS:
sage: for a in range(-10, 50):
... for b in range(-10, 50):
... assert list(primes(a,b)) == list(filter(is_prime, xrange(a,b)))
...
sage: sum(primes(-10, 9973, proof=False)) == sum(filter(is_prime, range(-10, 9973)))
True
sage: for p in primes(10, infinity):
... if p > 20: break
... print p
...
11
13
17
19
sage: next(p for p in primes(10,oo)) # checks alternate infinity notation
11
Return the first \(n\) primes.
INPUT:
OUTPUT:
EXAMPLES:
sage: primes_first_n(10)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
sage: len(primes_first_n(1000))
1000
sage: primes_first_n(0)
[]
Return a positive integer that generates the multiplicative group of integers modulo \(n\), if one exists; otherwise, raise a ValueError.
A primitive root exists if \(n=4\) or \(n=p^k\) or \(n=2p^k\), where \(p\) is an odd prime and \(k\) is a nonnegative number.
INPUT:
OUTPUT:
A primitive root of \(n\). If \(n\) is prime, this is the smallest primitive root.
EXAMPLES:
sage: primitive_root(23)
5
sage: primitive_root(-46)
5
sage: primitive_root(25)
2
sage: print [primitive_root(p) for p in primes(100)]
[1, 2, 2, 3, 2, 2, 3, 2, 5, 2, 3, 2, 6, 3, 5, 2, 2, 2, 2, 7, 5, 3, 2, 3, 5]
sage: primitive_root(8)
Traceback (most recent call last):
...
ValueError: no primitive root
Note
It takes extra work to check if \(n\) has a primitive root; to avoid this, use check=False, which may slightly speed things up (but could also result in undefined behavior). For example, the second call below is an order of magnitude faster than the first:
sage: n = 10^50 + 151 # a prime
sage: primitive_root(n)
11
sage: primitive_root(n, check=False)
11
TESTS:
Various special cases:
sage: primitive_root(-1)
0
sage: primitive_root(0)
Traceback (most recent call last):
...
ValueError: no primitive root
sage: primitive_root(1)
0
sage: primitive_root(2)
1
sage: primitive_root(3)
2
sage: primitive_root(4)
3
We test that various numbers without primitive roots give an error - see trac ticket #10836:
sage: primitive_root(15)
Traceback (most recent call last):
...
ValueError: no primitive root
sage: primitive_root(16)
Traceback (most recent call last):
...
ValueError: no primitive root
sage: primitive_root(1729)
Traceback (most recent call last):
...
ValueError: no primitive root
sage: primitive_root(4*7^8)
Traceback (most recent call last):
...
ValueError: no primitive root
Return a sorted list of all squares modulo the integer \(n\) in the range \(0\leq x < |n|\).
EXAMPLES:
sage: quadratic_residues(11)
[0, 1, 3, 4, 5, 9]
sage: quadratic_residues(1)
[0]
sage: quadratic_residues(2)
[0, 1]
sage: quadratic_residues(8)
[0, 1, 4]
sage: quadratic_residues(-10)
[0, 1, 4, 5, 6, 9]
sage: v = quadratic_residues(1000); len(v);
159
Return the product of the prime divisors of n.
This calls n.radical(*args, **kwds). If that doesn’t work, it does n.factor(*args, **kwds) and returns the product of the prime factors in the resulting factorization.
EXAMPLES:
sage: radical(2 * 3^2 * 5^5)
30
sage: radical(0)
Traceback (most recent call last):
...
ArithmeticError: Radical of 0 not defined.
sage: K.<i> = QuadraticField(-1)
sage: radical(K(2))
i + 1
The next example shows how to compute the radical of a number, assuming no prime > 100000 has exponent > 1 in the factorization:
sage: n = 2^1000-1; n / radical(n, limit=100000)
125
Returns a random prime p between \(lbound\) and n (i.e. \(lbound <= p <= n\)). The returned prime is chosen uniformly at random from the set of prime numbers less than or equal to n.
INPUT:
EXAMPLES:
sage: random_prime(100000)
88237
sage: random_prime(2)
2
Here we generate a random prime between 100 and 200:
sage: random_prime(200, lbound=100)
149
If all we care about is finding a pseudo prime, then we can pass in proof=False
sage: random_prime(200, proof=False, lbound=100)
149
TESTS:
sage: type(random_prime(2))
<type 'sage.rings.integer.Integer'>
sage: type(random_prime(100))
<type 'sage.rings.integer.Integer'>
sage: random_prime(1, lbound=-2) #caused Sage hang #10112
Traceback (most recent call last):
...
ValueError: n must be greater than or equal to 2
sage: random_prime(126, lbound=114)
Traceback (most recent call last):
...
ValueError: There are no primes between 114 and 126 (inclusive)
AUTHORS:
This function tries to compute \(x/y\), where \(x/y\) is a rational number in lowest terms such that the reduction of \(x/y\) modulo \(m\) is equal to \(a\) and the absolute values of \(x\) and \(y\) are both \(\le \sqrt{m/2}\). If such \(x/y\) exists, that pair is unique and this function returns it. If no such pair exists, this function raises ZeroDivisionError.
An efficient algorithm for computing rational reconstruction is very similar to the extended Euclidean algorithm. For more details, see Knuth, Vol 2, 3rd ed, pages 656-657.
INPUT:
OUTPUT:
Numerator and denominator \(n\), \(d\) of the unique rational number \(r=n/d\), if it exists, with \(n\) and \(|d| \le \sqrt{N/2}\). Return \((0,0)\) if no such number exists.
The algorithm for rational reconstruction is described (with a complete nontrivial proof) on pages 656-657 of Knuth, Vol 2, 3rd ed. as the solution to exercise 51 on page 379. See in particular the conclusion paragraph right in the middle of page 657, which describes the algorithm thus:
This discussion proves that the problem can be solved efficiently by applying Algorithm 4.5.2X with \(u=m\) and \(v=a\), but with the following replacement for step X2: If \(v3 \le \sqrt{m/2}\), the algorithm terminates. The pair \((x,y)=(|v2|,v3*\mathrm{sign}(v2))\) is then the unique solution, provided that \(x\) and \(y\) are coprime and \(x \le \sqrt{m/2}\); otherwise there is no solution. (Alg 4.5.2X is the extended Euclidean algorithm.)
Knuth remarks that this algorithm is due to Wang, Kornerup, and Gregory from around 1983.
EXAMPLES:
sage: m = 100000
sage: (119*inverse_mod(53,m))%m
11323
sage: rational_reconstruction(11323,m)
119/53
sage: rational_reconstruction(400,1000)
Traceback (most recent call last):
...
ArithmeticError: rational reconstruction of 400 (mod 1000) does not exist
sage: rational_reconstruction(3, 292393)
3
sage: a = Integers(292393)(45/97); a
204977
sage: rational_reconstruction(a, 292393, algorithm='fast')
45/97
sage: rational_reconstruction(293048, 292393)
Traceback (most recent call last):
...
ArithmeticError: rational reconstruction of 655 (mod 292393) does not exist
sage: rational_reconstruction(0, 0)
Traceback (most recent call last):
...
ZeroDivisionError: rational reconstruction with zero modulus
sage: rational_reconstruction(0, 1, algorithm="foobar")
Traceback (most recent call last):
...
ValueError: unknown algorithm 'foobar'
Returns the rising factorial \((x)^a\).
The notation in the literature is a mess: often \((x)^a\), but there are many other notations: GKP: Concrete Mathematics uses \(x^{\overline{a}}\).
The rising factorial is also known as the Pochhammer symbol, see Maple and Mathematica.
Definition: for integer \(a \ge 0\) we have \(x(x+1) \cdots (x+a-1)\). In all other cases we use the GAMMA-function: \(\frac {\Gamma(x+a)} {\Gamma(x)}\).
INPUT:
OUTPUT: the rising factorial
EXAMPLES:
sage: rising_factorial(10,3)
1320
sage: rising_factorial(10,RR('3.0'))
1320.00000000000
sage: rising_factorial(10,RR('3.3'))
2826.38895824964
sage: a = rising_factorial(1+I, I); a
gamma(2*I + 1)/gamma(I + 1)
sage: CC(a)
0.266816390637832 + 0.122783354006372*I
sage: a = rising_factorial(I, 4); a
-10
See falling_factorial(I, 4).
sage: x = polygen(ZZ)
sage: rising_factorial(x, 4)
x^4 + 6*x^3 + 11*x^2 + 6*x
TESTS:
Check that trac ticket #14858 is fixed:
sage: bool(rising_factorial(-4, 2) ==
....: rising_factorial(-4, SR(2)) ==
....: rising_factorial(SR(-4), SR(2)))
True
Check that trac ticket #16770 is fixed:
sage: d = var('d')
sage: type(rising_factorial(d, 0))
<type 'sage.symbolic.expression.Expression'>
AUTHORS:
Given a list of complex numbers (or a list of tuples, where the first element of each tuple is a complex number), we sort the list in a “pretty” order. First come the real numbers (with zero imaginary part), then the complex numbers sorted according to their real part. If two complex numbers have a real part which is sufficiently close, then they are sorted according to their imaginary part.
This is not a useful function mathematically (not least because there’s no principled way to determine whether the real components should be treated as equal or not). It is called by various polynomial root-finders; its purpose is to make doctest printing more reproducible.
We deliberately choose a cumbersome name for this function to discourage use, since it is mathematically meaningless.
EXAMPLES:
sage: import sage.rings.arith
sage: sort_c = sort_complex_numbers_for_display
sage: nums = [CDF(i) for i in range(3)]
sage: for i in range(3):
....: nums.append(CDF(i + RDF.random_element(-3e-11, 3e-11),
....: RDF.random_element()))
....: nums.append(CDF(i + RDF.random_element(-3e-11, 3e-11),
....: RDF.random_element()))
sage: shuffle(nums)
sage: sort_c(nums)
[0.0, 1.0, 2.0, -2.862406201002009e-11 - 0.7088740263015161*I, 2.2108362706985576e-11 - 0.43681052967509904*I, 1.0000000000138833 - 0.7587654737635712*I, 0.9999999999760288 - 0.7238965893336062*I, 1.9999999999874383 - 0.4560801012073723*I, 1.9999999999869107 + 0.6090836283134269*I]
Iterator over the squarefree divisors (up to units) of the element x.
Depends on the output of the prime_divisors function.
INPUT:
EXAMPLES:
sage: list(squarefree_divisors(7))
[1, 7]
sage: list(squarefree_divisors(6))
[1, 2, 3, 6]
sage: list(squarefree_divisors(12))
[1, 2, 3, 6]
TESTS:
Check that the first divisor (i.e. \(1\)) is a Sage integer (see trac ticket #17852):
sage: a = next(squarefree_divisors(14))
sage: a
1
sage: type(a)
<type 'sage.rings.integer.Integer'>
Subfactorial or rencontres numbers, or derangements: number of permutations of \(n\) elements with no fixed points.
INPUT:
OUTPUT:
EXAMPLES:
sage: subfactorial(0)
1
sage: subfactorial(1)
0
sage: subfactorial(8)
14833
AUTHORS:
Write the integer \(n\) as a sum of \(k\) integer squares if possible; otherwise raise a ValueError.
INPUT:
OUTPUT: a tuple \((x_1, ..., x_k)\) of non-negative integers such that their squares sum to \(n\).
EXAMPLES:
sage: sum_of_k_squares(2, 9634)
(15, 97)
sage: sum_of_k_squares(3, 9634)
(0, 15, 97)
sage: sum_of_k_squares(4, 9634)
(1, 2, 5, 98)
sage: sum_of_k_squares(5, 9634)
(0, 1, 2, 5, 98)
sage: sum_of_k_squares(6, 11^1111-1)
(19215400822645944253860920437586326284, 37204645194585992174252915693267578306, 3473654819477394665857484221256136567800161086815834297092488779216863122, 5860191799617673633547572610351797996721850737768032876360978911074629287841061578270832330322236796556721252602860754789786937515870682024273948, 20457423294558182494001919812379023992538802203730791019728543439765347851316366537094696896669915675685581905102118246887673397020172285247862426612188418787649371716686651256443143210952163970564228423098202682066311189439731080552623884051737264415984619097656479060977602722566383385989, 311628095411678159849237738619458396497534696043580912225334269371611836910345930320700816649653412141574887113710604828156159177769285115652741014638785285820578943010943846225597311231847997461959204894255074229895666356909071243390280307709880906261008237873840245959883405303580405277298513108957483306488193844321589356441983980532251051786704380984788999660195252373574924026139168936921591652831237741973242604363696352878914129671292072201700073286987126265965322808664802662993006926302359371379531571194266134916767573373504566621665949840469229781956838744551367172353)
sage: sum_of_k_squares(7, 0)
(0, 0, 0, 0, 0, 0, 0)
sage: sum_of_k_squares(30,999999)
(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 7, 44, 999)
sage: sum_of_k_squares(1, 9)
(3,)
sage: sum_of_k_squares(1, 10)
Traceback (most recent call last):
...
ValueError: 10 is not a sum of 1 square
sage: sum_of_k_squares(1, -10)
Traceback (most recent call last):
...
ValueError: -10 is not a sum of 1 square
sage: sum_of_k_squares(0, 9)
Traceback (most recent call last):
...
ValueError: 9 is not a sum of 0 squares
sage: sum_of_k_squares(0, 0)
()
sage: sum_of_k_squares(7, -1)
Traceback (most recent call last):
...
ValueError: -1 is not a sum of 7 squares
sage: sum_of_k_squares(-1, 0)
Traceback (most recent call last):
...
ValueError: k = -1 must be non-negative
Write the integer \(n\) as a sum of three integer squares if possible; otherwise raise a ValueError.
INPUT:
OUTPUT: a tuple \((a,b,c)\) of non-negative integers such that \(n = a^2 + b^2 + c^2\) with \(a <= b <= c\).
EXAMPLES:
sage: three_squares(389)
(1, 8, 18)
sage: three_squares(946)
(9, 9, 28)
sage: three_squares(2986)
(3, 24, 49)
sage: three_squares(7^100)
(0, 0, 1798465042647412146620280340569649349251249)
sage: three_squares(11^111-1)
(616274160655975340150706442680, 901582938385735143295060746161, 6270382387635744140394001363065311967964099981788593947233)
sage: three_squares(7 * 2^41)
(1048576, 2097152, 3145728)
sage: three_squares(7 * 2^42)
Traceback (most recent call last):
...
ValueError: 30786325577728 is not a sum of 3 squares
sage: three_squares(0)
(0, 0, 0)
sage: three_squares(-1)
Traceback (most recent call last):
...
ValueError: -1 is not a sum of 3 squares
TESTS:
sage: for _ in xrange(100):
....: a = ZZ.random_element(2**16, 2**20)
....: b = ZZ.random_element(2**16, 2**20)
....: c = ZZ.random_element(2**16, 2**20)
....: n = a**2 + b**2 + c**2
....: aa,bb,cc = three_squares(n)
....: assert aa**2 + bb**2 + cc**2 == n
ALGORITHM:
Return the smallest prime divisor <= bound of the positive integer n, or n if there is no such prime. If the optional argument bound is omitted, then bound <= n.
INPUT:
OUTPUT:
EXAMPLES:
sage: trial_division(15)
3
sage: trial_division(91)
7
sage: trial_division(11)
11
sage: trial_division(387833, 300)
387833
sage: # 300 is not big enough to split off a
sage: # factor, but 400 is.
sage: trial_division(387833, 400)
389
Write the integer \(n\) as a sum of two integer squares if possible; otherwise raise a ValueError.
INPUT:
OUTPUT: a tuple \((a,b)\) of non-negative integers such that \(n = a^2 + b^2\) with \(a <= b\).
EXAMPLES:
sage: two_squares(389)
(10, 17)
sage: two_squares(21)
Traceback (most recent call last):
...
ValueError: 21 is not a sum of 2 squares
sage: two_squares(21^2)
(0, 21)
sage: a,b = two_squares(100000000000000000129); a,b
(4418521500, 8970878873)
sage: a^2 + b^2
100000000000000000129
sage: two_squares(2^222+1)
(253801659504708621991421712450521, 2583712713213354898490304645018692)
sage: two_squares(0)
(0, 0)
sage: two_squares(-1)
Traceback (most recent call last):
...
ValueError: -1 is not a sum of 2 squares
TESTS:
sage: for _ in xrange(100):
....: a = ZZ.random_element(2**16, 2**20)
....: b = ZZ.random_element(2**16, 2**20)
....: n = a**2 + b**2
....: aa,bb = two_squares(n)
....: assert aa**2 + bb**2 == n
ALGORITHM:
Return the valuation of m.
This function simply calls the m.valuation() method. See the documentation of m.valuation() for a more precise description.
Note that the use of this functions is discouraged as it is better to use m.valuation() directly.
Note
This is not always a valuation in the mathematical sense. For more information see: sage.rings.finite_rings.integer_mod.IntegerMod_int.valuation
EXAMPLES:
sage: valuation(512,2)
9
sage: valuation(1,2)
0
sage: valuation(5/9, 3)
-2
Valuation of 0 is defined, but valuation with respect to 0 is not:
sage: valuation(0,7)
+Infinity
sage: valuation(3,0)
Traceback (most recent call last):
...
ValueError: You can only compute the valuation with respect to a integer larger than 1.
Here are some other examples:
sage: valuation(100,10)
2
sage: valuation(200,10)
2
sage: valuation(243,3)
5
sage: valuation(243*10007,3)
5
sage: valuation(243*10007,10007)
1
sage: y = QQ['y'].gen()
sage: valuation(y^3, y)
3
sage: x = QQ[['x']].gen()
sage: valuation((x^3-x^2)/(x-4))
2
sage: valuation(4r,2r)
2
sage: valuation(1r,1r)
Traceback (most recent call last):
...
ValueError: You can only compute the valuation with respect to a integer larger than 1.
Return a triple (g,s,t) such that \(g = s\cdot a+t\cdot b = \gcd(a,b)\).
Note
One exception is if \(a\) and \(b\) are not in a principal ideal domain (see Wikipedia article Principal_ideal_domain), e.g., they are both polynomials over the integers. Then this function can’t in general return (g,s,t) as above, since they need not exist. Instead, over the integers, we first multiply \(g\) by a divisor of the resultant of \(a/g\) and \(b/g\), up to sign.
INPUT:
OUTPUT:
Note
There is no guarantee that the returned cofactors (s and t) are minimal.
EXAMPLES:
sage: xgcd(56, 44)
(4, 4, -5)
sage: 4*56 + (-5)*44
4
sage: g, a, b = xgcd(5/1, 7/1); g, a, b
(1, 3, -2)
sage: a*(5/1) + b*(7/1) == g
True
sage: x = polygen(QQ)
sage: xgcd(x^3 - 1, x^2 - 1)
(x - 1, 1, -x)
sage: K.<g> = NumberField(x^2-3)
sage: g.xgcd(g+2)
(1, 1/3*g, 0)
sage: R.<a,b> = K[]
sage: S.<y> = R.fraction_field()[]
sage: xgcd(y^2, a*y+b)
(1, a^2/b^2, ((-a)/b^2)*y + 1/b)
sage: xgcd((b+g)*y^2, (a+g)*y+b)
(1, (a^2 + (2*g)*a + 3)/(b^3 + (g)*b^2), ((-a + (-g))/b^2)*y + 1/b)
Here is an example of a xgcd for two polynomials over the integers, where the linear combination is not the gcd but the gcd multiplied by the resultant:
sage: R.<x> = ZZ[]
sage: gcd(2*x*(x-1), x^2)
x
sage: xgcd(2*x*(x-1), x^2)
(2*x, -1, 2)
sage: (2*(x-1)).resultant(x)
2
This function is similar to the xgcd function, but behaves in a completely different way.
INPUT:
OUTPUT:
This function outputs nothing it just prints something. Note that this function does not feel itself at ease in a html deprived environment.
EXAMPLES:
sage: xkcd(353) # optional - internet
<html><font color='black'><h1>Python</h1><img src="http://imgs.xkcd.com/comics/python.png" title="I wrote 20 short programs in Python yesterday. It was wonderful. Perl, I'm leaving you."><div>Source: <a href="http://xkcd.com/353" target="_blank">http://xkcd.com/353</a></div></font></html>
Extended lcm function: given two positive integers \(m,n\), returns a triple \((l,m_1,n_1)\) such that \(l=\mathop{\mathrm{lcm}}(m,n)=m_1 \cdot n_1\) where \(m_1|m\), \(n_1|n\) and \(\gcd(m_1,n_1)=1\), all with no factorization.
Used to construct an element of order \(l\) from elements of orders \(m,n\) in any group: see sage/groups/generic.py for examples.
EXAMPLES:
sage: xlcm(120,36)
(360, 40, 9)