qsieve does not return for some numbers

[kruoli]

qsieve_factor will not return for some numbers, e.g. 500000000000000000000000000000000000000017711.

Minimal example:

#include <stdio.h>
#include <stdlib.h>
#include <flint/fmpz.h>
#include <flint/fmpz_factor.h>
#include <flint/qsieve.h>

int main(int argc, char *argv[])
{
	fmpz_t x;
	fmpz_init(x);
	fmpz_set_str(x, "500000000000000000000000000000000000000017711", 10);
	fmpz_factor_t f;
	fmpz_factor_init(f);
	qsieve_factor(f, x);
	fmpz_print(f->p);
	unsigned long exp = f->exp[0];
	if (exp != 1)
	{
		printf("^%lu", exp);
	}
	for (long i = 1; i < f->num; i++)
	{
		printf("*");
		fmpz_print(f->p + i);
		exp = f->exp[i];
		if (exp != 1)
		{
			printf("^%lu", exp);
		}
	}
	fmpz_factor_clear(f);
	fmpz_clear(x);
	return EXIT_SUCCESS;
}

The example was compiled with gcc -o example -Og -g example.c -Wall -pedantic -lflint.

I would expect qsieve to return at some point.

System:

• OS: Debian
• CPU: tested on 5950X, some 11th-gen Intel CPU, 3800X
• version of FLINT: in all cases after the last commit to qsieve (2 months ago)
• options pushed to configure: --use-avx2, on the 11th-gen Intel also --use-avx512

[albinahlback]

Do you have more examples of where it fails to return in a reasonable amount of time?

I've not looked into it that much, but fmpz_factor also gets stuck (which it shouldn't). Perhaps some threshold should be changed.

[kruoli]

I did some tests from 25 digits to 65 digits in steps of 5 digits. It started occuring at 45 digits (when testing 1e5 numbers each). Is a C45 reasonable in your opinion? If so, I can surely find some more.

Edit: I misunderstood. I will generate some more examples.

[fredrik-johansson]

FWIW, this example succeeds if I toggle the alternative qsieve_tune table in qsieve.h.

Specifically, it works if one changes the 2 * 65536 in

{150, 100,  1500, 10,   2 *  65536, 70}, /* 46 digit */

to various smaller values, eg 2 * 30000.

I guess we could just do a big tuneup for qsieve where we test the new parameters very extensively for termination. Though it would be nice to understand why there is an issue with the larger sieve size too...

[kruoli]

Here is a list of composites that will get stuck in fmpz_factor (some of them with small factors; they will eventually get stuck in the QS step):

1000000000000000000000000000000000000000035422
1000000000000000000000000000000000000000042223
1000000000000000000000000000000000000000062878
1000000000000000000000000000000000000000078147
1000000000000000000000000000000000000000148937
1000000000000000000000000000000000000000169537
1000000000000000000000000000000000000000170202
1000000000000000000000000000000000000000186843
1000000000000000000000000000000000000000221741
1000000000000000000000000000000000000000252886
1000000000000000000000000000000000000000253483
1000000000000000000000000000000000000000280147
1000000000000000000000000000000000000000287517
1000000000000000000000000000000000000000322347
1000000000000000000000000000000000000000356278
1000000000000000000000000000000000000000368621
1000000000000000000000000000000000000000414187
1000000000000000000000000000000000000000415947
1000000000000000000000000000000000000000420217
1000000000000000000000000000000000000000449343
1000000000000000000000000000000000000000484747
1000000000000000000000000000000000000000489307
1000000000000000000000000000000000000000496267

When enabling QS_DEBUG and comparing e.g. 500000000000000000000000000000000000000021111 to 500000000000000000000000000000000000000017711, the latter runs out of A's really quickly:

<snip>

Initializing Relations and Linear Algebra

Polynomial Initialisation and Sieving
second prime index = 1500
s = 6
high = 266
low = 92
span = 174
First A_ind = (93, 94, 95, 97, 98, 100)
number of factors in hypercube = 6
factor base size = 1500 max prime = 27647
possible candidate in range [ 92, 266 ]
optimal value of hypercube = 1276641886190911175
lower bound       = 638320943095455587
upper bound       = 2553283772381822350
initial hypercube = 2540818837309778947
0 relations
full relations = 60, num cycles = 1, ks_primes = 100, extra rels = 64, poly_count = 32, num_primes = 1500
A_ind = (93, 94, 95, 99, 101, 92)
100 relations
full relations = 132, num cycles = 3, ks_primes = 100, extra rels = 64, poly_count = 64, num_primes = 1500
A_ind = (93, 94, 95, 102, 103, 92)
full relations = 194, num cycles = 4, ks_primes = 100, extra rels = 64, poly_count = 96, num_primes = 1500
A_ind = (93, 94, 98, 99, 101, 92)
200 relations
full relations = 261, num cycles = 7, ks_primes = 100, extra rels = 64, poly_count = 128, num_primes = 1500
--> compute_poly_data.c:443 sets 0 as return value for qsieve_next_A <--
Increasing factor base.

