gcd (f(n), g(n))>1 ?

[CRGreathouse]

sm88 and I have been working on a still-faster method; I'll let him post the details once he's pulled them together.

[science_man_88]

Originally Posted by CRGreathouse

sm88 and I have been working on a still-faster method; I'll let him post the details once he's pulled them together.

I think you have most of it ( you did most of the speeding up of my suggestion at least) the main part if I understand you correct is now putting a for loop around it to check each one or using apply.

[axn]

Originally Posted by CRGreathouse

(2n^2+n+9) - 2*(n^2+4n-5) = n^2-3n+14

You shouldn't have the multiplier 2. That will directly give 7n-19. Should be (2n^2+n+9) - (n^2+4n-5) = n^2-3n+14

[Dr Sardonicus]

The resultant approach suggested by CRGreathouse works for any polynomials of positive degree with integer coefficients, and allows a precise, detailed determination of gcd's. I am alas too lazy to work out all the details. However, I note the following:

1) If the two polynomials f1 = f1(x) and f2 = f2(x) have a gcd of positive degree, the resultant is 0. I suppose this possibility could be checked before proceeding further.

2) If the resultant is 1 or -1, f1(n) and f2(n) are always relatively prime.

3) For integers n, gcd(f1(n), f2(n)) always divides r = polresultant(f1,f2).

Note: If |r| is small, simply checking from 0 to |r| - 1 might be the way to go.

4) For each prime factor p of r, p divides gcd(f1(n), f2(n)) precisely when x - Mod(n,p) is a linear factor of gcd(Mod(1,p)*f1, Mod(1,p)*f2).

That is, the reductions of f1 and f2 (mod p) have to have a linear common factor.

5) If p^e is a prime-power factor of r, the common linear factor of f1 and f2 (mod p) might lead to a common linear factor modulo p^k for k up to e. I am not exactly sure how to do this. However, I note that M = factor(r) gives the prime factors of r, along with the exponents.

6) Solutions for the various prime-power moduli can be combined using the Chinese Reamainder Theorem.

Taking the example f1=3*x^2+5*x+4 and f2=5*x^2+6*x+13 provided by science_man_88, we find r = 648 = 2^3 * 3^4. For the prime factor p = 2, we find that 2 divides gcd(f1(n), f2(n)) when n == 1 (mod 2).

Also, 4 divides the gcd when x == 1 (mod 4) and 8 divides the gcd when x == 5 (mod 8).

For the prime factor 3, we find that 3 divides gcd(f1(n), f2(n)) when x == 1 (mod 3).

Also, 9 divides the gcd when x == 4 (mod 9), 27 when x == 22 (mod 27), and 81 when x == 49 (mod 81).

[CRGreathouse]

Originally Posted by Dr Sardonicus

4) For each prime factor p of r, p divides gcd(f1(n), f2(n)) precisely when x - Mod(n,p) is a linear factor of gcd(Mod(1,p)*f1, Mod(1,p)*f2).

That is, the reductions of f1 and f2 (mod p) have to have a linear common factor.

You can also check polrootsmod(f, p) and polrootsmod(g, p).

[CRGreathouse]

Originally Posted by axn

You shouldn't have the multiplier 2. That will directly give 7n-19. Should be (2n^2+n+9) - (n^2+4n-5) = n^2-3n+14

Yes, you can see I worked and re-worked a few times.

[Dr Sardonicus]

Originally Posted by CRGreathouse

You can also check polrootsmod(f, p) and polrootsmod(g, p).

Good! Meanwhile, back at the problem of extending solutions mod p to solutions mod p^k. I thought of using polrootspadic(), but alas this doesn't always work. The problem is, if a power of p divides the discriminant, factorpadic() won't make sense except for larger powers of p.

For example, the polynomial f = 5*x^2+6*x+13 has discriminant D = 36-260 = -224 = -7*2^5. Feeding it to Pari-GP, we find

? f = 5*x^2+6*x+13;

? subst(f,x,Mod(1,4))
%2 = Mod(0, 4)

? subst(f,x,Mod(5,8))
%3 = Mod(0, 8)

OK, so we have solutions to f == 0 (mod 8). However, factorpadic() doesn't fare too well, even to order 2^2:

? factorpadic(f,2,2)
%4 =
[(1 + O(2^2))*x^2 + (2 + O(2^2))*x + (1 + O(2^2)) 1]

[CRGreathouse]

Originally Posted by Dr Sardonicus

Meanwhile, back at the problem of extending solutions mod p to solutions mod p^k.

Of course this is not needed for the OP's problem, which only concerns finding where the gcd is > 1. But isn't this just Hensel lifting?

[science_man_88]

gcdpol(f,g)=forprime(n=2,100,if(setminus(Vec(polrootsmod(f,n,1)),Vec(polrootsmod(g,n,1)))!=[],print(n)))

