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

is_perfect_power implementation #137

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open
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
38 changes: 20 additions & 18 deletions ulong_extras.h
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ FLINT_DLL mp_limb_t n_flog(mp_limb_t n, mp_limb_t b);

FLINT_DLL mp_limb_t n_clog(mp_limb_t n, mp_limb_t b);

static __inline__
static __inline__
double n_precompute_inverse(mp_limb_t n)
{
return (double) 1 / (double) n;
Expand All @@ -176,29 +176,29 @@ FLINT_DLL mp_limb_t n_mod_precomp(mp_limb_t a, mp_limb_t n, double ninv);

FLINT_DLL mp_limb_t n_mod2_precomp(mp_limb_t a, mp_limb_t n, double ninv);

FLINT_DLL mp_limb_t n_divrem2_precomp(mp_limb_t * q, mp_limb_t a,
FLINT_DLL mp_limb_t n_divrem2_precomp(mp_limb_t * q, mp_limb_t a,
mp_limb_t n, double npre);

FLINT_DLL mp_limb_t n_mod2_preinv(mp_limb_t a, mp_limb_t n, mp_limb_t ninv);

FLINT_DLL mp_limb_t n_ll_mod_preinv(mp_limb_t a_hi, mp_limb_t a_lo,
FLINT_DLL mp_limb_t n_ll_mod_preinv(mp_limb_t a_hi, mp_limb_t a_lo,
mp_limb_t n, mp_limb_t ninv);

FLINT_DLL mp_limb_t n_lll_mod_preinv(mp_limb_t a_hi, mp_limb_t a_mi,
FLINT_DLL mp_limb_t n_lll_mod_preinv(mp_limb_t a_hi, mp_limb_t a_mi,
mp_limb_t a_lo, mp_limb_t n, mp_limb_t ninv);

FLINT_DLL mp_limb_t n_mulmod_precomp(mp_limb_t a, mp_limb_t b,
FLINT_DLL mp_limb_t n_mulmod_precomp(mp_limb_t a, mp_limb_t b,
mp_limb_t n, double ninv);

FLINT_DLL mp_limb_t n_mulmod2_preinv(mp_limb_t a, mp_limb_t b,
FLINT_DLL mp_limb_t n_mulmod2_preinv(mp_limb_t a, mp_limb_t b,
mp_limb_t n, mp_limb_t ninv);

FLINT_DLL mp_limb_t n_mulmod_preinv(mp_limb_t a, mp_limb_t b,
FLINT_DLL mp_limb_t n_mulmod_preinv(mp_limb_t a, mp_limb_t b,
mp_limb_t n, mp_limb_t ninv, ulong norm);

FLINT_DLL mp_limb_t n_powmod_ui_precomp(mp_limb_t a, mp_limb_t exp, mp_limb_t n, double npre);

FLINT_DLL mp_limb_t n_powmod_precomp(mp_limb_t a,
FLINT_DLL mp_limb_t n_powmod_precomp(mp_limb_t a,
mp_limb_signed_t exp, mp_limb_t n, double npre);

static __inline__
Expand All @@ -209,13 +209,13 @@ mp_limb_t n_powmod(mp_limb_t a, mp_limb_signed_t exp, mp_limb_t n)
return n_powmod_precomp(a, exp, n, npre);
}

FLINT_DLL mp_limb_t n_powmod2_preinv(mp_limb_t a,
FLINT_DLL mp_limb_t n_powmod2_preinv(mp_limb_t a,
mp_limb_signed_t exp, mp_limb_t n, mp_limb_t ninv);

FLINT_DLL mp_limb_t n_powmod2_ui_preinv(mp_limb_t a, mp_limb_t exp,
mp_limb_t n, mp_limb_t ninv);

FLINT_DLL mp_limb_t n_powmod_ui_preinv(mp_limb_t a, mp_limb_t exp, mp_limb_t n,
FLINT_DLL mp_limb_t n_powmod_ui_preinv(mp_limb_t a, mp_limb_t exp, mp_limb_t n,
mp_limb_t ninv, ulong norm);

static __inline__
Expand Down Expand Up @@ -246,9 +246,9 @@ mp_limb_t n_negmod(mp_limb_t x, mp_limb_t n)

FLINT_DLL mp_limb_t n_sqrtmod(mp_limb_t a, mp_limb_t p);

FLINT_DLL slong n_sqrtmod_2pow(mp_limb_t ** sqrt, mp_limb_t a, slong exp);
FLINT_DLL slong n_sqrtmod_2pow(mp_limb_t ** sqrt, mp_limb_t a, slong exp);

FLINT_DLL slong n_sqrtmod_primepow(mp_limb_t ** sqrt, mp_limb_t a,
FLINT_DLL slong n_sqrtmod_primepow(mp_limb_t ** sqrt, mp_limb_t a,
mp_limb_t p, slong exp);

FLINT_DLL slong n_sqrtmodn(mp_limb_t ** sqrt, mp_limb_t a, n_factor_t * fac);
Expand Down Expand Up @@ -289,6 +289,8 @@ FLINT_DLL mp_limb_t n_cbrtrem(mp_limb_t* remainder, mp_limb_t n);

FLINT_DLL int n_is_perfect_power235(mp_limb_t n);

FLINT_DLL int n_is_perfect_power(mp_limb_t * r, mp_limb_t x);

FLINT_DLL int n_is_oddprime_small(mp_limb_t n);

FLINT_DLL int n_is_oddprime_binary(mp_limb_t n);
Expand All @@ -301,10 +303,10 @@ FLINT_DLL int n_is_probabprime_lucas(mp_limb_t n);

FLINT_DLL int n_is_probabprime_BPSW(mp_limb_t n);

FLINT_DLL int n_is_strong_probabprime_precomp(mp_limb_t n,
FLINT_DLL int n_is_strong_probabprime_precomp(mp_limb_t n,
double npre, mp_limb_t a, mp_limb_t d);

FLINT_DLL int n_is_strong_probabprime2_preinv(mp_limb_t n,
FLINT_DLL int n_is_strong_probabprime2_preinv(mp_limb_t n,
mp_limb_t ninv, mp_limb_t a, mp_limb_t d);

FLINT_DLL int n_is_probabprime(mp_limb_t n);
Expand Down Expand Up @@ -335,16 +337,16 @@ void n_factor_init(n_factor_t * factors)

FLINT_DLL void n_factor_insert(n_factor_t * factors, mp_limb_t p, ulong exp);

FLINT_DLL mp_limb_t n_factor_trial_range(n_factor_t * factors,
FLINT_DLL mp_limb_t n_factor_trial_range(n_factor_t * factors,
mp_limb_t n, ulong start, ulong num_primes);

FLINT_DLL mp_limb_t n_factor_trial_partial(n_factor_t * factors, mp_limb_t n,
FLINT_DLL mp_limb_t n_factor_trial_partial(n_factor_t * factors, mp_limb_t n,
mp_limb_t * prod, ulong num_primes, mp_limb_t limit);

FLINT_DLL mp_limb_t n_factor_trial(n_factor_t * factors,
FLINT_DLL mp_limb_t n_factor_trial(n_factor_t * factors,
mp_limb_t n, ulong num_primes);

FLINT_DLL mp_limb_t n_factor_partial(n_factor_t * factors,
FLINT_DLL mp_limb_t n_factor_partial(n_factor_t * factors,
mp_limb_t n, mp_limb_t limit, int proved);

FLINT_DLL mp_limb_t n_factor_power235(ulong *exp, mp_limb_t n);
Expand Down
6 changes: 6 additions & 0 deletions ulong_extras/doc/ulong_extras.txt
Original file line number Diff line number Diff line change
Expand Up @@ -908,6 +908,12 @@ int n_is_perfect_power235(mp_limb_t n)
root can be taken, if indicated, to determine whether the power
of that root is exactly equal to $n$.

int n_is_perfect_power(mp_limb_t *r, mp_limb_t x)

Returns maximum $k$, where $k > 1$, such that $x = y^k$, for some
positive integer $y$. $r$ contains the value of $y$. if no such $k$ exist,
$0$ is returned.

mp_limb_t n_rootrem(mp_limb_t* remainder, mp_limb_t n, mp_limb_t root)

This function uses the Newton iteration method to calculate the nth root of
Expand Down
154 changes: 154 additions & 0 deletions ulong_extras/is_perfect_power.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
/*=============================================================================
This file is part of FLINT.
FLINT is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
FLINT is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with FLINT; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Nitin Kumar
******************************************************************************/

#include <stdlib.h>
#include <flint.h>
#include <ulong_extras.h>


/*
check for the bit size of mp_limb_t, which is equal to FLINT_BITS,
q[i] - 1 congruent class modulo q[i] (about q[i] line no. 105)
are stored as follows:
qclass[i][j]th word contain FLINT_BITS Number of bit and
kth bit (from LSB To MSB) represents the,
(say y = FLINT_BITS * max(j - 1, 0) + k), y mod q[i] congruent
class of q[i] and it is set to 1 if y satisfy the property mentioned
in the line no. 106 otherwise it is set to 0.
for total q[i] congruent class modulo q[i],
q[i] / FLINT_BITS + (q[i] % FLINT_BITS) ? 0 : 1;
word are used, i.e. the size of qclass[i]th row.
*/

#if FLINT_BITS == 32
mp_limb_t qclass[18][25] = {{ 0x3 }, { 0x43 }, { 0x403 }, { 0x10021003 }
, { 0x400003 },{ 0x40800002 , 0x100001 },
{ 0x2 , 0x300c000 , 0x0 , 0x41 },
{ 0x82 , 0x20080 , 0x40000 , 0x2000 ,
0x1004000 , 0x41000001 },
{ 0x2 , 0x4001 }, { 0x2 , 0x4000001 },
{ 0x42 , 0x100010 , UWORD(0x80000000) , 0x0
, 0x0 , 0x0 , 0x1000000 , 0x0 , 0x80008 ,
0x420001 }, { 0x2 , 0x1000 , 0x0 , 0x200 ,
0x100001 }, { 0x2 , 0x0 , 0x40001 },
{ 0x2 , 0x0 , 0x20010000 , 0x0 , 0x0 ,
0x1001 }, { 0x2 , 0x3000 , 0x0 , 0x0 , 0x0 ,
0x0 , 0x0 , 0xc000 , 0x4000001 },
{ 0x2 , 0x0 , 0x0 , 0x401 }, { 0x2 , 0x0 ,
0x8000001 , 0x0 , 0x0 , 0x8000000 , 0x0 ,
0x18 , 0x0 , 0x0 , 0x0 , 0x0 , 0x0 , 0x0 ,
0x0 , 0x6 , 0x400 , 0x0 , 0x0 , 0x420 , 0x0
, 0x0 , 0x11 }, { 0x2 , 0x0 , 0x180000 ,
0x0 , 0x0 , 0x0 , 0x0 , 0x0 , 0x18000000 ,
0x0 , 0x0 , 0x4001 }
};
#else
mp_limb_t qclass[18][15] = { { 0x3 }, { 0x43 }, { 0x403 },
{ 0x10021003 } , { 0x400003 },
{ UWORD(0x10000040800003) },
{ UWORD(0x300c00000000002) ,
UWORD(0x4000000001) },
{ UWORD(0x2008000000082) ,
UWORD(0x200000040000) ,
UWORD(0x4100000001004001) },
{ UWORD(0x400000000003) },
{ UWORD(0x400000000000003) },
{ UWORD(0x10001000000042)
, UWORD(0x80000000) , 0x0 , 0x1000000 ,
UWORD(0x42000000080009) },
{ UWORD(0x100000000002) ,
UWORD(0x20000000000) , 0x100001 },{ 0x2 ,
0x40001 }, { 0x2 , 0x20010000 ,
UWORD(0x100000000001) }, {
UWORD(0x300000000002) , 0x0 ,
0x0 , UWORD(0xc00000000000) , 0x4000001 },
{ 0x2 , UWORD(0x40000000001) }, { 0x2 ,
UWORD(0x108000000) ,
UWORD(0x800000000000000) ,
UWORD(0x1800000000) , 0x0 , 0x0 , 0x0 ,
UWORD(0x600000000) , 0x400 ,
UWORD(0x42000000000) ,
0x0 , 0x11 }, { 0x2 , 0x180000 , 0x0 ,
0x0 , 0x18000000 , UWORD(0x400000000001) }
};
#endif


int n_is_perfect_power(mp_limb_t * r, mp_limb_t x)
{
/*
p[i], is prime less than 64.
q[i], is prime such that q[i] mod p[i] = 1.
if x = y ^ p[i] , for some y than,
x ^ ((q[i] - 1) / p[i]) mod q[i] <= 1.
see this following link for info:
("http://citeseerx.ist.psu.edu/viewdoc/
download?doi=10.1.1.108.458&rep=rep1&type=pdf")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

AKA
Eric Bach, Jonathan Sorenson,
Sieve algorithms for perfect power testing,
Algorithmica 9 (1993), 313–328. MR 94d:11103

z[i] = maximum k such that i = j ^ k for some positive integer j.
v[i] = minimum j such that i = j ^ k for some positive integer k.
*/

int p[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61};
int q[] = {3, 7, 11, 29, 23, 53, 103, 191, 47, 59,
311, 149, 83, 173, 283, 107, 709, 367};
int z[] = {1, 1, 1, 1, 2, 1, 1, 1, 3, 2,
1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1};
int v[] = {0, 1, 2, 3, 2, 5, 6, 7, 2, 3,
10, 11, 12, 13, 14, 15, 2, 17, 18, 19, 20};

int i, power = 1, pos, k;
mp_limb_t prev;

for (i = 0; i < 18; i++)
{
if (x <= 20)
{
power *= z[x];
x = v[x];
break;
}

/* k is congruent class of x mod q[i] */
k = x % q[i];
/* pos is the position of above class, in the table*/
pos = k / FLINT_BITS + (k % FLINT_BITS) ? 0 : 1;

if (pos == 0) pos = 1; /* in case q[i] divides x */
/* check if bit corresponding to above position is set */
if ((qclass[i][pos - 1] & (1 << (k % FLINT_BITS))))
{
prev = x;
x = n_rootrem(r, x, p[i]); /* test if x is p[i]th power */
if (*r == 0)
{
while (*r == 0) /* test also for subsequent power of p[i] */
{
power = power * p[i];
prev = x;
x = n_rootrem(r, x, p[i]);
}
x = prev;
}
else x = prev;
}
}

*r = x; /* store root in r */
return (power == 1 ? 0 : power); /* return power */
}
77 changes: 77 additions & 0 deletions ulong_extras/test/t-is_perfect_power.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,77 @@
/*=============================================================================
This file is part of FLINT.
FLINT is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
FLINT is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with FLINT; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Nitin Kumar
******************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <flint.h>
#include <ulong_extras.h>

int main()
{
int i, result, j, k;
mp_limb_t d, r, t;
FLINT_TEST_INIT(state);

flint_printf("n_is_perfect_power....");
fflush(stdout);

/* check for random number */
for (i = 0; i < 100000; i++)
{
d = n_randlimb(state);
result = n_is_perfect_power(&r, d);
if (result && n_pow(r, result) != d)
{
flint_printf("FAIL:\n");
flint_printf("%wu ** %wu != %wu\n", r, result, d);
abort();
}
}

/* check for possible power of random numbers, which are not
perfect powers.
*/

for (j = 0; j < 100000; j++)
{
d = n_randtest_not_zero(state);

if (d == 1) continue;

while (d < UWORD_MAX && n_is_perfect_power(&r, d) != 0) d += 1;

i = FLINT_BIT_COUNT(d);

for (k = 2; k <= FLINT_BITS / i; k++)
{
t = n_pow(d, k);
result = n_is_perfect_power(&r, t);
if (result != k)
{
flint_printf("FAIL:\n");
flint_printf("%wu ** %wu != %wu ( = %wu ** %wu)\n", r, result, t, d, k);
abort();
}
}
}

FLINT_TEST_CLEANUP(state);

flint_printf("PASS\n");
return 0;
}