Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adapt to precision changes in pari 2.17 #166

Merged
merged 17 commits into from
Jan 6, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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
Loading