Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*- """ Miscellaneous arithmetic functions """
#***************************************************************************** # Copyright (C) 2006 William Stein <wstein@gmail.com> # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # http://www.gnu.org/licenses/ #*****************************************************************************
from __future__ import absolute_import, print_function from six.moves import range from six import integer_types
import math
from sage.misc.misc import powerset from sage.misc.misc_c import prod
from sage.libs.pari.all import pari import sage.libs.flint.arith as flint_arith
from sage.structure.element import parent, Element from sage.structure.coerce import py_scalar_to_element
from sage.rings.rational_field import QQ from sage.rings.integer_ring import ZZ from sage.rings.integer import Integer, GCD_list from sage.rings.rational import Rational from sage.rings.real_mpfr import RealNumber from sage.rings.complex_number import ComplexNumber
import sage.rings.fast_arith as fast_arith prime_range = fast_arith.prime_range
from sage.misc.lazy_import import lazy_import lazy_import('sage.arith.all', 'lcm', deprecation=22630) lazy_import('sage.arith.all', 'lcm', '__LCM_sequence', deprecation=22630) lazy_import('sage.arith.all', 'lcm', 'LCM', deprecation=22630)
################################################################## # Elementary Arithmetic ##################################################################
def algdep(z, degree, known_bits=None, use_bits=None, known_digits=None, use_digits=None, height_bound=None, proof=False): """ Returns an irreducible polynomial of degree at most `degree` which is approximately satisfied by the number `z`.
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
- ``height_bound`` - an integer (default: ``None``) specifying the maximum 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: algdep(z, 6) x^2 - x + 1
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 precision::
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
We get an irreducible polynomial even if PARI returns a reducible one::
sage: z = CDF(1, RR(3).sqrt())/2 sage: pari(z).algdep(5) x^5 + x^2 sage: algdep(z, 5) x^2 - x + 1
Check that cases where a constant polynomial might look better get handled correctly::
sage: z=CC(-1)**(1/3) sage: algdep(z,1) x """ raise ValueError("height_bound must be given for proof=True")
if height_bound and abs(z) >= height_bound: return None return x - z
if height_bound and max(abs(z.denominator()), abs(z.numerator())) >= height_bound: return None return z.denominator()*x - z.numerator()
else: #we're supposed to find an irreducible polynomial, so we cannot #return a constant one. If the first LLL basis vector gives #a constant polynomial, use the next one.
# norm on an integer vector invokes Integer.sqrt() which tries to factor... # Given an LLL reduced basis $b_1, ..., b_n$, we only # know that $|b_1| <= 2^((n-1)/2) |x|$ for non-zero $x \in L$.
raise NotImplementedError("proof and height bound only implemented for real and complex numbers")
else:
# f might be reducible. Find the best fitting irreducible factor
algebraic_dependency = algdep
def bernoulli(n, algorithm='default', num_threads=1): r""" Return the n-th Bernoulli number, as a rational number.
INPUT:
- ``n`` - an integer - ``algorithm``:
- ``'default'`` -- use 'flint' for n <= 300000, and 'bernmm' otherwise (this is just a heuristic, and not guaranteed to be optimal on all hardware) - ``'arb'`` -- use the arb library - ``'flint'`` -- use the FLINT library - ``'pari'`` -- use the PARI C library - ``'gap'`` -- use GAP - ``'gp'`` -- use PARI/GP interpreter - ``'magma'`` -- use MAGMA (optional) - ``'bernmm'`` -- use bernmm package (a multimodular algorithm)
- ``num_threads`` - positive integer, number of threads to use (only used for bernmm algorithm)
EXAMPLES::
sage: bernoulli(12) -691/2730 sage: bernoulli(50) 495057205241079648212477525/66
We demonstrate each of the alternative algorithms::
sage: bernoulli(12, algorithm='arb') -691/2730 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 = ['arb','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:
- David Joyner and William Stein """
import sage.interfaces.magma x = sage.interfaces.magma.magma('Bernoulli(%s)'%n) return Rational(x) else: raise ValueError("invalid choice of algorithm")
def factorial(n, algorithm='gmp'): r""" Compute the factorial of `n`, which is the product `1\cdot 2\cdot 3 \cdots (n-1)\cdot n`.
INPUT:
- ``n`` - an integer
- ``algorithm`` - string (default: 'gmp'):
- ``'gmp'`` - use the GMP C-library factorial function
- ``'pari'`` - use PARI's factorial function
OUTPUT: an integer
EXAMPLES::
sage: from sage.arith.misc 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.
- It takes less than a minute to compute the factorial of `10^7` using the GMP algorithm, and the factorial of `10^6` takes less than 4 seconds.
- The GMP algorithm is faster and more memory efficient than the PARI algorithm. E.g., PARI computes `10^7` factorial in 100 seconds on the core duo 2Ghz.
- For comparison, computation in Magma `\leq` 2.12-10 of `n!` is best done using ``*[1..n]``. It takes 113 seconds to compute the factorial of `10^7` and 6 seconds to compute the factorial of `10^6`. Mathematica V5.2 compute the factorial of `10^7` in 136 seconds and the factorial of `10^6` in 7 seconds. (Mathematica is notably very efficient at memory usage when doing factorial calculations.) """ elif algorithm == 'pari': return pari.factorial(n) else: raise ValueError('unknown algorithm')
def is_prime(n): r""" 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 :mod:`arithmetic proof flag <sage.structure.proof.proof>`.
INPUT:
- ``n`` - the object for which to determine primality
.. SEEALSO::
- :meth:`is_pseudoprime` - :meth:`sage.rings.integer.Integer.is_prime`
AUTHORS:
- Kevin Stueve kstueve@uw.edu (2010-01-17): delegated calculation to ``n.is_prime()``
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) """
def is_pseudoprime(n): r""" Test whether ``n`` is a pseudo-prime
The result is *NOT* proven correct - *this is a pseudo-primality test!*.
INPUT:
- ``n`` -- an integer
.. 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 """
def is_prime_power(n, get_data=False): r""" Test whether ``n`` is a positive power of a prime number
This function simply calls the method :meth:`Integer.is_prime_power() <sage.rings.integer.Integer.is_prime_power>` of Integers.
INPUT:
- ``n`` -- an integer
- ``get_data`` -- if set to ``True``, return a pair ``(p,k)`` such that this integer equals ``p^k`` instead of ``True`` or ``(self,0)`` instead of ``False``
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 """
def is_pseudoprime_power(n, get_data=False): r""" 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:
- ``n`` - an integer
- ``get_data`` - (boolean) instead of a boolean return a pair `(p,k)` so that ``n`` equals `p^k` and `p` is a pseudoprime or `(n,0)` otherwise.
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) """
def valuation(m, *args, **kwds): """ 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. """
def prime_powers(start, stop=None): r""" 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:
- ``start`` - an integer. If two inputs are given, a lower bound for the returned set of prime powers. If this is the only input, then it is an upper bound.
- ``stop`` - an integer (default: ``None``). An upper bound for the returned set of prime powers.
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:`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:`17852`)::
sage: prime_powers(6l) [2, 3, 4, 5] sage: prime_powers(6l,10l) [7, 8, 9] """
else:
def primes_first_n(n, leave_pari=False): r""" Return the first `n` primes.
INPUT:
- `n` - a nonnegative integer
OUTPUT:
- a list of the first `n` prime numbers.
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) [] """ raise ValueError("n must be nonnegative")
# # This is from # http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/366178 # It's impressively fast given that it's in Pure Python. # def eratosthenes(n): r""" Return a list of the primes `\leq n`.
This is extremely slow and is for educational purposes only.
INPUT:
- ``n`` - a positive integer
OUTPUT:
- a list of primes less than or equal to n.
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 """
return [] return [ZZ(2)]
def primes(start, stop=None, proof=None): r""" 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 :func:`next_prime`, the optional argument proof controls whether the numbers returned are guaranteed to be prime or not.
This command is like the Python 2 ``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 :func:`next_prime` function repeatedly, and :func:`next_prime` is slow.
INPUT:
- ``start`` - an integer - lower bound for the primes
- ``stop`` - an integer (or infinity) optional argument - giving upper (open) bound for the primes
- ``proof`` - bool or None (default: None) If True, the function yields only proven primes. If False, the function uses a pseudo-primality test, which is much faster for really big numbers but does not provide a proof of primality. If None, uses the global default (see :mod:`sage.structure.proof.proof`)
OUTPUT:
- an iterator over primes from start to stop-1, inclusive
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, range(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 """
else:
def next_prime_power(n): """ 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 :meth:`Integer.next_prime_power() <sage.rings.integer.Integer.next_prime_power>` of Integers.
.. SEEALSO::
- :func:`is_prime_power` (and :meth:`Integer.is_prime_power() <sage.rings.integer.Integer.is_prime_power()>`) - :func:`previous_prime_power` (and :meth:`Integer.previous_prime_power() <sage.rings.integer.Integer.next_prime_power>`)
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 """
def next_probable_prime(n): """ Returns the next probable prime after self, as determined by PARI.
INPUT:
- ``n`` - an integer
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 """
def next_prime(n, proof=None): """ 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 :mod:`sage.structure.proof.proof`)
INPUT:
- ``n`` - integer
- ``proof`` - bool or None (default: None)
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 """
def previous_prime(n): """ 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 """
def previous_prime_power(n): r""" 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 :meth:`Integer.previous_prime_power() <sage.rings.integer.Integer.previous_prime_power>` of Integers.
.. SEEALSO::
- :func:`is_prime_power` (and :meth:`Integer.is_prime_power() <sage.rings.integer.Integer.is_prime_power>`) - :func:`next_prime_power` (and :meth:`Integer.next_prime_power() <sage.rings.integer.Integer.next_prime_power>`)
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 """
def random_prime(n, proof=None, lbound=2): """ 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:
- ``n`` - an integer >= 2.
- ``proof`` - bool or None (default: None) If False, the function uses a pseudo-primality test, which is much faster for really big numbers but does not provide a proof of primality. If None, uses the global default (see :mod:`sage.structure.proof.proof`)
- ``lbound`` - an integer >= 2 lower bound for the chosen primes
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:
- Jon Hanke (2006-08-08): with standard Stein cleanup
- Jonathan Bober (2007-03-17) """ # since we don't want current_randstate to get # pulled when you say "from sage.arith.misc import *". raise ValueError("n must be at least lbound: %s"%(lbound)) # check for Betrand's postulate (proved by Chebyshev) # see J. Nagura, Proc. Japan Acad. 28, (1952). 177-181. # see L. Schoenfeld, Math. Comp. 30 (1976), no. 134, 337-360. else: smallest_prime = ZZ(lbound-1).next_probable_prime()
else: # In order to ensure that the returned prime is chosen # uniformly from the set of primes it is necessary to # choose a random number and then test for primality. # The method of choosing a random number and then returning # the closest prime smaller than it would typically not, # for example, return the first of a pair of twin primes.
def divisors(n): """ Returns a list of all positive integer divisors of the nonzero integer n.
INPUT:
- ``n`` - the element
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] """
class Sigma: """ Return the sum of the k-th powers of the divisors of n.
INPUT:
- ``n`` - integer
- ``k`` - integer (default: 1)
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:
- William Stein: original implementation
- Craig Citro (2007-06-01): rewrote for huge speedup
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 """ def __repr__(self): """ A description of this class, which computes the sum of the k-th powers of the divisors of n.
EXAMPLES::
sage: Sigma().__repr__() 'Function that adds up (k-th powers of) the divisors of n' """
def __call__(self, n, k=1): """ Computes the sum of (the k-th powers of) the divisors of n.
EXAMPLES::
sage: q = Sigma() sage: q(10) 18 sage: q(10,2) 130 """
for p, expt in factor(n)) else: for p,expt in factor(n))
def plot(self, xmin=1, xmax=50, k=1, pointsize=30, rgbcolor=(0,0,1), join=True, **kwds): """ Plot the sigma (sum of k-th powers of divisors) function.
INPUT:
- ``xmin`` - default: 1
- ``xmax`` - default: 50
- ``k`` - default: 1
- ``pointsize`` - default: 30
- ``rgbcolor`` - default: (0,0,1)
- ``join`` - default: True; whether to join the points.
- ``**kwds`` - passed on
EXAMPLES::
sage: p = Sigma().plot() sage: p.ymax() 124.0 """
sigma = Sigma()
def gcd(a, b=None, **kwargs): r""" 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:
- ``a,b`` - two elements of a ring with gcd or
- ``a`` - a list or tuple of elements of a ring with gcd
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 :trac:`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 ... 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:`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:`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
""" # Most common use case first:
GCD = gcd
def __GCD_sequence(v, **kwargs): """ Internal function returning the gcd of the elements of a sequence
INPUT:
- ``v`` - A sequence (possibly empty)
OUTPUT: The gcd of the elements of the sequence as an element of the sequence's universe, or the integer 0 if the sequence is empty.
EXAMPLES::
sage: from sage.arith.misc import __GCD_sequence sage: from sage.structure.sequence import Sequence sage: l = () sage: __GCD_sequence(l) 0 sage: __GCD_sequence(Sequence(srange(10))) 1 sage: X=polygen(QQ) sage: __GCD_sequence(Sequence((2*X+4,2*X^2,2))) 1 sage: X=polygen(ZZ) sage: __GCD_sequence(Sequence((2*X+4,2*X^2,2))) 2 """ else: g = ZZ(0)
def xlcm(m, n): r""" 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) """ # higher power than m
def xgcd(a, b): r""" 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:`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:
- ``a, b`` - integers or more generally, element of a ring for which the xgcd make sense (e.g. a field or univariate polynomials).
OUTPUT:
- ``g, s, t`` - such that `g = s\cdot a + t\cdot b`
.. 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 """
XGCD = xgcd
## def XGCD_python(a, b): ## """ ## Returns triple (g,p,q) such that g = p*a+b*q = GCD(a,b). ## This function should behave exactly the same as XGCD, ## but is implemented in pure python. ## """ ## if a == 0 and b == 0: ## return (0,0,1) ## if a == 0: ## return (abs(b), 0, b/abs(b)) ## if b == 0: ## return (abs(a), a/abs(a), 0) ## psign = 1 ## qsign = 1 ## if a < 0: ## a = -a ## psign = -1 ## if b < 0: ## b = -b ## qsign = -1 ## p = 1; q = 0; r = 0; s = 1 ## while b != 0: ## c = a % b ## quot = a/b ## a = b; b = c ## new_r = p - quot*r ## new_s = q - quot*s ## p = r; q = s ## r = new_r; s = new_s ## return (a, p*psign, q*qsign)
def xkcd(n=""): r""" This function is similar to the xgcd function, but behaves in a completely different way.
INPUT:
- ``n`` -- an integer (optional)
OUTPUT: a fragment of HTML
EXAMPLES::
sage: xkcd(353) # optional - internet <h1>Python</h1><img src="https://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> """ import contextlib import json from sage.misc.html import html
# import compatible with py2 and py3 from six.moves.urllib.request import urlopen from six.moves.urllib.error import HTTPError, URLError
data = None url = "http://dynamic.xkcd.com/api-0/jsonp/comic/{}".format(n)
try: with contextlib.closing(urlopen(url)) as f: data = f.read() except HTTPError as error: if error.getcode() == 400: # this error occurs when asking for a non valid comic number raise RuntimeError("Could not obtain comic data from {}. Maybe you should enable time travel!".format(url)) except URLError: pass
if n == 1024: data = None
if data: data = json.loads(data) img = data['img'] alt = data['alt'] title = data['safe_title'] link = "http://xkcd.com/{}".format(data['num']) return html('<h1>{}</h1><img src="{}" title="{}">'.format(title, img, alt) + '<div>Source: <a href="{0}" target="_blank">{0}</a></div>'.format(link))
# TODO: raise this error in such a way that it's not clear that # it is produced by sage, see http://xkcd.com/1024/ return html('<script> alert("Error: -41"); </script>')
def inverse_mod(a, m): """ 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 """
####################################################### # Functions to find the fastest available commands # for gcd and inverse_mod #######################################################
def get_gcd(order): """ Return the fastest gcd function for integers of size no larger than order.
EXAMPLES::
sage: sage.arith.misc.get_gcd(4000) <built-in method gcd_int of sage.rings.fast_arith.arith_int object at ...> sage: sage.arith.misc.get_gcd(400000) <built-in method gcd_longlong of sage.rings.fast_arith.arith_llong object at ...> sage: sage.arith.misc.get_gcd(4000000000) <function gcd at ...> """ else:
def get_inverse_mod(order): """ Return the fastest inverse_mod function for integers of size no larger than order.
EXAMPLES::
sage: sage.arith.misc.get_inverse_mod(6000) <built-in method inverse_mod_int of sage.rings.fast_arith.arith_int object at ...> sage: sage.arith.misc.get_inverse_mod(600000) <built-in method inverse_mod_longlong of sage.rings.fast_arith.arith_llong object at ...> sage: sage.arith.misc.get_inverse_mod(6000000000) <function inverse_mod at ...> """ else:
# def sqrt_mod(a, m): # """A square root of a modulo m."""
# def xxx_inverse_mod(a, m): # """The inverse of a modulo m.""" # g,s,t = XGCD(a,m) # if g != 1: # raise "inverse_mod(a=%s,m=%s), error since GCD=%s"%(a,m,g) # return s
def power_mod(a,n,m): """ 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. """ return 0 return 1
def rational_reconstruction(a, m, algorithm='fast'): r""" 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:
- ``a`` -- an integer
- ``m`` -- a modulus
- ``algorithm`` -- (default: 'fast')
- ``'fast'`` - a fast implementation using direct MPIR calls in Cython.
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' """ else:
def mqrr_rational_reconstruction(u, m, T): r""" Maximal Quotient Rational Reconstruction.
For research purposes only - this is pure Python, so slow.
INPUT:
- ``u, m, T`` - integers such that `m > u \ge 0`, `T > 0`.
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) """ if m > T: return (0,1) else: return None return None
######################
def trial_division(n, bound=None): """ 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:
- ``n`` - a positive integer
- ``bound`` - (optional) a positive integer
OUTPUT:
- ``int`` - a prime p=bound that divides n, or n if there is no such prime.
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 """ else:
def factor(n, proof=None, int_=False, algorithm='pari', verbose=0, **kwds): """ 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.
EXAMPLES::
sage: f(n)=n^2 sage: is_prime(f(3)) False sage: factor(f(3)) 9
INPUT:
- ``n`` - an nonzero integer
- ``proof`` - bool or None (default: None)
- ``int_`` - bool (default: False) whether to return answers as Python ints
- ``algorithm`` - string
- ``'pari'`` - (default) use the PARI c library
- ``'kash'`` - use KASH computer algebra system (requires the optional kash package be installed)
- ``'magma'`` - use Magma (requires magma be installed)
- ``verbose`` - integer (default: 0); PARI's debug variable is set to this; e.g., set to 4 or 8 to see lots of output during factorization.
OUTPUT:
- factorization of n
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 :class:`~sage.structure.factorization.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: factorization of 0 is 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 :mod:`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]
We can factor Python and numpy numbers::
sage: factor(math.pi) 3.141592653589793 sage: import numpy sage: factor(numpy.int8(30)) 2 * 3 * 5
TESTS::
sage: factor(Mod(4, 100)) Traceback (most recent call last): ... TypeError: unable to factor 4 sage: factor("xyz") Traceback (most recent call last): ... TypeError: unable to factor 'xyz' """ # Maybe n is not a Sage Element, try to convert it # Either n was a Sage Element without a factor() method # or we cannot it convert it to Sage
verbose=verbose, **kwds)
# Polynomial or other factorable object # Maybe the factor() method doesn't have a proof option
def radical(n, *args, **kwds): """ 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 """
def prime_divisors(n): """ The prime divisors of ``n``.
INPUT:
- ``n`` -- any object which can be factored
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] """
prime_factors = prime_divisors
def odd_part(n): r""" 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 """ n = ZZ(n)
def prime_to_m_part(n,m): """ Returns the prime-to-m part of n, i.e., the largest divisor of n that is coprime to m.
INPUT:
- ``n`` - Integer (nonzero)
- ``m`` - Integer
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 """
def is_square(n, root=False): """ 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:
- ``n`` - an integer
- ``root`` - whether or not to also return a square root (default: False)
OUTPUT:
- ``bool`` - whether or not a square
- ``object`` - (optional) an actual square if found, and None otherwise.
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) """ else: return False, None
def is_squarefree(n): """ 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 """
return False
################################################################# # Euler phi function ################################################################# class Euler_Phi: r""" 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:
- ``n`` - an integer
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:
- William Stein
- Alex Clemesha (2006-01-10): some examples """ def __repr__(self): """ Returns a string describing this class.
EXAMPLES::
sage: Euler_Phi().__repr__() 'Number of positive integers <=n but relatively prime to n' """
def __call__(self, n): """ Calls the euler_phi function.
EXAMPLES::
sage: Euler_Phi()(10) 4 sage: Euler_Phi()(720) 192 """
def plot(self, xmin=1, xmax=50, pointsize=30, rgbcolor=(0,0,1), join=True, **kwds): """ Plot the Euler phi function.
INPUT:
- ``xmin`` - default: 1
- ``xmax`` - default: 50
- ``pointsize`` - default: 30
- ``rgbcolor`` - default: (0,0,1)
- ``join`` - default: True; whether to join the points.
- ``**kwds`` - passed on
EXAMPLES::
sage: p = Euler_Phi().plot() sage: p.ymax() 46.0 """
euler_phi = Euler_Phi()
def crt(a,b,m=None,n=None): r""" Returns a solution to a Chinese Remainder Theorem problem.
INPUT:
- ``a``, ``b`` - two residues (elements of some ring for which extended gcd is available), or two lists, one of residues and one of moduli.
- ``m``, ``n`` - (default: ``None``) two moduli, or ``None``.
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.
.. SEEALSO::
- :func:`CRT_list`
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 """
CRT = crt
def CRT_list(v, moduli): r""" 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.
.. SEEALSO::
- :func:`crt`
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
"""
def CRT_basis(moduli): r""" Returns a CRT basis for the given moduli.
INPUT:
- ``moduli`` - list of pairwise coprime moduli `m` which admit an extended Euclidean algorithm
OUTPUT:
- a list of elements `a_i` of the same length as `m` such that `a_i` is congruent to 1 modulo `m_i` and to 0 modulo `m_j` for `j\not=i`.
.. 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]] """ return []
def CRT_vectors(X, moduli): r""" 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 :func:`CRT` to each entry.
INPUT:
- ``X`` - list or tuple, consisting of lists/tuples/vectors/etc of integers of the same length - ``moduli`` - list of len(X) moduli
OUTPUT:
- ``list`` - application of CRT componentwise.
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] """ # First find the CRT basis: return [] raise ValueError("number of moduli must equal length of X")
def binomial(x, m, **kwds): r""" Return the binomial coefficient
.. MATH::
\binom{x}{m} = x (x-1) \cdots (x-m+1) / m!
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
.. MATH::
\binom{x}{m} = \binom{x}{x-m}
If `m < 0`, return `0`.
INPUT:
- ``x``, ``m`` - numbers or symbolic expressions. Either ``m`` or ``x-m`` must be an integer.
OUTPUT: number or symbolic expression (if input is symbolic)
EXAMPLES::
sage: from sage.arith.misc 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:`3309`::
sage: a = binomial(RR(1140000.78), 23310000)
We test conversion of arguments to Integers -- see :trac:`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:`7562`, :trac:`9633`, and :trac:`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) <... 'float'> sage: binomial(float(1000), 1001) 0.0
Test more output types::
sage: type(binomial(5r, 2)) <... 'int'> sage: type(binomial(5, 2r)) <type 'sage.rings.integer.Integer'>
sage: type(binomial(5.0r, 2)) <... '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) 6*y^3 + 3*y^2 + 5*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 :func:`~sage.functions.other.binomial` from the module :mod:`sage.functions.other`::
sage: from sage.functions.other import binomial sage: binomial(k, i) binomial(k, i) """
# case 1: native binomial implemented on x
# case 2: conversion to integers else: # Check invertibility of factorial(m) in P # Assume that P has characteristic zero (can be int, float, ...) else:
# case 3: rational, real numbers, complex numbers -> use pari
# case 4: naive method return P(0)
def multinomial(*ks): r""" Return the multinomial coefficient
INPUT:
- An arbitrary number of integer arguments `k_1,\dots,k_n` - A list of integers `[k_1,\dots,k_n]`
OUTPUT:
Returns the integer:
.. MATH::
\binom{k_1 + \cdots + k_n}{k_1, \cdots, k_n} =\frac{\left(\sum_{i=1}^n k_i\right)!}{\prod_{i=1}^n k_i!} = \prod_{i=1}^n \binom{\sum_{j=1}^i k_j}{k_i}
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:
- Gabriel Ebner """ raise ValueError("multinomial takes only one list argument")
def binomial_coefficients(n): r""" 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:
- ``n`` - an integer
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:
- Fredrik Johansson """
def multinomial_coefficients(m, n): r""" 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:
- ``m`` - integer - ``n`` - integer
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) {}
""" else: else: j = m # enumerate tuples in co-lex order # compute next tuple else: # compute the value # NB: the initialization of v was done above
def kronecker_symbol(x,y): """ The Kronecker symbol `(x|y)`.
INPUT:
- ``x`` -- integer
- ``y`` -- integer
OUTPUT:
- an integer
EXAMPLES::
sage: kronecker_symbol(13,21) -1 sage: kronecker_symbol(101,4) 1
This is also available as :func:`kronecker`::
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 """
kronecker = kronecker_symbol
def legendre_symbol(x,p): r""" The Legendre symbol `(x|p)`, for `p` prime.
.. note::
The :func:`kronecker_symbol` command extends the Legendre symbol to composite moduli and `p=2`.
INPUT:
- ``x`` - integer
- ``p`` - an odd prime number
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 """
def jacobi_symbol(a,b): r""" The Jacobi symbol of integers a and b, where b is odd.
.. note::
The :func:`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:
- ``a`` - an integer
- ``b`` - an odd integer
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 """
def primitive_root(n, check=True): """ 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:
- ``n`` -- a non-zero integer - ``check`` -- bool (default: True); if False, then `n` is assumed to be a positive integer possessing a primitive root, and behavior is undefined otherwise.
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:`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 """ # n-1 is a primitive root for n in {1,2,3,4} else: # n even
def nth_prime(n): """
Return the n-th prime number (1-indexed, so that 2 is the 1st prime.)
INPUT:
- ``n`` -- a positive integer
OUTPUT:
- the n-th prime number
EXAMPLES::
sage: nth_prime(3) 5 sage: nth_prime(10) 29 sage: nth_prime(10^7) 179424673
::
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 """
def quadratic_residues(n): r""" 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 """
class Moebius: r""" Returns the value of the Möbius 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:
- ``n`` - anything that can be factored.
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 """ def __call__(self, n): """ EXAMPLES::
sage: Moebius().__call__(7) -1 """ # Use a generic algorithm. n = -n
# Use fast PARI algorithm
def __repr__(self): """ Returns a description of this function.
EXAMPLES::
sage: q = Moebius() sage: q.__repr__() 'The Moebius function' """
def plot(self, xmin=0, xmax=50, pointsize=30, rgbcolor=(0,0,1), join=True, **kwds): """ Plot the Möbius function.
INPUT:
- ``xmin`` - default: 0
- ``xmax`` - default: 50
- ``pointsize`` - default: 30
- ``rgbcolor`` - default: (0,0,1)
- ``join`` - default: True; whether to join the points (very helpful in seeing their order).
- ``**kwds`` - passed on
EXAMPLES::
sage: p = Moebius().plot() sage: p.ymax() 1.0 """
def range(self, start, stop=None, step=None): """ Return the Möbius function evaluated at the given range of values, i.e., the image of the list range(start, stop, step) under the Möbius 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 """ start, stop = 1, int(start) else: else:
self.range(step, stop, step)
stop-start, start)) else: n, step, start))
moebius = Moebius()
## Note: farey, convergent, continued_fraction_list and convergents have been moved to ## sage.rings.continued_fraction
def continuant(v, n=None): r""" 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
- `K_0() = 1` - `K_1(x_1) = x_1` - `K_n(x_1, \cdots, x_n) = K_{n-1}(x_n, \cdots x_{n-1})x_n + K_{n-2}(x_1, \cdots, x_{n-2})`
If ``n = None`` or ``n > len(v)`` the default ``n = len(v)`` is used.
INPUT:
- ``v`` - list or tuple of elements of a ring - ``n`` - optional integer
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
.. MATH::
K_n(z,z,\cdots,z) = \sum_{k=0}^n \binom{n-k}{k} z^{n-2k}
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:
- Jaap Spies (2007-02-06) """ return 1 return v[0]
def number_of_divisors(n): """ Return the number of divisors of the integer n.
INPUT:
- ``n`` - a nonzero integer
OUTPUT:
- an integer, the number of divisors of n
EXAMPLES::
sage: number_of_divisors(100) 9 sage: number_of_divisors(-720) 30 """ raise ValueError("input must be nonzero")
def hilbert_symbol(a, b, p, algorithm="pari"): """ 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:
- ``a, b`` - integers
- ``p`` - integer; either prime or -1 (which represents the archimedean place)
- ``algorithm`` - string
- ``'pari'`` - (default) use the PARI C library
- ``'direct'`` - use a Python implementation
- ``'all'`` - use both PARI and direct and check that the results agree, then return the common answer
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:
- William Stein and David Kohel (2006-01-05) """ raise ValueError("p must be prime or -1")
elif p == 2 and (b%4) == 3: if kronecker(a+b,p) == -1: return -one elif kronecker(b,p) == -1: return -one if p == 2 and (a%4) == 3: if kronecker(a+b,p) == -1: return -one elif kronecker(a,p) == -1: return -one raise RuntimeError("There is a bug in hilbert_symbol; two ways of computing the Hilbert symbol (%s,%s)_%s disagree"%(a,b,p)) else: raise ValueError("Algorithm %s not defined"%algorithm)
def hilbert_conductor(a, b): """ 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:
- ``a``, ``b`` -- integers
OUTPUT:
- squarefree positive integer
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:
- Gonzalo Tornaria (2009-03-02) """
def hilbert_conductor_inverse(d): """ 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:
- ``d`` -- square-free positive integer
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:
- Gonzalo Tornaria (2009-03-02)
TESTS::
sage: for i in range(100): ....: d = ZZ.random_element(2**32).squarefree_part() ....: if hilbert_conductor(*hilbert_conductor_inverse(d)) != d: ....: print("hilbert_conductor_inverse failed for d = {}".format(d)) """ q = next_prime(q) else: else:
############################################################################## ## falling and rising factorials ## By Jaap Spies ## ## Copyright (C) 2006 Jaap Spies <j.spies@hccnet.nl> ## Copyright (C) 2006 William Stein <wstein@gmail.com> ## ## Distributed under the terms of the GNU General Public License (GPL) ## http://www.gnu.org/licenses/ ##############################################################################
def falling_factorial(x, a): r""" 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:
- ``x`` - element of a ring
- ``a`` - a non-negative integer or
OR
- ``x and a`` - any numbers
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:`14858` is fixed::
sage: falling_factorial(-4, SR(2)) 20
Check that :trac:`16770` is fixed::
sage: d = var('d') sage: parent(falling_factorial(d, 0)) Symbolic Ring
Check that :trac:`20075` is fixed::
sage: bool(falling_factorial(int(4), int(2)) == falling_factorial(4,2)) True
AUTHORS:
- Jaap Spies (2006-03-05) """ (isinstance(a, Expression) and a.is_integer())) and a >= 0:
def rising_factorial(x, a): r""" 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:
- ``x`` - element of a ring
- ``a`` - a non-negative integer or
- ``x and a`` - any numbers
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:`14858` is fixed::
sage: bool(rising_factorial(-4, 2) == ....: rising_factorial(-4, SR(2)) == ....: rising_factorial(SR(-4), SR(2))) True
Check that :trac:`16770` is fixed::
sage: d = var('d') sage: parent(rising_factorial(d, 0)) Symbolic Ring
Check that :trac:`20075` is fixed::
sage: bool(rising_factorial(int(4), int(2)) == rising_factorial(4,2)) True
AUTHORS:
- Jaap Spies (2006-03-05) """ (isinstance(a, Expression) and a.is_integer())) and a >= 0:
def integer_ceil(x): """ 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 """
def integer_floor(x): r""" Return the largest integer `\leq x`.
INPUT:
- ``x`` - an object that has a floor method or is coercible to int
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 """
def two_squares(n): """ Write the integer `n` as a sum of two integer squares if possible; otherwise raise a ``ValueError``.
INPUT:
- ``n`` -- an integer
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 range(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:
See http://www.schorn.ch/howto.html """
# Start by factoring n (which seems to be unavoidable)
# First check whether it is possible to write n as a sum of two # squares: all prime powers p^e must have p = 2 or p = 1 mod 4 # or e even. raise ValueError("%s is not a sum of 2 squares"%n)
# We run over all factors of n, write each factor p^e as # a sum of 2 squares and accumulate the product # (using multiplication in Z[I]) in a^2 + b^2. # (a + bi) *= (1 + I) else: # p = 1 mod 4 # Find a square root of -1 mod p. # If y is a non-square, then y^((p-1)/4) is a square root of -1. # Apply Cornacchia's algorithm to write p as r^2 + s^2.
# Multiply (a + bI) by (r + sI)
else:
def three_squares(n): """ Write the integer `n` as a sum of three integer squares if possible; otherwise raise a ``ValueError``.
INPUT:
- ``n`` -- an integer
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 range(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:
See http://www.schorn.ch/howto.html """
# First, remove all factors 4 from n
# Let x be the largest integer at most sqrt(N) # We need to check for this special case, # otherwise N - x^2 will always factor.
# Consider different cases to find an x such that N - x^2 is easily # written as the sum of 2 squares, because it is either p or 2p, # with p a prime which is 1 mod 4. # Write N = x^2 + p with x even, p = 1 mod 4 prime # Write N = x^2 + p with x odd, p = 1 mod 4 prime # Write N = x^2 + 2p with x odd, p = 1 mod 4 prime else: # 7 mod 8
# We found no good x, brute force instead. # Normally, this should only happen for small values of N. if N > 10000: from warnings import warn warn("Brute forcing sum of 3 squares for large N = %s"%N, RuntimeWarning) x = N.isqrt()
# In the usual case, this loop will only be executed once, since # we already know the "right" value of x. # This will only really loop if we hit the "x < 0" case above. except ValueError: x -= 1 assert x >= 0
elif x >= a: return (a*m, x*m, b*m) else: return (x*m, a*m, b*m)
def four_squares(n): """ Write the integer `n` as a sum of four integer squares.
INPUT:
- ``n`` -- an integer
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 range(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 """
raise ValueError("%s is not a sum of 4 squares"%n)
# First, remove all factors 4 from n
# Subtract a suitable x^2 such that N - x^2 is 1,2,3,5,6 mod 8, # which can then be written as a sum of 3 squares.
# Correct sorting is guaranteed by construction
def sum_of_k_squares(k,n): """ Write the integer `n` as a sum of `k` integer squares if possible; otherwise raise a ``ValueError``.
INPUT:
- ``k`` -- a non-negative integer
- ``n`` -- an integer
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 """
# Recursively subtract the largest square
def subfactorial(n): r""" Subfactorial or rencontres numbers, or derangements: number of permutations of `n` elements with no fixed points.
INPUT:
- ``n`` - non negative integer
OUTPUT:
- ``integer`` - function value
EXAMPLES::
sage: subfactorial(0) 1 sage: subfactorial(1) 0 sage: subfactorial(8) 14833
AUTHORS:
- Jaap Spies (2007-01-23) """
def is_power_of_two(n): r""" This function returns True if and only if ``n`` is a power of 2
INPUT:
- ``n`` - integer
OUTPUT:
- ``True`` - if n is a power of 2
- ``False`` - if not
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 """
def differences(lis, n=1): """ 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:
- Timothy Clemans (2008-03-09) """ raise ValueError('n must be greater than 0')
def _key_complex_for_display(a): """ Key function to sort complex numbers for display only.
Real numbers (with a zero imaginary part) come before complex numbers, and are sorted. Complex numbers are sorted by their real part unless their real parts are quite close, in which case they are sorted by their imaginary part.
EXAMPLES::
sage: import sage.arith.misc sage: key_c = sage.arith.misc._key_complex_for_display
sage: key_c(CC(5)) (0, 5.00000000000000) sage: key_c(CC(5, 5)) (1, 5.00000000, 5.00000000000000)
sage: CIF200 = ComplexIntervalField(200) sage: key_c(CIF200(5)) (0, 5) sage: key_c(CIF200(5, 5)) (1, 5.00000000, 5) """ else:
def sort_complex_numbers_for_display(nums): r""" 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 the same real part, then they are sorted according to their imaginary part.
This is not a useful function mathematically (not least because there is 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.arith.misc 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] """
else:
def fundamental_discriminant(D): r""" 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:
- ``D`` - an integer
OUTPUT:
- an integer, the fundamental discriminant
EXAMPLES::
sage: fundamental_discriminant(102) 408 sage: fundamental_discriminant(720) 5 sage: fundamental_discriminant(2) 8
"""
def squarefree_divisors(x): """ Iterator over the squarefree divisors (up to units) of the element x.
Depends on the output of the prime_divisors function.
INPUT:
- x -- an element of any ring for which the prime_divisors function works.
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:`17852`)::
sage: a = next(squarefree_divisors(14)) sage: a 1 sage: type(a) <type 'sage.rings.integer.Integer'> """
def dedekind_sum(p, q, algorithm='default'): r""" Return the Dedekind sum `s(p,q)` defined for integers `p`, `q` as
.. MATH::
s(p,q) = \sum_{i=0}^{q-1} \left(\!\left(\frac{i}{q}\right)\!\right) \left(\!\left(\frac{pi}{q}\right)\!\right)
where
.. MATH::
((x))=\begin{cases} x-\lfloor x \rfloor - \frac{1}{2} &\mbox{if } x \in \QQ \setminus \ZZ \\ 0 & \mbox{if } x \in \ZZ. \end{cases}
.. 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:
- ``p``, ``q`` -- integers - ``algorithm`` -- must be one of the following
- ``'default'`` - (default) use FLINT - ``'flint'`` - use FLINT - ``'pari'`` - use PARI (gives different results if `p` and `q` are not coprime)
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:
- [Ap1997]_
- :wikipedia:`Dedekind\_sum` """
raise ValueError('unknown algorithm')
def gauss_sum(char_value, finite_field): r""" Return the Gauss sums for a general finite field.
INPUT:
- ``char_value`` -- choice of multiplicative character, given by its value on the ``finite_field.multiplicative_generator()``
- ``finite_field`` -- a finite field
OUTPUT:
an element of the parent ring of ``char_value``, that can be any field containing enough roots of unity, for example the ``UniversalCyclotomicField``, ``QQbar`` or ``ComplexField``
For a finite field `F` of characteristic `p`, the Gauss sum associated to a multiplicative character `\chi` (with values in a ring `K`) is defined as
.. MATH::
\sum_{x \in F^{\times}} \chi(x) \zeta_p^{\operatorname{Tr} x},
where `\zeta_p \in K` is a primitive root of unity of order `p` and Tr is the trace map from `F` to its prime field `\GF{p}`.
For more info on Gauss sums, see :wikipedia:`Gauss_sum`.
.. TODO::
Implement general Gauss sums for an arbitrary pair ``(multiplicative_character, additive_character)``
EXAMPLES::
sage: from sage.arith.misc import gauss_sum sage: F = GF(5); q = 5 sage: zq = UniversalCyclotomicField().zeta(q-1) sage: L = [gauss_sum(zq**i,F) for i in range(5)]; L [-1, E(20)^4 + E(20)^13 - E(20)^16 - E(20)^17, E(5) - E(5)^2 - E(5)^3 + E(5)^4, E(20)^4 - E(20)^13 - E(20)^16 + E(20)^17, -1] sage: [g*g.conjugate() for g in L] [1, 5, 5, 5, 1]
sage: F = GF(11**2); q = 11**2 sage: zq = UniversalCyclotomicField().zeta(q-1) sage: g = gauss_sum(zq**4,F) sage: g*g.conjugate() 121
TESTS::
sage: F = GF(11); q = 11 sage: zq = UniversalCyclotomicField().zeta(q-1) sage: gauss_sum(zq**2,F).n(60) 2.6361055643248352 + 2.0126965627574471*I
sage: zq = QQbar.zeta(q-1) sage: gauss_sum(zq**2,F) 2.636105564324836? + 2.012696562757447?*I
sage: zq = ComplexField(60).zeta(q-1) sage: gauss_sum(zq**2,F) 2.6361055643248352 + 2.0126965627574471*I
sage: F = GF(7); q = 7 sage: zq = QQbar.zeta(q-1) sage: D = DirichletGroup(7, QQbar) sage: all(D[i].gauss_sum()==gauss_sum(zq**i,F) for i in range(6)) True
sage: gauss_sum(1,QQ) Traceback (most recent call last): ... ValueError: second input must be a finite field
.. SEEALSO::
- :func:`sage.rings.padics.misc.gauss_sum` for a `p`-adic version - :meth:`sage.modular.dirichlet.DirichletCharacter.gauss_sum` for prime finite fields - :meth:`sage.modular.dirichlet.DirichletCharacter.gauss_sum_numerical` for prime finite fields """
|