Skip to content

Commit

Permalink
Merge pull request #166 from antonio-rojas/prec-pari-2.17
Browse files Browse the repository at this point in the history
Adapt to precision changes in pari 2.17
  • Loading branch information
dimpase authored Jan 6, 2025
2 parents 7408a05 + db4b182 commit 01a5d62
Show file tree
Hide file tree
Showing 9 changed files with 53 additions and 114 deletions.
10 changes: 3 additions & 7 deletions autogen/args.py
Original file line number Diff line number Diff line change
Expand Up @@ -304,11 +304,11 @@ def _typerepr(self):
def ctype(self):
return "long"
def always_default(self):
return "0"
return "DEFAULT_BITPREC"
def get_argument_name(self, namesiter):
return "precision"
def c_convert_code(self):
s = " {name} = prec_bits_to_words({name})\n"
s = " {name} = nbits2prec({name})\n"
return s.format(name=self.name)

class PariArgumentBitprec(PariArgumentClass):
Expand All @@ -317,13 +317,9 @@ def _typerepr(self):
def ctype(self):
return "long"
def always_default(self):
return "0"
return "DEFAULT_BITPREC"
def get_argument_name(self, namesiter):
return "precision"
def c_convert_code(self):
s = " if not {name}:\n"
s += " {name} = default_bitprec()\n"
return s.format(name=self.name)

class PariArgumentSeriesPrec(PariArgumentClass):
def _typerepr(self):
Expand Down
16 changes: 8 additions & 8 deletions autogen/doc.py
Original file line number Diff line number Diff line change
Expand Up @@ -289,7 +289,7 @@ def get_rest_doc(function):
>>> from autogen.doc import get_rest_doc
>>> print(get_rest_doc("teichmuller"))
Teichmüller character of the :math:`p`-adic number :math:`x`, i.e. the unique
:math:`(p-1)`-th root of unity congruent to :math:`x / p^{v_p(x)}` modulo :math:`p`...
:math:`(p-1)`-th root of unity congruent to :math:`x / p^{v_...(x)}` modulo :math:`p`...
::
Expand All @@ -300,24 +300,24 @@ def get_rest_doc(function):
.. MATH::
<BLANKLINE>
f(x) = \exp (-i\pi/24).\eta ((x+1)/2)/\eta (x) {such that}
j = (f^{24}-16)^3/f^{24},
j = (f^{24}-16)^.../f^{24},
<BLANKLINE>
where :math:`j` is the elliptic :math:`j`-invariant (see the function :literal:`ellj`).
If :math:`flag = 1`, returns
<BLANKLINE>
.. MATH::
<BLANKLINE>
f_1(x) = \eta (x/2)/\eta (x) {such that}
j = (f_1^{24}+16)^3/f_1^{24}.
f_...(x) = \eta (x/2)/\eta (x) {such that}
j = (f_...^{24}+16)^.../f_...^{24}.
<BLANKLINE>
Finally, if :math:`flag = 2`, returns
<BLANKLINE>
.. MATH::
<BLANKLINE>
f_2(x) = \sqrt{2}\eta (2x)/\eta (x) {such that}
j = (f_2^{24}+16)^3/f_2^{24}.
f_...(x) = \sqrt{2}\eta (2x)/\eta (x) {such that}
j = (f_...^{24}+16)^.../f_...^{24}.
<BLANKLINE>
Note the identities :math:`f^8 = f_1^8+f_2^8` and :math:`ff_1f_2 = \sqrt2`.
Note the identities :math:`f^... = f_...^...+f_...^...` and :math:`ff_...f_... = \sqrt2`.
::
Expand All @@ -333,7 +333,7 @@ def get_rest_doc(function):
.. MATH::
<BLANKLINE>
\sum
(x_i or y_i) 2^i
(x_... or y_...) 2^...
<BLANKLINE>
See ``bitand`` (in the PARI manual) for the behavior for negative arguments.
"""
Expand Down
4 changes: 2 additions & 2 deletions autogen/generator.py
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ def handle_pari_function(self, function, cname, prototype="", help="", obsolete=
... help=r"bnfinit(P,{flag=0},{tech=[]}): compute...",
... **{"class":"basic", "section":"number_fields"})
GEN bnfinit0(GEN, long, GEN, long)
def bnfinit(P, long flag=0, tech=None, long precision=0):
def bnfinit(P, long flag=0, tech=None, long precision=DEFAULT_BITPREC):
...
cdef bint _have_tech = (tech is not None)
if _have_tech:
Expand All @@ -149,7 +149,7 @@ def bnfinit(P, long flag=0, tech=None, long precision=0):
cdef GEN _tech = NULL
if _have_tech:
_tech = (<Gen>tech).g
precision = prec_bits_to_words(precision)
precision = nbits2prec(precision)
cdef GEN _ret = bnfinit0(_P, flag, _tech, precision)
return new_gen(_ret)
<BLANKLINE>
Expand Down
2 changes: 1 addition & 1 deletion autogen/parser.py
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ def parse_prototype(proto, help, initial_args=[]):
>>> parse_prototype(proto, help)
([GEN x, GEN* r=NULL], GEN)
>>> parse_prototype("lp", "foo()", [str("TEST")])
(['TEST', prec precision=0], long)
(['TEST', prec precision=DEFAULT_BITPREC], long)
"""
# Use the help string just for the argument names.
# "names" should be an iterator over the argument names.
Expand Down
35 changes: 19 additions & 16 deletions cypari2/gen.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -66,8 +66,7 @@ from .types cimport *
from .string_utils cimport to_string, to_bytes
from .paripriv cimport *
from .convert cimport PyObject_AsGEN, gen_to_integer
from .pari_instance cimport (prec_bits_to_words,
default_bitprec, get_var)
from .pari_instance cimport DEFAULT_BITPREC, get_var
from .stack cimport (new_gen, new_gens2, new_gen_noclear,
clone_gen, clear_stack, reset_avma,
remove_from_pari_stack, move_gens_to_heap)
Expand Down Expand Up @@ -647,7 +646,7 @@ cdef class Gen(Gen_base):
if m is not None:
t0 = t0.Mod(m)
sig_on()
return new_gen(gpow(t0.g, t1.g, prec_bits_to_words(0)))
return new_gen(gpow(t0.g, t1.g, nbits2prec(DEFAULT_BITPREC)))

