Skip to content
This repository has been archived by the owner on Jan 30, 2023. It is now read-only.

Commit

Permalink
Merge branch 't/20205/get_rid_of_factorint_withproof_sage_in_pari_int…
Browse files Browse the repository at this point in the history
…erface' into HEAD
  • Loading branch information
jdemeyer committed Mar 15, 2016
2 parents 4bb8337 + c0ed97a commit edc5ce2
Show file tree
Hide file tree
Showing 6 changed files with 57 additions and 137 deletions.
3 changes: 3 additions & 0 deletions src/sage/libs/pari/decl.pxi
Original file line number Diff line number Diff line change
@@ -1 +1,4 @@
from sage.misc.superseded import deprecation
deprecation(20205, '''pari/decl.pxi is deprecated, use "from sage.libs.pari.paridecl cimport *" instead''')

from sage.libs.pari.paridecl cimport *
106 changes: 39 additions & 67 deletions src/sage/libs/pari/gen.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,6 @@ include "cysignals/signals.pxi"

cimport cython

cdef extern from "misc.h":
int factorint_withproof_sage(GEN* ans, GEN x, GEN cutoff)

from sage.libs.gmp.mpz cimport *
from sage.libs.gmp.pylong cimport mpz_set_pylong
from sage.libs.pari.closure cimport objtoclosure
Expand Down Expand Up @@ -5578,10 +5575,10 @@ cdef class gen(gen_auto):
0
sage: pari(500509).primepi()
41581
sage: pari(10^7).primepi()
664579
"""
pari_catch_sig_on()
if self > P._primelimit():
P.init_primes(self + 10)
if signe(self.g) != 1:
pari_catch_sig_off()
return P.PARI_ZERO
Expand Down Expand Up @@ -9045,35 +9042,30 @@ cdef class gen(gen_auto):
pari_catch_sig_on()
return P.new_gen(matfrobenius(self.g, flag, 0))


###########################################
# polarit2.c
###########################################
def factor(gen self, limit=-1, bint proof=1):
def factor(self, long limit=-1, proof=None):
"""
Return the factorization of x.
INPUT:
- ``limit`` -- (default: -1) is optional and can be set
whenever x is of (possibly recursive) rational type. If limit is
set return partial factorization, using primes up to limit (up to
primelimit if limit=0).
set, return partial factorization, using primes up to limit.
- ``proof`` -- (default: True) optional. If False (not the default),
- ``proof`` -- optional flag. If ``False`` (not the default),
returned factors larger than `2^{64}` may only be pseudoprimes.
.. note::
In the standard PARI/GP interpreter and C-library the
factor command *always* has proof=False, so beware!
If ``True``, always check primality. If not given, use the
global PARI default ``factor_proven`` which is ``True`` by
default in Sage.
EXAMPLES::
sage: pari('x^10-1').factor()
[x - 1, 1; x + 1, 1; x^4 - x^3 + x^2 - x + 1, 1; x^4 + x^3 + x^2 + x + 1, 1]
sage: pari(2^100-1).factor()
[3, 1; 5, 3; 11, 1; 31, 1; 41, 1; 101, 1; 251, 1; 601, 1; 1801, 1; 4051, 1; 8101, 1; 268501, 1]
sage: pari(2^100-1).factor(proof=True)
[3, 1; 5, 3; 11, 1; 31, 1; 41, 1; 101, 1; 251, 1; 601, 1; 1801, 1; 4051, 1; 8101, 1; 268501, 1]
sage: pari(2^100-1).factor(proof=False)
[3, 1; 5, 3; 11, 1; 31, 1; 41, 1; 101, 1; 251, 1; 601, 1; 1801, 1; 4051, 1; 8101, 1; 268501, 1]
Expand All @@ -9082,30 +9074,45 @@ cdef class gen(gen_auto):
sage: pari(next_prime(10^50)*next_prime(10^60)*next_prime(10^4)).factor(10^5)
[10007, 1; 100000000000000000000000000000000000000000000000151000000000700000000000000000000000000000000000000000000001057, 1]
Setting a limit is invalid when factoring polynomials::
sage: pari('x^11 + 1').factor(limit=17)
Traceback (most recent call last):
...
PariError: incorrect type in boundfact (t_POL)
PARI doesn't have an algorithm for factoring multivariate
polynomials::
sage: pari('x^3 - y^3').factor()
Traceback (most recent call last):
...
PariError: sorry, factor for general polynomials is not yet implemented
TESTS::
sage: pari(2^1000+1).factor(limit=0)
doctest:...: DeprecationWarning: factor(..., lim=0) is deprecated, use an explicit limit instead
See http://trac.sagemath.org/20205 for details.
[257, 1; 1601, 1; 25601, 1; 76001, 1; 133842787352016..., 1]
"""
cdef int r
cdef GEN t0
cdef GEN cutoff
if limit == -1 and typ(self.g) == t_INT and proof:
cdef GEN g
if limit == 0:
deprecation(20205, "factor(..., lim=0) is deprecated, use an explicit limit instead")
limit = maxprime()
global factor_proven
cdef int saved_factor_proven = factor_proven
try:
if proof is not None:
factor_proven = 1 if proof else 0
pari_catch_sig_on()
# cutoff for checking true primality: 2^64 according to the
# PARI documentation ??ispseudoprime.
cutoff = mkintn(3, 1, 0, 0) # expansion of 2^64 in base 2^32: (1,0,0)
r = factorint_withproof_sage(&t0, self.g, cutoff)
z = P.new_gen(t0)
if not r:
return z
if limit >= 0:
g = boundfact(self.g, limit)
else:
return _factor_int_when_pari_factor_failed(self, z)
pari_catch_sig_on()
return P.new_gen(factor0(self.g, limit))
g = factor(self.g)
return P.new_gen(g)
finally:
factor_proven = saved_factor_proven


###########################################
Expand Down Expand Up @@ -9952,38 +9959,3 @@ cdef GEN _Vec_append(GEN v, GEN a, long n):
return w
else:
return v


cdef _factor_int_when_pari_factor_failed(x, failed_factorization):
"""
This is called by factor when PARI's factor tried to factor, got
the failed_factorization, and it turns out that one of the factors
in there is not proved prime. At this point, we don't care too much
about speed (so don't write everything below using the PARI C
library), since the probability this function ever gets called is
infinitesimal. (That said, we of course did test this function by
forcing a fake failure in the code in misc.h.)
"""
P = failed_factorization[0] # 'primes'
E = failed_factorization[1] # exponents
if len(P) == 1 and E[0] == 1:
# Major problem -- factor can't split the integer at all, but it's composite. We're stuffed.
print "BIG WARNING: The number %s wasn't split at all by PARI, but it's definitely composite."%(P[0])
print "This is probably an infinite loop..."
w = []
for i in range(len(P)):
p = P[i]
e = E[i]
if not p.isprime():
# Try to factor further -- assume this works.
F = p.factor(proof=True)
for j in range(len(F[0])):
w.append((F[0][j], F[1][j]))
else:
w.append((p, e))
m = P.matrix(len(w), 2)
for i in range(len(w)):
m[i,0] = w[i][0]
m[i,1] = w[i][1]
return m

43 changes: 0 additions & 43 deletions src/sage/libs/pari/misc.h

This file was deleted.

32 changes: 10 additions & 22 deletions src/sage/libs/pari/pari_instance.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -195,10 +195,6 @@ from sage.libs.pari.gen cimport gen, objtogen
from sage.libs.pari.handle_error cimport _pari_init_error_handling
from sage.misc.superseded import deprecated_function_alias

# so Galois groups are represented in a sane way
# See the polgalois section of the PARI users manual.
new_galois_format = 1

# real precision in decimal digits: see documentation for
# get_real_precision() and set_real_precision(). This variable is used
# in gp to set the precision of input quantities (e.g. sqrt(2)), and for
Expand Down Expand Up @@ -503,6 +499,16 @@ cdef class PariInstance(PariInstance_auto):
# (which is want we want for the PARI library interface).
GP_DATA.flags = gpd_TEST

# Ensure that Galois groups are represented in a sane way,
# see the polgalois section of the PARI users manual.
global new_galois_format
new_galois_format = 1

# By default, factor() should prove primality of returned
# factors. This not only influences the factor() function, but
# also many functions indirectly using factoring.
global factor_proven
factor_proven = 1

def debugstack(self):
r"""
Expand Down Expand Up @@ -1335,24 +1341,6 @@ cdef class PariInstance(PariInstance_auto):
pari_catch_sig_on()
return self.new_gen(gp_read_file(filename))


