Skip to content

Commit 374e7e9

Browse files
dimpasedimpase
committed
fix integer factoring regression from #35219 
1 parent 517657d commit 374e7e9

File tree

2 files changed

+194
-0
lines changed

2 files changed

+194
-0
lines changed

‎build/pkgs/pari/patches/Keep-product-of-prime-numbers-to-use-in-Z_factor_limit.patch

Lines changed: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
1+
From ea6992e38d651643982fd88777d5d5a186b0f54e Mon Sep 17 00:00:00 2001
2+
From: Bill Allombert <[email protected]>
3+
Date: Wed, 16 Nov 2022 13:33:21 +0100
4+
Subject: [PATCH] Keep product of prime numbers to use in Z_factor_limit
5+
6+
---
7+
src/basemath/ifactor1.c | 67 ++++++++++++++++++++++++++++++++---------
8+
src/headers/paridecl.h | 2 ++
9+
src/language/forprime.c | 26 +++++++++++++++-
10+
3 files changed, 80 insertions(+), 15 deletions(-)
11+
12+
diff --git a/src/basemath/ifactor1.c b/src/basemath/ifactor1.c
13+
index 4478fc416..5dd04c966 100644
14+
--- a/src/basemath/ifactor1.c
15+
+++ b/src/basemath/ifactor1.c
16+
@@ -3652,29 +3652,68 @@ ifactor_sign(GEN n, ulong all, long hint, long sn, GEN *pU)
17+
av = avma; affii(shifti(n,-i), n); set_avma(av);
18+
}
19+
if (is_pm1(n)) return aux_end(M,n,nb);
20+
- /* trial division */
21+
maxp = maxprime();
22+
- av = avma; u_forprime_init(&T, 3, minss(lim, maxp)); av2 = avma;
23+
- /* first pass: known to fit in private prime table */
24+
- while ((p = u_forprime_next_fast(&T)))
25+
+ if (lim >= maxprimelim()>>2)
26+
{
27+
- pari_sp av3 = avma;
28+
- int stop;
29+
- long k = Z_lvalrem_stop(&n, p, &stop);
30+
- if (k)
31+
+ GEN nr, NR;
32+
+ /* fast trial division */
33+
+ av = avma;
34+
+ NR = nr = gcdii(prodprimes(),n);
35+
+ u_forprime_init(&T, 3, minss(lim, maxp)); av2 = avma;
36+
+ /* first pass: known to fit in private prime table */
37+
+ while ((p = u_forprime_next_fast(&T)))
38+
{
39+
- affii(n, N); n = N; set_avma(av3);
40+
- STOREu(&nb, p, k);
41+
+ pari_sp av3 = avma;
42+
+ int stop;
43+
+ long k;
44+
+ if (!dvdiu(nr, p)) continue;
45+
+ nr = diviuexact(nr, p);
46+
+ affii(nr, NR); nr = NR; set_avma(av3);
47+
+ k = Z_lvalrem_stop(&n, p, &stop);
48+
+ if (k)
49+
+ {
50+
+ affii(n, N); n = N; set_avma(av3);
51+
+ STOREu(&nb, p, k);
52+
+ }
53+
+ if (is_pm1(n))
54+
+ {
55+
+ stackdummy(av, av2);
56+
+ return aux_end(M,n,nb);
57+
+ }
58+
+ if (is_pm1(nr)) break;
59+
}
60+
- if (p == 16381 && bit_accuracy(lgefint(n)) < 2048)
61+
- { stop = ifac_isprime(n); nb0 = nb; }
62+
- if (stop)
63+
+ if (ifac_isprime(n))
64+
{
65+
- if (!is_pm1(n)) STOREi(&nb, n, 1);
66+
+ STOREi(&nb, n, 1);
67+
stackdummy(av, av2);
68+
return aux_end(M,n,nb);
69+
}
70+
}
71+
+ else
72+
+ {
73+
+ /* trial division */
74+
+ av = avma; u_forprime_init(&T, 3, minss(lim, maxp)); av2 = avma;
75+
+ /* first pass: known to fit in private prime table */
76+
+ while ((p = u_forprime_next_fast(&T)))
77+
+ {
78+
+ pari_sp av3 = avma;
79+
+ int stop;
80+
+ long k = Z_lvalrem_stop(&n, p, &stop);
81+
+ if (k)
82+
+ {
83+
+ affii(n, N); n = N; set_avma(av3);
84+
+ STOREu(&nb, p, k);
85+
+ }
86+
+ if (p == 16381 && bit_accuracy(lgefint(n)) < 2048)
87+
+ { stop = ifac_isprime(n); nb0 = nb; }
88+
+ if (stop)
89+
+ {
90+
+ if (!is_pm1(n)) STOREi(&nb, n, 1);
91+
+ stackdummy(av, av2);
92+
+ return aux_end(M,n,nb);
93+
+ }
94+
+ }
95+
+ }
96+
stackdummy(av, av2);
97+
if (lim > maxp)
98+
{ /* second pass, usually empty: outside private prime table */
99+
diff --git a/src/headers/paridecl.h b/src/headers/paridecl.h
100+
index be375c67d..1d5041e7f 100644
101+
--- a/src/headers/paridecl.h
102+
+++ b/src/headers/paridecl.h
103+
@@ -3565,7 +3565,9 @@ ulong init_primepointer_lt(ulong a, byteptr *pd);
104+
ulong maxprime(void);
105+
ulong maxprimeN(void);
106+
void maxprime_check(ulong c);
107+
+ulong maxprimelim(void);
108+
void pari_init_primes(ulong maxprime);
109+
+GEN prodprimes(void);
110+
ulong u_forprime_next(forprime_t *T);
111+
int u_forprime_init(forprime_t *T, ulong a, ulong b);
112+
void u_forprime_restrict(forprime_t *T, ulong c);
113+
diff --git a/src/language/forprime.c b/src/language/forprime.c
114+
index 22be44726..1059900df 100644
115+
--- a/src/language/forprime.c
116+
+++ b/src/language/forprime.c
117+
@@ -22,7 +22,9 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
118+
/**********************************************************************/
119+
120+
static ulong _maxprime = 0;
121+
+static ulong _maxprimelim = 0;
122+
static ulong diffptrlen;
123+
+static GEN _prodprimes;
124+
125+
/* Building/Rebuilding the diffptr table. The actual work is done by the
126+
* following two subroutines; the user entry point is the function
127+
@@ -352,6 +354,18 @@ sieve_chunk(byteptr known_primes, ulong s, byteptr data, ulong n)
128+
}
129+
}
130+
131+
+static void
132+
+set_prodprimes(void)
133+
+{
134+
+ pari_sp av = avma;
135+
+ GEN p = zv_prod_Z(primes_zv(diffptrlen-1));
136+
+ long i, l = lgefint(p);
137+
+ _prodprimes = (GEN) pari_malloc(l*sizeof(long));
138+
+ for (i = 0; i < l; i++)
139+
+ _prodprimes[i] = p[i];
140+
+ set_avma(av);
141+
+}
142+
+
143+
/* assume maxnum <= 436273289 < 2^29 */
144+
static void
145+
initprimes0(ulong maxnum, long *lenp, ulong *lastp, byteptr p1)
146+
@@ -426,7 +440,11 @@ initprimes0(ulong maxnum, long *lenp, ulong *lastp, byteptr p1)
147+
ulong
148+
maxprime(void) { return diffptr ? _maxprime : 0; }
149+
ulong
150+
+maxprimelim(void) { return diffptr ? _maxprimelim : 0; }
151+
+ulong
152+
maxprimeN(void) { return diffptr ? diffptrlen-1: 0; }
153+
+GEN
154+
+prodprimes(void) { return diffptr ? _prodprimes: NULL; }
155+
156+
void
157+
maxprime_check(ulong c) { if (_maxprime < c) pari_err_MAXPRIME(c); }
158+
@@ -445,6 +463,7 @@ initprimes(ulong maxnum, long *lenp, ulong *lastp)
159+
maxnum = 436273289;
160+
t = (byteptr)pari_malloc((size_t) (1.09 * maxnum/log((double)maxnum)) + 146);
161+
initprimes0(maxnum, lenp, lastp, t);
162+
+ _maxprimelim = maxnum;
163+
return (byteptr)pari_realloc(t, *lenp);
164+
}
165+
166+
@@ -457,6 +476,7 @@ initprimetable(ulong maxnum)
167+
diffptrlen = minss(diffptrlen, len);
168+
_maxprime = minss(_maxprime,last); /*Protect against ^C*/
169+
diffptr = p; diffptrlen = len; _maxprime = last;
170+
+ set_prodprimes();
171+
if (old) free(old);
172+
}
173+
174+
@@ -770,7 +790,11 @@ pari_init_primes(ulong maxprime)
175+
void
176+
pari_close_primes(void)
177+
{
178+
- pari_free(diffptr);
179+
+ if (diffptr)
180+
+ {
181+
+ pari_free(diffptr);
182+
+ pari_free(_prodprimes);
183+
+ }
184+
pari_free(pari_sieve_modular.sieve);
185+
}
186+
187+
--
188+
2.30.2
189+

‎src/sage/arith/misc.py

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2688,6 +2688,11 @@ def factor(n, proof=None, int_=False, algorithm='pari', verbose=0, **kwds):
26882688
Traceback (most recent call last):
26892689
...
26902690
TypeError: unable to factor 'xyz'
2691+
2692+
test that :issue:`35219` is fixed::
2693+
2694+
sage: len(factor(2^2203-1,proof=false))
2695+
1
26912696
"""
26922697
try:
26932698
m = n.factor

0 commit comments

Comments
 (0)