failed on me so I sent it to CRG ( since he's the smart one).

[CRGreathouse]

commonRoots(f,g,p)=setintersect(polrootsmod(f,p)~, polrootsmod(g,p)~);
residues(f,g)=
{
  my(fac=factor(polresultant(f,g))[,1]);
  concat(apply(p->commonRoots(f,g,p), fac)~);
}
findResiduesInRange(v, start, end)=
{
  my(u=List());
  for(i=1,#v,
    forstep(n=(ceil(start)-lift(v[i]))\v[i].mod*v[i].mod+lift(v[i]), end, v[i].mod,
      listput(u,n)
    )
  );
  Set(u);
}
commonFactor(f,g,start,end)=
{
  findResiduesInRange(residues(f,g), start, end);
}
addhelp(commonFactor, "commonFactor(f, g, start, end): Given univariate integer polynomials f and g, return a vector containing all the integers start <= n <= end on which f(n) and g(n) have a common factor.");

\\ New example
P=32*x^2-17*x+4;Q=17*x^2+8*x+9;
v=commonFactor(P,Q,1,10^6); #v

[bhelmes]

A peaceful evening,

i tried to change the following gcd procedure from gmp 6.1.2 in order that the calculation of gcd stops and return 1 when a minimum limit is reached.

Does anybody understand where to add the if statement and
return the value equal 1

Perhaps something like that

mpz_gcd_min_break (mpz_ptr g, mpz_srcptr u, mpz_srcptr v, mpz_srcptr min)
{
mp_size_t minsize;
minsize = ABSIZ (min);

if (minsize > vsize)
{
SIZ (g) = 1;
PTR (g)[0] = mpn_gcd_1 (vp, vsize, up[0]);
return;
}
}

The original source code is :

/* mpz/gcd.c: Calculate the greatest common divisor of two integers.

Copyright 1991, 1993, 1994, 1996, 2000-2002, 2005, 2010 Free Software
Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

* the GNU Lesser General Public License as published by the Free
Software Foundation; either version 3 of the License, or (at your
option) any later version.

or

* 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.

or both in parallel, as here.

The GNU MP Library 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 copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library. If not,
see https://www.gnu.org/licenses/. */

#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"


void
mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
{
unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
mp_ptr tp;
mp_ptr up;
mp_size_t usize;
mp_ptr vp;
mp_size_t vsize;
mp_size_t gsize;
TMP_DECL;

up = PTR(u);
usize = ABSIZ (u);
vp = PTR(v);
vsize = ABSIZ (v);
/* GCD(0, V) == V. */
if (usize == 0)
{
SIZ (g) = vsize;
if (g == v)
return;
MPZ_REALLOC (g, vsize);
MPN_COPY (PTR (g), vp, vsize);
return;
}

/* GCD(U, 0) == U. */
if (vsize == 0)
{
SIZ (g) = usize;
if (g == u)
return;
MPZ_REALLOC (g, usize);
MPN_COPY (PTR (g), up, usize);
return;
}

if (usize == 1)
{
SIZ (g) = 1;
PTR (g)[0] = mpn_gcd_1 (vp, vsize, up[0]);
return;
}

if (vsize == 1)
{
SIZ(g) = 1;
PTR (g)[0] = mpn_gcd_1 (up, usize, vp[0]);
return;
}

TMP_MARK;

/* Eliminate low zero bits from U and V and move to temporary storage. */
while (*up == 0)
up++;
u_zero_limbs = up - PTR(u);
usize -= u_zero_limbs;
count_trailing_zeros (u_zero_bits, *up);
tp = up;
up = TMP_ALLOC_LIMBS (usize);
if (u_zero_bits != 0)
{
mpn_rshift (up, tp, usize, u_zero_bits);
usize -= up[usize - 1] == 0;
}
else
MPN_COPY (up, tp, usize);

while (*vp == 0)
vp++;
v_zero_limbs = vp - PTR (v);
vsize -= v_zero_limbs;
count_trailing_zeros (v_zero_bits, *vp);
tp = vp;
vp = TMP_ALLOC_LIMBS (vsize);
if (v_zero_bits != 0)
{
mpn_rshift (vp, tp, vsize, v_zero_bits);
vsize -= vp[vsize - 1] == 0;
}
else
MPN_COPY (vp, tp, vsize);

if (u_zero_limbs > v_zero_limbs)
{
g_zero_limbs = v_zero_limbs;
g_zero_bits = v_zero_bits;
}
else if (u_zero_limbs < v_zero_limbs)
{
g_zero_limbs = u_zero_limbs;
g_zero_bits = u_zero_bits;
}
else /* Equal. */
{
g_zero_limbs = u_zero_limbs;
g_zero_bits = MIN (u_zero_bits, v_zero_bits);
}

/* Call mpn_gcd. The 2nd argument must not have more bits than the 1st. */
vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
? mpn_gcd (vp, vp, vsize, up, usize)
: mpn_gcd (vp, up, usize, vp, vsize);

/* Here G <-- V << (g_zero_limbs*GMP_LIMB_BITS + g_zero_bits). */
gsize = vsize + g_zero_limbs;
if (g_zero_bits != 0)
{
mp_limb_t cy_limb;
gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
MPZ_REALLOC (g, gsize);
MPN_ZERO (PTR (g), g_zero_limbs);

tp = PTR(g) + g_zero_limbs;
cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
if (cy_limb != 0)
tp[vsize] = cy_limb;
}
else
{
MPZ_REALLOC (g, gsize);
MPN_ZERO (PTR (g), g_zero_limbs);
MPN_COPY (PTR (g) + g_zero_limbs, vp, vsize);
}

SIZ (g) = gsize;
TMP_FREE;
}

which use the function:
/* mpn/gcd.c: mpn_gcd for gcd of two odd integers.

Copyright 1991, 1993-1998, 2000-2005, 2008, 2010, 2012 Free Software
Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

* the GNU Lesser General Public License as published by the Free
Software Foundation; either version 3 of the License, or (at your
option) any later version.

or

* 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.

or both in parallel, as here.

The GNU MP Library 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 copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library. If not,
see https://www.gnu.org/licenses/. */

#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"

/* Uses the HGCD operation described in

N. Möller, On Schönhage's algorithm and subquadratic integer gcd
computation, Math. Comp. 77 (2008), 589-607.

to reduce inputs until they are of size below GCD_DC_THRESHOLD, and
then uses Lehmer's algorithm.
*/

/* Some reasonable choices are n / 2 (same as in hgcd), and p = (n +
* 2)/3, which gives a balanced multiplication in
* mpn_hgcd_matrix_adjust. However, p = 2 n/3 gives slightly better
* performance. The matrix-vector multiplication is then
* 4:1-unbalanced, with matrix elements of size n/6, and vector
* elements of size p = 2n/3. */

/* From analysis of the theoretical running time, it appears that when
* multiplication takes time O(n^alpha), p should be chosen so that
* the ratio of the time for the mpn_hgcd call, and the time for the
* multiplication in mpn_hgcd_matrix_adjust, is roughly 1/(alpha -
* 1). */
#ifdef TUNE_GCD_P
#define P_TABLE_SIZE 10000
mp_size_t p_table[P_TABLE_SIZE];
#define CHOOSE_P(n) ( (n) < P_TABLE_SIZE ? p_table[n] : 2*(n)/3)
#else
#define CHOOSE_P(n) (2*(n) / 3)
#endif

struct gcd_ctx
{
mp_ptr gp;
mp_size_t gn;
};

static void
gcd_hook (void *p, mp_srcptr gp, mp_size_t gn,
mp_srcptr qp, mp_size_t qn, int d)
{
struct gcd_ctx *ctx = (struct gcd_ctx *) p;
MPN_COPY (ctx->gp, gp, gn);
ctx->gn = gn;
}

#if GMP_NAIL_BITS > 0
/* Nail supports should be easy, replacing the sub_ddmmss with nails
* logic. */
#error Nails not supported.
#endif

/* Use binary algorithm to compute G <-- GCD (U, V) for usize, vsize == 2.
Both U and V must be odd. */
static inline mp_size_t
gcd_2 (mp_ptr gp, mp_srcptr up, mp_srcptr vp)
{
mp_limb_t u0, u1, v0, v1;
mp_size_t gn;

u0 = up[0];
u1 = up[1];
v0 = vp[0];
v1 = vp[1];

ASSERT (u0 & 1);
ASSERT (v0 & 1);

/* Check for u0 != v0 needed to ensure that argument to
* count_trailing_zeros is non-zero. */
while (u1 != v1 && u0 != v0)
{
unsigned long int r;
if (u1 > v1)
{
sub_ddmmss (u1, u0, u1, u0, v1, v0);
count_trailing_zeros (r, u0);
u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
u1 >>= r;
}
else /* u1 < v1. */
{
sub_ddmmss (v1, v0, v1, v0, u1, u0);
count_trailing_zeros (r, v0);
v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
v1 >>= r;
}
}

gp[0] = u0, gp[1] = u1, gn = 1 + (u1 != 0);

/* If U == V == GCD, done. Otherwise, compute GCD (V, |U - V|). */
if (u1 == v1 && u0 == v0)
return gn;

v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
gp[0] = mpn_gcd_1 (gp, gn, v0);

return 1;
}

mp_size_t
mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
{
mp_size_t talloc;
mp_size_t scratch;
mp_size_t matrix_scratch;

struct gcd_ctx ctx;
mp_ptr tp;
TMP_DECL;

ASSERT (usize >= n);
ASSERT (n > 0);
ASSERT (vp[n-1] > 0);

/* FIXME: Check for small sizes first, before setting up temporary
storage etc. */
talloc = MPN_GCD_SUBDIV_STEP_ITCH(n);

/* For initial division */
scratch = usize - n + 1;
if (scratch > talloc)
talloc = scratch;

#if TUNE_GCD_P
if (CHOOSE_P (n) > 0)
#else
if (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD))
#endif
{
mp_size_t hgcd_scratch;
mp_size_t update_scratch;
mp_size_t p = CHOOSE_P (n);
mp_size_t scratch;
#if TUNE_GCD_P
/* Worst case, since we don't guarantee that n - CHOOSE_P(n)
is increasing */
matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n);
hgcd_scratch = mpn_hgcd_itch (n);
update_scratch = 2*(n - 1);
#else
matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
hgcd_scratch = mpn_hgcd_itch (n - p);
update_scratch = p + n - 1;
#endif
scratch = matrix_scratch + MAX(hgcd_scratch, update_scratch);
if (scratch > talloc)
talloc = scratch;
}