<snip>

[fredrik-johansson]

Can you generate these at various other bit sizes as well?

[fredrik-johansson]

Pinging @wbhart in case he has an idea how to fix this as robustly as possible.

[wbhart]

Fiddling with the tuning values is probably the only option here. Unfortunately at the small sizes, the sieve is a bit flakey. There's not a heap can be done about it.

I tuned it as best I could, but it looks like I still didn't get it right.

There's just not enough polynomials of the right size to factor the number comfortably, so there's always some tradeoff.

[wbhart]

I'm a bit surprised that making the sieve size smaller actually helps. That shows you just how marginal the behaviour is.

[kruoli]

Can you generate these at various other bit sizes as well?

A collection of composites in different sizes:

126582278481012658227848101265822806309
322580645161290322580645161290322654733
555555555555555555555555555555555693427
561482313307130825379000561482313307499
720461095100864553314121037463976947183
1628664495114006514657980456026058632341
1650516364044491319109183307997907145251
1666666666666666666666666666666666845161
1666666666666666666666666666666666979801
1672240802675585284280936454849498335537
95662720772928579912788805098723834088195937
97349439977949650929026562827088794528412033
111111111111111111111111111111111111111135749
277315585135884636716583471991125901275651731
333333333333333333333333333333333333333382979
500000000000000000000000000000000000000017711
500000000000000000000000000000000000000031439
500000000000000000000000000000000000000085101
500000000000000000000000000000000000000126443
500000000000000000000000000000000000000178139
571755288736420811892510005717552887364208149
633713561470215462610899873257287705956907549
970931288163667945383173178217350736305742379
1000000000000000000000000000000000000000042223
1000000000000000000000000000000000000000078147
1000000000000000000000000000000000000000169537
1000000000000000000000000000000000000000186843
1000000000000000000000000000000000000000253483
1000000000000000000000000000000000000000280147
1000000000000000000000000000000000000000287517
1000000000000000000000000000000000000000322347
1000000000000000000000000000000000000000414187
1000000000000000000000000000000000000000415947
1000000000000000000000000000000000000000420217
1000000000000000000000000000000000000000449343
1000000000000000000000000000000000000000484747
1000000000000000000000000000000000000000489307
1000000000000000000000000000000000000000496267
1030927835051546391752577319587628865979381499
1441266465979343215563261116966752565379363587
1537794372287715175877542358545984665114519547
8312897460409825844798204414148551477617523587847
12195121951219512195121951219512195121951219512199
23809523809523809523809523809523809523809523811539

[wbhart]

Hmm, I don't have an explanation for the failure at larger sizes. There should not be pressure on the number of available polynomials there.

I have to wonder if something hasn't become corrupted at some point.

I wish I had time to dig into it myself, but I absolutely don't.

[fredrik-johansson]

We should test these with some ancient versions of Flint to see if there has been a regression.

[wbhart]

Perhaps turn on some of the debugging code and see if you can narrow down why it is running out of polynomials, if that is what is going on. If you start with the larger numbers, that should be a clearer case, as there really ought to be enough.

From memory, from 40 digits on, you shouldn't run into much of a problem.

This whole implementation was meant to solve this problem, and in my testing it worked pretty well. But apparently I didn't test well enough.

I mean, the quadratic sieve can fail, but it shouldn't do so much, if at all, in practice.

[kruoli]

Something I saw in qsieve_next_A: First, some polynomials are found when curr_subset[0] == 0. After it gets increased, it does not find anything. This again repeats after each factor base increment. Because of that, the factor base is incremented over and over.

[kruoli]

It is possible that there gets something corrupted. When executing the example in OP with 1000000000000000000000000000000000000000420217, within I minute I got:

munmap_chunk(): invalid pointer

Program received signal SIGABRT, Aborted.
__pthread_kill_implementation (threadid=<optimized out>, signo=signo@entry=6, no_tid=no_tid@entry=0) at ./nptl/pthread_kill.c:44
44      ./nptl/pthread_kill.c: Datei oder Verzeichnis nicht gefunden.
(gdb) bt
#0  __pthread_kill_implementation (threadid=<optimized out>, signo=signo@entry=6, no_tid=no_tid@entry=0) at ./nptl/pthread_kill.c:44
#1  0x00007ffff70a9f1f in __pthread_kill_internal (signo=6, threadid=<optimized out>) at ./nptl/pthread_kill.c:78
#2  0x00007ffff705afb2 in __GI_raise (sig=sig@entry=6) at ../sysdeps/posix/raise.c:26
#3  0x00007ffff7045472 in __GI_abort () at ./stdlib/abort.c:79
#4  0x00007ffff709e430 in __libc_message (action=action@entry=do_abort, fmt=fmt@entry=0x7ffff71b8459 "%s\n") at ../sysdeps/posix/libc_fatal.c:155
#5  0x00007ffff70b383a in malloc_printerr (str=str@entry=0x7ffff71babb8 "munmap_chunk(): invalid pointer") at ./malloc/malloc.c:5660
#6  0x00007ffff70b39fc in munmap_chunk (p=<optimized out>) at ./malloc/malloc.c:3054
#7  0x00007ffff70b7f68 in __GI___libc_free (mem=<optimized out>) at ./malloc/malloc.c:3375
#8  0x00007ffff7917bc8 in qsieve_poly_clear () from /usr/local/lib/libflint.so.20
#9  0x00007ffff7915104 in qsieve_factor () from /usr/local/lib/libflint.so.20
#10 0x0000555555555203 in main (argc=<optimized out>, argv=<optimized out>) at example.c:14

