1 /* factor.c -- factorization of very smooth numbers.
2 *
3 * Copyright (C) 2012, 2013, 2021 INRIA
4 *
5 * This file is part of CMH.
6 *
7 * CMH is free software; you can redistribute it and/or modify it under
8 * the terms of the GNU General Public License as published by the
9 * Free Software Foundation; either version 3 of the License, or (at your
10 * option) any later version.
11 *
12 * CMH is distributed in the hope that it will be useful, but WITHOUT ANY
13 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for
15 * more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see http://www.gnu.org/licenses/ .
19 */
20
21 #include <stdlib.h>
22 #include <stdio.h>
23 #include <limits.h>
24 #include <string.h>
25 #include <pari/pari.h>
26 #include "parsepari.h"
27 #include "factor.h"
28
factor_matrix_init(int size)29 static factor_matrix factor_matrix_init (int size)
30 /* Allocate a matrix f to hold size factors, where f [i, 0] is the i-th
31 prime factor and f [i, 1] its multiplicity; the last entry is
32 indicated by a prime factor 0 */
33 {
34 factor_matrix f;
35
36 f = (factor_matrix) malloc ((size + 1) * 2 * sizeof (unsigned long));
37
38 return f;
39 }
40
41
factor_matrix_dup(factor_matrix f)42 factor_matrix factor_matrix_dup(factor_matrix f)
43 {
44 int size;
45 for (size = 0; f [size][0] ; size++);
46 factor_matrix g = factor_matrix_init (size);
47 memcpy(g, f, (size + 1) * 2 * sizeof (unsigned long));
48 return g;
49 }
50
51
factor_matrix_free(factor_matrix f)52 void factor_matrix_free (factor_matrix f)
53 {
54 free (f);
55 }
56
57
factor_trialdiv(mpz_srcptr n,unsigned long int bound)58 factor_matrix factor_trialdiv (mpz_srcptr n, unsigned long int bound)
59 /* factors n by trial division with numbers <= bound, or indefinitely
60 if bound is 0.
61 The return value is a newly allocated factor list if the number is
62 smooth, NULL otherwise. */
63 {
64 pari_sp av;
65 factor_matrix f;
66 GEN n_gen, fact;
67 mpz_t tmp;
68 unsigned long int last;
69 int i, m, ok;
70
71 if (mpz_cmp_ui (n, 1) == 0) {
72 f = factor_matrix_init (0);
73 f [0][0] = 0;
74 return f;
75 }
76
77 if (bound == 0)
78 bound = ULONG_MAX;
79
80 pari_init (500000, 0);
81 av = avma;
82
83 n_gen = mpz_get_GEN (n);
84 fact = boundfact (n_gen, bound);
85 m = lg (gel (fact, 1)) - 1;
86 last = itou_or_0 (gel (gel (fact, 1), m));
87 if (last > bound || last == 0)
88 ok = 0;
89 else {
90 ok = 1;
91 mpz_init (tmp);
92 f = factor_matrix_init (m);
93 for (i = 0; i < m; i++) {
94 mpz_set_GEN (tmp, gel (gel (fact, 1), i+1));
95 f [i][0] = mpz_get_ui (tmp);
96 mpz_set_GEN (tmp, gel (gel (fact, 2), i+1));
97 f [i][1] = mpz_get_ui (tmp);
98 }
99 f [m][0] = 0;
100 mpz_clear (tmp);
101 }
102
103 avma = av;
104 pari_close ();
105
106 return (ok ? f : NULL);
107 }
108
109
factor_matrix_asprintf(char ** pres,factor_matrix f)110 void factor_matrix_asprintf (char ** pres, factor_matrix f)
111 {
112 size_t alloc = 10;
113 size_t size = 0;
114 char * res = malloc (alloc);
115 int d, i;
116
117 #define PUSH(fmt__, arg__) do { \
118 for( ; (d = snprintf(res + size, alloc-size, fmt__, arg__)) >= (int) (alloc-size) ; ) { \
119 alloc = size+d+1 + alloc / 4; \
120 res = realloc(res, alloc); \
121 } \
122 size+=d; \
123 } while (0)
124
125 if (f[0][0] == 0) {
126 snprintf(res, alloc, "1");
127 *pres = res;
128 return;
129 }
130 for (i = 0; f [i][0] != 0; i++) {
131 if (i > 0)
132 PUSH ("%c", '*');
133 PUSH ("%lu", f [i][0]);
134 if (f [i][1] > 1)
135 PUSH ("^%lu", f [i][1]);
136 }
137 #undef PUSH
138
139 *pres = res;
140 }
141
142
factor_matrix_fprintf(FILE * file,factor_matrix f)143 int factor_matrix_fprintf (FILE *file, factor_matrix f)
144 {
145 int ok, i, size;
146
147 for (size = 0; f [size][0] != 0; size++);
148 ok = fprintf (file, "%i", size);
149 for (i = 0; i < size; i++)
150 ok += fprintf (file, " %lu %lu", f [i][0], f [i][1]);
151
152 return (ok > 0);
153 }
154
155
factor_matrix_fold_inner(mpz_ptr z,factor_matrix f,int i0,int i1)156 static void factor_matrix_fold_inner(mpz_ptr z, factor_matrix f, int i0, int i1)
157 {
158 if (i1 - i0 == 1) {
159 mpz_ui_pow_ui(z, f[i0][0], f[i0][1]);
160 } else {
161 int j = (i0 + i1) / 2;
162 mpz_t l,r;
163 mpz_init(l);
164 mpz_init(r);
165 factor_matrix_fold_inner(l, f, i0, j);
166 factor_matrix_fold_inner(r, f, j, i1);
167 mpz_mul(z, l, r);
168 mpz_clear(l);
169 mpz_clear(r);
170 }
171 }
172
factor_matrix_fold(mpz_ptr z,factor_matrix f)173 void factor_matrix_fold(mpz_ptr z, factor_matrix f)
174 {
175 int size;
176 for (size = 0; f [size][0] != 0; size++);
177 if (size == 0) {
178 mpz_set_ui(z, 1);
179 }
180 else
181 factor_matrix_fold_inner(z, f, 0, size);
182 }
183
184
mpz_remove_squares(mpz_ptr D)185 void mpz_remove_squares(mpz_ptr D)
186 {
187 factor_matrix f;
188 int i;
189
190 f = factor_trialdiv (D, 0);
191 for (i = 0; f [i][0] != 0; i++)
192 f [i][1] %= 2;
193
194 factor_matrix_fold (D, f);
195
196 factor_matrix_free (f);
197 }
198
199