def __neg__(self):
sig_on()
Expand Down Expand Up @@ -902,6 +901,7 @@ cdef class Gen(Gen_base):
>>> x = pari('x')
>>> pari.setrand(1)
>>> (x**2 - 65).bnfinit().bnf_get_fu()
[Mod(x - 8, x^2 - 65)]
>>> (x**4 - x**2 + 1).bnfinit().bnf_get_fu()
Expand Down Expand Up @@ -951,6 +951,7 @@ cdef class Gen(Gen_base):
>>> import warnings
>>> with warnings.catch_warnings(record=True) as w:
... warnings.simplefilter('always')
... pari.setrand(1)
... funits = (x**2 - 65).bnfinit().bnfunit()
... assert len(w) == 1
... assert issubclass(w[0].category, DeprecationWarning)
Expand Down Expand Up @@ -2911,7 +2912,7 @@ cdef class Gen(Gen_base):
sig_on()
return new_gen(bernfrac(self))

def bernreal(self, unsigned long precision=0):
def bernreal(self, unsigned long precision=DEFAULT_BITPREC):
r"""
The Bernoulli number `B_x`, as for the function bernfrac,
but `B_x` is returned as a real number (with the current
Expand All @@ -2926,9 +2927,9 @@ cdef class Gen(Gen_base):
54.9711779448622
"""
sig_on()
return new_gen(bernreal(self, prec_bits_to_words(precision)))
return new_gen(bernreal(self, nbits2prec(precision)))

def besselk(nu, x, unsigned long precision=0):
def besselk(nu, x, unsigned long precision=DEFAULT_BITPREC):
"""
nu.besselk(x): K-Bessel function (modified Bessel function
of the second kind) of index nu, which can be complex, and argument
Expand Down Expand Up @@ -2963,9 +2964,9 @@ cdef class Gen(Gen_base):
"""
cdef Gen t0 = objtogen(x)
sig_on()
return new_gen(kbessel(nu.g, t0.g, prec_bits_to_words(precision)))
return new_gen(kbessel(nu.g, t0.g, nbits2prec(precision)))

def eint1(x, long n=0, unsigned long precision=0):
def eint1(x, long n=0, unsigned long precision=DEFAULT_BITPREC):
r"""
x.eint1(n): exponential integral E1(x):
Expand All @@ -2991,13 +2992,14 @@ cdef class Gen(Gen_base):
"""
sig_on()
if n <= 0:
return new_gen(eint1(x.g, prec_bits_to_words(precision)))
return new_gen(eint1(x.g, nbits2prec(precision)))
else:
return new_gen(veceint1(x.g, stoi(n), prec_bits_to_words(precision)))
return new_gen(veceint1(x.g, stoi(n), nbits2prec(precision)))

log_gamma = Gen_base.lngamma

def polylog(x, long m, long flag=0, unsigned long precision=0):
def polylog(x, long m, long flag=0,
unsigned long precision=DEFAULT_BITPREC):
"""
x.polylog(m,flag=0): m-th polylogarithm of x. flag is optional, and
can be 0: default, 1: D_m -modified m-th polylog of x, 2:
Expand Down Expand Up @@ -3026,9 +3028,9 @@ cdef class Gen(Gen_base):
-0.400459056163451
"""
sig_on()
return new_gen(polylog0(m, x.g, flag, prec_bits_to_words(precision)))
return new_gen(polylog0(m, x.g, flag, nbits2prec(precision)))

