Pari/GP compare integer with real

Question

I want to test if an integer is a perfect power in Pari-gp. The test sqrt(n)==floor(sqrt(n)) works fine for testing squares, but it fails for every other power: sqrtn(n,k)==floor(sqrtn(n,k)) with k >=3.

I think it maybe since one number is real and the other one is integer. Still the test works for squares. What am I doing wrong?

K.B. · Answer 1 · 2023-04-03T20:44:55.780

You should directly use ispower(), and ispower(, k) for perfect k-th powers.
The problem with the factor approach is that it will be exceedingly slow for large inputs, whereas ispower runs in polynomial time.
The test sqrtn(n,k)==floor(sqrtn(n,k)) can't be reliable since PARI doesn't guarantee exact rounding: even if the value of sqrtn() is mathematically an exact integer, PARI may return any real number which is a little larger or a little less than the exact answer. You will be a little better off with round than with floor. It's still not reliable though, here's a solution that would more or less work

y = round(sqrtn(n, k)); y^k == n

(provided realprecision is large enough). But ispower starts with modular tests that will make it immensely faster whenever the number is not a k-th power.

Great answer! Just note that you can have some list styling by putting `1.` instead of `1)`! — Hugo Dozois, Apr 05 '13 at 19:18

Andrew Lamoureux · Accepted Answer · 2012-04-17T22:22:41.093

2

Prime powers have factorizations involving only one base (the prime factor itself). So a better way to perform the test would be:

isPrimePower(x) = {
  matsize(factor(x))[1]==1
}

Here's a test for the first 10 numbers:

for(i=0,10,print(i,"->",isPrimePower(i)))
0->1  (yes, p^0)
1->0
2->1  (yes, 2^1)
3->1  (yes, 3^1)
4->1  (yes, 2^2)
5->1  (yes, 5^1)
6->0
7->1  (yes, 7^1)
8->1  (yes, 2^3)
9->1  (yes, 3^3)
10->0

For composite bases, I'll have to assume you mean that a perfect power factors into some single base raised to an exponent e >= 2. Otherwise any n = n^1. Even now we have corner case of 1 since 1=1^k.

isPerfectPower(x) = { 
  local(e);
  if(x==1, return(0));
  factors = factor(x);
  e = factors[1,2];
  if(e < 2, return(0));
  for(i=2,matsize(factors)[1],
      if(factors[i,2] != e, return(0));
  );
  return(1);
}

And the test again:

> for(i=1,20,print(i,"->",isPerfectPower(i)))
1->0
2->0
3->0
4->1
5->0
6->0
7->0
8->1
9->1
10->0
11->0
12->0
13->0
14->0
15->0
16->1
17->0
18->0
19->0
20->0

edited Apr 17 '12 at 22:22

answered Apr 17 '12 at 19:55

Andrew Lamoureux

70
7

Thank you for that, but it only works for prime powers. I still don't get why the test with roots I wrote doesn't work. – Beni Bogosel Apr 17 '12 at 20:00
Testing with 2.5.0 here, your test does work for k>3: `(18:13) gp > sqrtn(8,3) == floor(sqrtn(8,3))` returns 1 and `(18:13) gp > sqrtn(16,4) == floor(sqrtn(16,4))` returns 1. Editing my answer to work for composite powers... – Andrew Lamoureux Apr 17 '12 at 22:13
Try `sqrtn(216,3)==floor(sqrtn(216,3))`, and greater cubes. I also have 2.5.0. – Beni Bogosel Apr 18 '12 at 08:21
If two differing types are being compared, pari attempts to add the first to the negation of the second and test for zero. Apparently the t_REAL result from sqrtn() can sometimes be a number with extremely small fractional parts in addition to the whole number part you are expecting. Try for example `sqrtn(216,3)-6` vs `sqrtn(8,3)-2`. – Andrew Lamoureux Apr 18 '12 at 21:51

score 0 · Answer 3 · answered Dec 22 '12 at 01:02

You could test frac(sqrtn(261,3)+epsilon) < 2*epsilon where epsilon is some very small number you think acceptable, such as 1E-15

I write it that why and not just "... < epsilon" because sometimes you get 4.0000000001 but you could also get 3.9999999999.

Pari/GP compare integer with real

3 Answers3