##############################################

def _primelimit(self):
"""
Return the number of primes already computed by PARI.
EXAMPLES::
sage: pari._primelimit()
499979
sage: pari.init_primes(600000)
sage: pari._primelimit()
599999
"""
from sage.rings.all import ZZ
return ZZ(maxprime())

def prime_list(self, long n):
"""
prime_list(n): returns list of the first n primes
Expand Down
8 changes: 4 additions & 4 deletions src/sage/rings/number_field/number_field.py
Original file line number Diff line number Diff line change
Expand Up @@ -5163,13 +5163,13 @@ def _pari_integral_basis(self, v=None, important=True):
elif self._assume_disc_small:
B = f.nfbasis(1)
elif not important:
# Trial divide the discriminant
m = self.pari_polynomial().poldisc().abs().factor(limit=0)
# Trial divide the discriminant with primes up to 10^6
m = self.pari_polynomial().poldisc().abs().factor(limit=10**6)
# Since we only need a *squarefree* factorization for
# primes with exponent 1, we need trial division up to D^(1/3)
# instead of D^(1/2).
trialdivlimit2 = pari(pari._primelimit()**2)
trialdivlimit3 = pari(pari._primelimit()**3)
trialdivlimit2 = pari(10**12)
trialdivlimit3 = pari(10**18)
if all([ p < trialdivlimit2 or (e == 1 and p < trialdivlimit3) or p.isprime() for p,e in zip(m[0],m[1]) ]):
B = f.nfbasis(fa = m)
else:
Expand Down
2 changes: 1 addition & 1 deletion src/sage/rings/polynomial/polynomial_element.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -3673,7 +3673,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
sage: f = (x+a)^50 - (a-1)^50
sage: len(factor(f))
6
sage: pari(K.discriminant()).factor(limit=0)
sage: pari(K.discriminant()).factor(limit=10^6)
[-1, 1; 3, 15; 23, 1; 887, 1; 12583, 1; 2354691439917211, 1]
sage: factor(K.discriminant())
-1 * 3^15 * 23 * 887 * 12583 * 6335047 * 371692813
Expand Down

0 comments on commit edc5ce2

Please sign in to comment.