Opened 8 years ago

Closed 7 years ago

#12116 closed enhancement (fixed)

perfect_power for integers

Reported by: roed Owned by: AlexGhitza
Priority: major Milestone: sage-5.8
Component: basic arithmetic Keywords:
Cc: jpflori Merged in: sage-5.8.beta3
Authors: David Roe Reviewers: David Loeffler, Aly Deines, Simon Spicer
Report Upstream: Fixed upstream, but not in a stable release. Work issues:
Branch: Commit:
Dependencies: #10596, #12363, #12638 Stopgaps:

Description (last modified by roed)

Currently there are two functions for testing whether an integer is a perfect power:

sage: 8.is_power()
True
sage: 8.is_perfect_power()
True
  • This patch deprecates is_power.
  • It implements a new function, perfect_power, wrapping PARI's ispower function, expressing an integer in the form ab with b maximal.
  • It adds cross references between is_perfect_power, is_power_of, is_prime_power and perfect_power.

This revealed three bugs in PARI's ispower() function:

  1. Two issues reported and fixed at http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1259
  2. One more issue reported and fixed at http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1302 (duplicate: http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1303)

Apply trac_12116-rebase4.patch

Attachments (4)

trac_12116-rebase.patch (9.1 KB) - added by roed 7 years ago.
Rebased against 5.8.beta0
trac_12116_rebase2.patch (9.1 KB) - added by spice 7 years ago.
Flips order of returned pair in self.perfect_power()
trac_12116-rebase3.patch (9.1 KB) - added by roed 7 years ago.
trac_12116-rebase4.patch (9.1 KB) - added by spice 7 years ago.
Changed a line to unbreak factor_aurifeuillian() in factorint.pyx

Download all attachments as: .zip

Change History (37)

comment:1 Changed 8 years ago by roed

  • Component changed from PLEASE CHANGE to basic arithmetic
  • Owner changed from tbd to AlexGhitza
  • Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. Little or no feedback.

The bug is being tracked by Pari: their ticket number is 1259.

comment:2 Changed 8 years ago by roed

  • Status changed from new to needs_review

comment:3 Changed 8 years ago by roed

  • Report Upstream changed from Reported upstream. Little or no feedback. to Fixed upstream, but not in a stable release.

Here's the response from the Pari developers:

The 2 problems turned out to be independant. They are now fixed in the
'testing' branch, in the following 2 commits:

commit f53eb2a7c79922ac51423e75044722628893159e
Author: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
Date:   Sun Dec 4 17:43:19 2011 +0100

   46- ispower(1 / n) return a wrong result [#1259]

commit 74aef498d380979211b99d64f4f4c92b8bac853f
Author: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
Date:   Sun Dec 4 17:28:14 2011 +0100

   45- ispower(x < 0) could return an even value ! [#1259]

That was quick. :-)

comment:4 Changed 8 years ago by jason

I'm curious why you didn't use the Sphinx seealso directive: http://sphinx.pocoo.org/markup/para.html#directive-seealso

comment:5 Changed 8 years ago by roed

No good reason.

Apply 12116.2.patch only.

comment:6 Changed 8 years ago by davidloeffler

  • Dependencies set to 10596
  • Reviewers set to David Loeffler

Sadly roed's patch conflicts with #10596 and thus doesn't apply to the current Sage beta. I've just uploaded a rebased patch based on 5.0.beta7. This makes no changes to roed's code other than some minor sphinx reformatting (including using the new :trac: directive to point to this ticket). David: if you're happy with my reworking of your patch, you can set this to positive review.

comment:7 Changed 8 years ago by davidloeffler

  • Dependencies changed from 10596 to #10596

comment:8 Changed 8 years ago by roed

  • Status changed from needs_review to positive_review

Thanks for rebasing! I looked it over and it looks good to me.

comment:9 Changed 8 years ago by jdemeyer

  • Dependencies changed from #10596 to #10596, #12363
  • Description modified (diff)
  • Report Upstream changed from Fixed upstream, but not in a stable release. to Completely fixed; Fix reported upstream
  • Status changed from positive_review to needs_work

The upstream fix is in Sage now (#12363). Could you remove the workaround code from the patch?

comment:10 Changed 8 years ago by davidloeffler

  • Report Upstream changed from Completely fixed; Fix reported upstream to None of the above - read trac for reasoning.
  • Status changed from needs_work to needs_review

Can we have a new option for the "Report Upstream" menu, "upstream claim to have fixed it but they apparently haven't"? PARI apparently believes that 2^3 = -8 :

                                      GP/PARI CALCULATOR Version 2.5.1 (development git-5c5e253)
                                     amd64 running linux (x86-64/GMP-5.0.1 kernel) 64-bit version
                                       compiled: Mar  5 2012, gcc-4.4.3 (Ubuntu 4.4.3-4ubuntu5) 
                                            (readline v6.2 enabled, extended help enabled)

                                                Copyright (C) 2000-2011 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?12 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 500509
? ispower(-8, , &n)
%1 = 3
? n
%2 = 2

So I suggest we merge this patch as it is, and send another bug report to PARI -- roed, can you do the honours?

comment:11 Changed 8 years ago by jdemeyer

  • Report Upstream changed from None of the above - read trac for reasoning. to Reported upstream. Little or no feedback.
  • Status changed from needs_review to needs_info

I reported the new issue upstream. I would postpone merging this patch until we hear from upstream. I prefer to have the issue fixed upstream and not apply a work-around in Sage.

comment:12 Changed 8 years ago by davidloeffler

Oops -- I also just filed a bug report, after hearing from David R that he's currently busy conference-ing. Oh well, no harm is done if they get two reports of the same bug.

comment:13 Changed 8 years ago by jdemeyer

  • Description modified (diff)

comment:14 Changed 8 years ago by davidloeffler

  • Description modified (diff)

comment:15 Changed 8 years ago by jdemeyer

  • Dependencies changed from #10596, #12363 to #10596, #12363, #12638
  • Report Upstream changed from Reported upstream. Little or no feedback. to Fixed upstream, but not in a stable release.
  • Status changed from needs_info to needs_review

The PARI patch has been, I will provide a patched spkg at #12638.

comment:16 Changed 8 years ago by jdemeyer

  • Status changed from needs_review to needs_work

comment:17 Changed 8 years ago by jdemeyer

  • Description modified (diff)

Patched PARI spkg at #12638 ready for review.

comment:18 Changed 8 years ago by SimonKing

Is still work needed? Or a review?

comment:19 Changed 8 years ago by roed

I think Jeroen suggested that this ticket should depend on #12638 and remove the PARI workaround. I'll try to put up a patch at some point.

comment:20 Changed 7 years ago by roed

  • Status changed from needs_work to needs_review

I've updated this to remove the PARI workaround and import deprecated_function_alias from sage.misc.superceded rather than sage.misc.misc.

comment:21 Changed 7 years ago by jpflori

  • Cc jpflori added

comment:22 follow-up: Changed 7 years ago by kcrisman

I like this a lot.

Silly questions:

  • cdef int exp = 1 in the second hunk - since
    sage: type(4.perfect_power()[0])
    sage.rings.integer.Integer
    
    and you now reassign exp to that, does that line make any sense at all? It's certainly not an int.
  • Do we need to deprecate base_exponent? It would be easy enough, basically same code as deprecating is_power except reverse the order of the output.

comment:23 in reply to: ↑ 22 Changed 7 years ago by roed

  • Status changed from needs_review to needs_work
  • cdef int exp = 1 in the second hunk - since
    sage: type(4.perfect_power()[0])
    sage.rings.integer.Integer
    
    and you now reassign exp to that, does that line make any sense at all? It's certainly not an int.

I think Cython is smart enough to figure out what's going on here since Integer has an __int__ method. You can see in the generated C code for exp, b = b.perfect_power() (leaving out lots of error checking and reference counting):

__pyx_t_1 = PyObject_GetAttr(__pyx_v_b, __pyx_n_s__perfect_power);
__pyx_t_7 = PyObject_Call(__pyx_t_1, ((PyObject *)__pyx_empty_tuple), NULL);
PyObject* sequence = __pyx_t_7;
Py_ssize_t size = Py_SIZE(sequence);
__pyx_t_1 = PyTuple_GET_ITEM(sequence, 0);
__pyx_t_11 = __Pyx_PyInt_AsInt(__pyx_t_1);
__pyx_v_exp = __pyx_t_11;
  • Do we need to deprecate base_exponent? It would be easy enough, basically same code as deprecating is_power except reverse the order of the output.

It's probably a good idea. I'll do it shortly.

Changed 7 years ago by roed

Rebased against 5.8.beta0

comment:24 Changed 7 years ago by roed

  • Status changed from needs_work to needs_review

Alright, fixed.

comment:25 Changed 7 years ago by aly.deines

  • Reviewers changed from David Loeffler to David Loeffler, Aly Deines
  • Status changed from needs_review to positive_review

Changed 7 years ago by spice

Flips order of returned pair in self.perfect_power()

comment:26 Changed 7 years ago by spice

  • Status changed from positive_review to needs_work

Everything checks out. Patch applies, tests pass. However:

def perfect_power(self): 
    r""" 	         
    Returns ``(a, b)``, where this integer is `b^a` and `a` is maximal.
    ...

So

sage: a = 11^14
sage: a.perfect_power()
(14, 11)

I know the code just wraps PARI which returns the integers in this order, but this seems an unnatural way to order the returned pair.

I've uploaded a modified version of the patch so that we get

sage: a = 11^14
sage: a.perfect_power()
(11, 14)

which feels a lot more natural.

The new patch also fixes a possible formatting error in the documentation of self.is_perfect_power():

Returns ``True`` if self is a perfect power, ie if there exist integers
a and b, `b > 1` with `\texttt{self} = a^b`.

becomes

Returns ``True`` if ``self`` is a perfect power, ie if there exist integers
`a` and `b`, `b > 1` with ``self`` `= a^b`.

comment:27 Changed 7 years ago by spice

  • Status changed from needs_work to needs_review

comment:28 Changed 7 years ago by spice

  • Reviewers changed from David Loeffler, Aly Deines to David Loeffler, Aly Deines, Simon Spicer

Changed 7 years ago by roed

comment:29 Changed 7 years ago by roed

Looks fine to me. I just updated the deprecated_function_alias.

Changed 7 years ago by spice

Changed a line to unbreak factor_aurifeuillian() in factorint.pyx

comment:30 Changed 7 years ago by spice

We missed another line. In self.factor_aurifeuillian() in factorint.pyx, we have

...
b = n + x
exp, b = b.perfect_power()
if exp > 1:
...

This must now become

...
b = n + x
b, exp = b.perfect_power()
if exp > 1:
...

otherwise this method breaks with the new ordering of values returned by perfect_power(). I've uploaded a new patch that fixes this.

Hopefully that should be it.

comment:31 Changed 7 years ago by roed

  • Description modified (diff)
  • Status changed from needs_review to positive_review

Looks good to me.

comment:32 Changed 7 years ago by roed

For patchbot:

Apply trac_12116-rebase4.patch

comment:33 Changed 7 years ago by jdemeyer

  • Merged in set to sage-5.8.beta3
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.