def sqrtn(x, n, unsigned long precision=0):
def sqrtn(x, n, unsigned long precision=DEFAULT_BITPREC):
r"""
x.sqrtn(n): return the principal branch of the n-th root of x,
i.e., the one such that
Expand Down Expand Up @@ -3090,7 +3092,7 @@ cdef class Gen(Gen_base):
cdef GEN ans, zetan
cdef Gen t0 = objtogen(n)
sig_on()
ans = gsqrtn(x.g, t0.g, &zetan, prec_bits_to_words(precision))
ans = gsqrtn(x.g, t0.g, &zetan, nbits2prec(precision))
return new_gens2(ans, zetan)

def ffprimroot(self):
Expand Down Expand Up @@ -4531,7 +4533,8 @@ cdef class Gen(Gen_base):
g = polint(self.g, t0.g, t1.g, &dy)
return new_gens2(g, dy)

def ellwp(self, z='z', long n=20, long flag=0, unsigned long precision=0):
def ellwp(self, z='z', long n=20, long flag=0,
unsigned long precision=DEFAULT_BITPREC):
"""
Return the value or the series expansion of the Weierstrass
`P`-function at `z` on the lattice `self` (or the lattice
Expand Down Expand Up @@ -4609,7 +4612,7 @@ cdef class Gen(Gen_base):
elif typ(g0) == t_RFRAC:
g0 = rfrac_to_ser(g0, n+4)

cdef GEN r = ellwp0(self.g, g0, flag, prec_bits_to_words(precision))
cdef GEN r = ellwp0(self.g, g0, flag, nbits2prec(precision))
if flag == 1 and have_ellwp_flag1_bug():
# Work around ellwp() bug: double the second element
set_gel(r, 2, gmulgs(gel(r, 2), 2))
Expand Down
8 changes: 7 additions & 1 deletion cypari2/pari_instance.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,16 @@ cimport cython

from .gen cimport Gen

cpdef long prec_bits_to_words(unsigned long prec_in_bits) noexcept
# DEPRECATED INTERNAL FUNCTION used (incorrectly) in sagemath < 10.5
cpdef long prec_words_to_bits(long prec_in_words) noexcept
cpdef long default_bitprec() noexcept

cdef extern from *:
"""
#define DEFAULT_BITPREC prec2nbits(DEFAULTPREC)
"""
long DEFAULT_BITPREC

cdef class Pari_auto:
pass

Expand Down
89 changes: 10 additions & 79 deletions cypari2/pari_instance.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -235,9 +235,9 @@ Verify that ``nfroots()`` (which has an unusual signature with a
non-default argument following a default argument) works:
>>> pari.nfroots(x='x^4 - 1')
[-1, 1]
[-1, 1]...
>>> pari.nfroots(pari.nfinit('t^2 + 1'), "x^4 - 1")
[-1, 1, Mod(-t, t^2 + 1), Mod(t, t^2 + 1)]
[-1, 1, Mod(-t, t^2 + 1), Mod(t, t^2 + 1)]...
Reset default precision for the following tests:
Expand Down Expand Up @@ -299,10 +299,6 @@ from .stack cimport (new_gen, new_gen_noclear, clear_stack,
from .handle_error cimport _pari_init_error_handling
from .closure cimport _pari_init_closure

# Default precision (in PARI words) for the PARI library interface,
# when no explicit precision is given and the inputs are exact.
cdef long prec = prec_bits_to_words(53)


#################################################################
# conversions between various real precision models
Expand Down Expand Up @@ -341,36 +337,10 @@ def prec_dec_to_bits(long prec_in_dec):
return int(prec_in_dec*log_10 + 1.0) # Add one to round up


cpdef long prec_bits_to_words(unsigned long prec_in_bits) noexcept:
r"""
Convert from precision expressed in bits to pari real precision
expressed in words. Note: this rounds up to the nearest word,
adjusts for the two codewords of a pari real, and is
architecture-dependent.
Examples:
>>> from cypari2.pari_instance import prec_bits_to_words
>>> import sys
>>> bitness = '64' if sys.maxsize > (1 << 32) else '32'
>>> prec_bits_to_words(70) == (5 if bitness == '32' else 4)
True
>>> ans32 = [(32, 3), (64, 4), (96, 5), (128, 6), (160, 7), (192, 8), (224, 9), (256, 10)]
>>> ans64 = [(32, 3), (64, 3), (96, 4), (128, 4), (160, 5), (192, 5), (224, 6), (256, 6)]
>>> [(32*n, prec_bits_to_words(32*n)) for n in range(1, 9)] == (ans32 if bitness == '32' else ans64)
True
"""
if not prec_in_bits:
return prec
cdef unsigned long wordsize = BITS_IN_LONG

# This equals ceil(prec_in_bits/wordsize) + 2
return (prec_in_bits - 1)//wordsize + 3


cpdef long prec_words_to_bits(long prec_in_words) noexcept:
r"""
Deprecated internal function. Used (incorrectly) in sagemath < 10.5.
Convert from pari real precision expressed in words to precision
expressed in bits. Note: this adjusts for the two codewords of a
pari real, and is architecture-dependent.
Expand All @@ -379,6 +349,8 @@ cpdef long prec_words_to_bits(long prec_in_words) noexcept:
>>> from cypari2.pari_instance import prec_words_to_bits
>>> import sys
>>> import warnings
>>> warnings.simplefilter("ignore")
>>> bitness = '64' if sys.maxsize > (1 << 32) else '32'
>>> prec_words_to_bits(10) == (256 if bitness == '32' else 512)
True
Expand All @@ -388,6 +360,9 @@ cpdef long prec_words_to_bits(long prec_in_words) noexcept:
>>> [(n, prec_words_to_bits(n)) for n in range(3, 10)] == (ans32 if bitness == '32' else ans64)
True
"""
from warnings import warn
warn("'prec_words_to_bits` in cypari2 is internal and deprecated",
DeprecationWarning)
# see user's guide to the pari library, page 10
return (prec_in_words - 2) * BITS_IN_LONG

Expand All @@ -402,51 +377,7 @@ cpdef long default_bitprec() noexcept:
>>> default_bitprec()
64
"""
return (prec - 2) * BITS_IN_LONG


def prec_dec_to_words(long prec_in_dec):
r"""
Convert from precision expressed in decimal to precision expressed
in words. Note: this rounds up to the nearest word, adjusts for the
two codewords of a pari real, and is architecture-dependent.
Examples:
>>> from cypari2.pari_instance import prec_dec_to_words
>>> import sys
>>> bitness = '64' if sys.maxsize > (1 << 32) else '32'
>>> prec_dec_to_words(38) == (6 if bitness == '32' else 4)
True
>>> ans32 = [(10, 4), (20, 5), (30, 6), (40, 7), (50, 8), (60, 9), (70, 10), (80, 11)]
>>> ans64 = [(10, 3), (20, 4), (30, 4), (40, 5), (50, 5), (60, 6), (70, 6), (80, 7)] # 64-bit
>>> [(n, prec_dec_to_words(n)) for n in range(10, 90, 10)] == (ans32 if bitness == '32' else ans64)
True
"""
return prec_bits_to_words(prec_dec_to_bits(prec_in_dec))


def prec_words_to_dec(long prec_in_words):
r"""
Convert from precision expressed in words to precision expressed in
decimal. Note: this adjusts for the two codewords of a pari real,
and is architecture-dependent.
Examples:
>>> from cypari2.pari_instance import prec_words_to_dec
>>> import sys
>>> bitness = '64' if sys.maxsize > (1 << 32) else '32'
>>> prec_words_to_dec(5) == (28 if bitness == '32' else 57)
True
>>> ans32 = [(3, 9), (4, 19), (5, 28), (6, 38), (7, 48), (8, 57), (9, 67)]
>>> ans64 = [(3, 19), (4, 38), (5, 57), (6, 77), (7, 96), (8, 115), (9, 134)]
>>> [(n, prec_words_to_dec(n)) for n in range(3, 10)] == (ans32 if bitness == '32' else ans64)
True
"""
return prec_bits_to_dec(prec_words_to_bits(prec_in_words))
return DEFAULT_BITPREC


# Callbacks from PARI to print stuff using sys.stdout.write() instead
Expand Down
2 changes: 2 additions & 0 deletions cypari2/paridecl.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -5107,6 +5107,7 @@ cdef extern from *: # PARI headers already included by types.pxd
void killblock(GEN x)
GEN leading_coeff(GEN x)
void lg_increase(GEN x)
long lg2prec(long x)
long lgcols(GEN x)
long lgpol(GEN x)
GEN matpascal(long n)
Expand Down Expand Up @@ -5196,6 +5197,7 @@ cdef extern from *: # PARI headers already included by types.pxd
GEN polx_zx(long sv)
GEN powii(GEN x, GEN n)
GEN powIs(long n)
long prec2lg(long x)
long prec2nbits(long x)
double prec2nbits_mul(long x, double y)
long prec2ndec(long x)
Expand Down
Loading

0 comments on commit 01a5d62

Please sign in to comment.