TMP_MARK;
tp = TMP_ALLOC_LIMBS(talloc);

if (usize > n)
{
mpn_tdiv_qr (tp, up, 0, up, usize, vp, n);

if (mpn_zero_p (up, n))
{
MPN_COPY (gp, vp, n);
ctx.gn = n;
goto done;
}
}

ctx.gp = gp;

#if TUNE_GCD_P
while (CHOOSE_P (n) > 0)
#else
while (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD))
#endif
{
struct hgcd_matrix M;
mp_size_t p = CHOOSE_P (n);
mp_size_t matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
mp_size_t nn;
mpn_hgcd_matrix_init (&M, n - p, tp);
nn = mpn_hgcd (up + p, vp + p, n - p, &M, tp + matrix_scratch);
if (nn > 0)
{
ASSERT (M.n <= (n - p - 1)/2);
ASSERT (M.n + p <= (p + n - 1) / 2);
/* Temporary storage 2 (p + M->n) <= p + n - 1. */
n = mpn_hgcd_matrix_adjust (&M, p + nn, up, vp, p, tp + matrix_scratch);
}
else
{
/* Temporary storage n */
n = mpn_gcd_subdiv_step (up, vp, n, 0, gcd_hook, &ctx, tp);
if (n == 0)
goto done;
}
}