[rburing]

Running the above 1000000000000000000000000000000000000000420217 example with

diff --git a/src/qsieve.h b/src/qsieve.h
index 68ab9cd52..680cc32be 100644
--- a/src/qsieve.h
+++ b/src/qsieve.h
@@ -27,7 +27,7 @@ extern "C" {
 /* Windows systems may define `small` macro, which leads to conflicts */
 #undef small
 
-#define QS_DEBUG 0 /* level of debug information printed, 0 = none */
+#define QS_DEBUG 64 + 128 /* level of debug information printed, 0 = none */
 
 #define BITS_ADJUST 25 /* add to sieve entries to compensate approximations */

on top of b36927c (current main branch) yields the debug output

start
factoring 1000000000000000000000000000000000000000420217 of 150 bits

Knuth_Schroeppel
Using Knuth-Schroeppel multiplier 3
kn bits = 152

Compute factor-base

Initializing Relations and Linear Algebra

Polynomial Initialisation and Sieving
second prime index = 1500
s = 6
high = 276
low = 87
span = 189
Increasing factor base.

factor base increment
munmap_chunk(): invalid pointer

Program received signal SIGABRT, Aborted.
__pthread_kill_implementation (no_tid=0, signo=6, threadid=<optimized out>) at ./nptl/pthread_kill.c:44
warning: 44     ./nptl/pthread_kill.c: No such file or directory
(gdb) bt
#0  __pthread_kill_implementation (no_tid=0, signo=6, threadid=<optimized out>) at ./nptl/pthread_kill.c:44
#1  __pthread_kill_internal (signo=6, threadid=<optimized out>) at ./nptl/pthread_kill.c:78
#2  __GI___pthread_kill (threadid=<optimized out>, signo=signo@entry=6) at ./nptl/pthread_kill.c:89
#3  0x00007ffff6a4527e in __GI_raise (sig=sig@entry=6) at ../sysdeps/posix/raise.c:26
#4  0x00007ffff6a288ff in __GI_abort () at ./stdlib/abort.c:79
#5  0x00007ffff6a297b6 in __libc_message_impl (fmt=fmt@entry=0x7ffff6bce8d7 "%s\n") at ../sysdeps/posix/libc_fatal.c:134
#6  0x00007ffff6aa8ff5 in malloc_printerr (str=str@entry=0x7ffff6bd1520 "munmap_chunk(): invalid pointer") at ./malloc/malloc.c:5772
#7  0x00007ffff6aa947c in munmap_chunk (p=<optimized out>) at ./malloc/malloc.c:3040
#8  0x00007ffff6aaddfa in __GI___libc_free (mem=0x5555555b78f0) at ./malloc/malloc.c:3388
#9  0x00007ffff6f4ea0a in _flint_free (ptr=0x5555555b78f0) at /home/rburing/src/flint/src/generic_files/memory_manager.c:144
#10 0x00007ffff6f4ea2d in flint_free (ptr=0x5555555b78f0) at /home/rburing/src/flint/src/generic_files/memory_manager.c:150
#11 0x00007ffff7734a59 in qsieve_poly_clear (qs_inf=0x7fffffffd300) at /home/rburing/src/flint/src/qsieve/poly.c:91
#12 0x00007ffff7730e9f in qsieve_factor (factors=0x7fffffffd560, n=0x7fffffffd558) at /home/rburing/src/flint/src/qsieve/factor.c:429
#13 0x00005555555552f6 in main (argc=1, argv=0x7fffffffd6b8) at example.c:14

[kruoli]

Regarding 1000000000000000000000000000000000000000420217, I found out that s will equal h at some point in compute_poly_data.c:273, which leads to accessing index -1 of an array. More values:

(gdb) p s
$1 = 6
(gdb) p h
$2 = 6
(gdb) p span
$3 = 189
(gdb) p m
$4 = 136
(gdb) p h
$5 = 6
(gdb) p i
$6 = 0
(gdb) p j
$7 = <optimized out>