while (n > 2)
{
struct hgcd_matrix1 M;
mp_limb_t uh, ul, vh, vl;
mp_limb_t mask;

mask = up[n-1] | vp[n-1];
ASSERT (mask > 0);

if (mask & GMP_NUMB_HIGHBIT)
{
uh = up[n-1]; ul = up[n-2];
vh = vp[n-1]; vl = vp[n-2];
}
else
{
int shift;

count_leading_zeros (shift, mask);
uh = MPN_EXTRACT_NUMB (shift, up[n-1], up[n-2]);
ul = MPN_EXTRACT_NUMB (shift, up[n-2], up[n-3]);
vh = MPN_EXTRACT_NUMB (shift, vp[n-1], vp[n-2]);
vl = MPN_EXTRACT_NUMB (shift, vp[n-2], vp[n-3]);
}

/* Try an mpn_hgcd2 step */
if (mpn_hgcd2 (uh, ul, vh, vl, &M))
{
n = mpn_matrix22_mul1_inverse_vector (&M, tp, up, vp, n);
MP_PTR_SWAP (up, tp);
}
else
{
/* mpn_hgcd2 has failed. Then either one of a or b is very
small, or the difference is very small. Perform one
subtraction followed by one division. */

/* Temporary storage n */
n = mpn_gcd_subdiv_step (up, vp, n, 0, &gcd_hook, &ctx, tp);
if (n == 0)
goto done;
}
}

ASSERT(up[n-1] | vp[n-1]);

if (n == 1)
{
*gp = mpn_gcd_1(up, 1, vp[0]);
ctx.gn = 1;
goto done;
}

/* Due to the calling convention for mpn_gcd, at most one can be
even. */

if (! (up[0] & 1))
MP_PTR_SWAP (up, vp);

ASSERT (up[0] & 1);

if (vp[0] == 0)
{
*gp = mpn_gcd_1 (up, 2, vp[1]);
ctx.gn = 1;
goto done;
}
else if (! (vp[0] & 1))
{
int r;
count_trailing_zeros (r, vp[0]);
vp[0] = ((vp[1] << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (vp[0] >> r);
vp[1] >>= r;
}

ctx.gn = gcd_2(gp, up, vp);

done:
TMP_FREE;
return ctx.gn;
}

Thanks in advance for every helpful advice
Greetings from the sinus and cosinus
Bernhard