Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /*********************************************************************/
16 : /** DIGITS / SUM OF DIGITS **/
17 : /*********************************************************************/
18 : #include "pari.h"
19 : #include "paripriv.h"
20 :
21 : /* set v[i] = 1 iff B^i is needed in the digits_dac algorithm */
22 : static void
23 3194262 : set_vexp(GEN v, long l)
24 : {
25 : long m;
26 3194262 : if (v[l]) return;
27 1149554 : v[l] = 1; m = l>>1;
28 1149554 : set_vexp(v, m);
29 1149554 : set_vexp(v, l-m);
30 : }
31 :
32 : /* return all needed B^i for DAC algorithm, for lz digits */
33 : static GEN
34 895154 : get_vB(GEN T, long lz, void *E, struct bb_ring *r)
35 : {
36 895154 : GEN vB, vexp = const_vecsmall(lz, 0);
37 895154 : long i, l = (lz+1) >> 1;
38 895154 : vexp[1] = 1;
39 895154 : vexp[2] = 1;
40 895154 : set_vexp(vexp, lz);
41 895154 : vB = zerovec(lz); /* unneeded entries remain = 0 */
42 895154 : gel(vB, 1) = T;
43 2204711 : for (i = 2; i <= l; i++)
44 1309557 : if (vexp[i])
45 : {
46 1149554 : long j = i >> 1;
47 1149554 : GEN B2j = r->sqr(E, gel(vB,j));
48 1149554 : gel(vB,i) = odd(i)? r->mul(E, B2j, T): B2j;
49 : }
50 895154 : return vB;
51 : }
52 :
53 : static void
54 369458 : gen_digits_dac(GEN x, GEN vB, long l, GEN *z,
55 : void *E, GEN div(void *E, GEN a, GEN b, GEN *r))
56 : {
57 : GEN q, r;
58 369458 : long m = l>>1;
59 369458 : if (l==1) { *z=x; return; }
60 175527 : q = div(E, x, gel(vB,m), &r);
61 175527 : gen_digits_dac(r, vB, m, z, E, div);
62 175527 : gen_digits_dac(q, vB, l-m, z+m, E, div);
63 : }
64 :
65 : static GEN
66 51212 : gen_fromdigits_dac(GEN x, GEN vB, long i, long l, void *E,
67 : GEN add(void *E, GEN a, GEN b),
68 : GEN mul(void *E, GEN a, GEN b))
69 : {
70 : GEN a, b;
71 51212 : long m = l>>1;
72 51212 : if (l==1) return gel(x,i);
73 23142 : a = gen_fromdigits_dac(x, vB, i, m, E, add, mul);
74 23142 : b = gen_fromdigits_dac(x, vB, i+m, l-m, E, add, mul);
75 23142 : return add(E, a, mul(E, b, gel(vB, m)));
76 : }
77 :
78 : static GEN
79 19006 : gen_digits_i(GEN x, GEN B, long n, void *E, struct bb_ring *r,
80 : GEN (*div)(void *E, GEN x, GEN y, GEN *r))
81 : {
82 : GEN z, vB;
83 19006 : if (n==1) retmkvec(gcopy(x));
84 18404 : vB = get_vB(B, n, E, r);
85 18404 : z = cgetg(n+1, t_VEC);
86 18404 : gen_digits_dac(x, vB, n, (GEN*)(z+1), E, div);
87 18404 : return z;
88 : }
89 :
90 : GEN
91 18949 : gen_digits(GEN x, GEN B, long n, void *E, struct bb_ring *r,
92 : GEN (*div)(void *E, GEN x, GEN y, GEN *r))
93 : {
94 18949 : pari_sp av = avma;
95 18949 : return gerepilecopy(av, gen_digits_i(x, B, n, E, r, div));
96 : }
97 :
98 : GEN
99 4928 : gen_fromdigits(GEN x, GEN B, void *E, struct bb_ring *r)
100 : {
101 4928 : pari_sp av = avma;
102 4928 : long n = lg(x)-1;
103 4928 : GEN vB = get_vB(B, n, E, r);
104 4928 : GEN z = gen_fromdigits_dac(x, vB, 1, n, E, r->add, r->mul);
105 4928 : return gerepilecopy(av, z);
106 : }
107 :
108 : static GEN
109 1911 : _addii(void *data /* ignored */, GEN x, GEN y)
110 1911 : { (void)data; return addii(x,y); }
111 : static GEN
112 1087359 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
113 : static GEN
114 319870 : _mulii(void *data /* ignored */, GEN x, GEN y)
115 319870 : { (void)data; return mulii(x,y); }
116 : static GEN
117 443 : _dvmdii(void *data /* ignored */, GEN x, GEN y, GEN *r)
118 443 : { (void)data; return dvmdii(x,y,r); }
119 :
120 : static struct bb_ring Z_ring = { _addii, _mulii, _sqri };
121 :
122 : /* does not affect stack unless B = NULL */
123 : static GEN
124 371 : check_basis(GEN B)
125 : {
126 371 : if (!B) return utoipos(10);
127 343 : if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
128 343 : if (abscmpiu(B,2)<0) pari_err_DOMAIN("digits","B","<",gen_2,B);
129 343 : return B;
130 : }
131 :
132 : /* x has l digits in base B, write them to z[0..l-1], vB[i] = B^i */
133 : static void
134 3110 : digits_dacsmall(GEN x, GEN vB, long l, ulong* z)
135 : {
136 3110 : pari_sp av = avma;
137 : GEN q,r;
138 : long m;
139 3110 : if (l==1) { *z=itou(x); return; }
140 1538 : m=l>>1;
141 1538 : q = dvmdii(x, gel(vB,m), &r);
142 1538 : digits_dacsmall(q,vB,l-m,z);
143 1538 : digits_dacsmall(r,vB,m,z+l-m);
144 1538 : set_avma(av);
145 : }
146 :
147 : /* x t_INT */
148 : static GEN
149 70 : digits_i(GEN x, GEN B)
150 : {
151 70 : pari_sp av = avma;
152 : long lz;
153 : GEN z;
154 70 : B = check_basis(B);
155 70 : if (signe(B) < 0) pari_err_DOMAIN("digits","B","<",gen_0,B);
156 70 : if (!signe(x)) {set_avma(av); return cgetg(1,t_VEC); }
157 63 : if (abscmpii(x,B)<0) {set_avma(av); retmkvec(absi(x)); }
158 56 : if (Z_ispow2(B))
159 : {
160 14 : long k = expi(B);
161 14 : if (k == 1) return binaire(x);
162 7 : if (k >= BITS_IN_LONG) return binary_2k(x, k);
163 0 : (void)new_chunk(4*(expi(x) + 2)); /* HACK */
164 0 : z = binary_2k_nv(x, k);
165 0 : set_avma(av); return Flv_to_ZV(z);
166 : }
167 42 : x = absi_shallow(x);
168 42 : lz = logint(x,B) + 1;
169 42 : if (lgefint(B) > 3)
170 : {
171 8 : z = gerepileupto(av, gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii));
172 8 : vecreverse_inplace(z); return z;
173 : }
174 : else
175 : {
176 34 : GEN vB = get_vB(B, lz, NULL, &Z_ring);
177 34 : (void)new_chunk(3*lz); /* HACK */
178 34 : z = zero_zv(lz);
179 34 : digits_dacsmall(x,vB,lz,(ulong*)(z+1));
180 34 : set_avma(av); return Flv_to_ZV(z);
181 : }
182 : }
183 : GEN
184 91 : digits(GEN x, GEN B)
185 : {
186 91 : pari_sp av = avma;
187 91 : long v = 0;
188 91 : if (typ(x) == t_INT) return digits_i(x, B);
189 35 : if (typ(x) != t_PADIC || (v = valp(x)) < 0 || (B && !gequal(B, gel(x,2))))
190 7 : pari_err_TYPE("digits",x);
191 28 : if (!signe(gel(x, 4))) return cgetg(1, t_VEC);
192 14 : x = digits_i(gel(x, 4), gel(x, 2));
193 14 : vecreverse_inplace(x);
194 14 : if (!v) return x;
195 0 : return gerepileupto(av, concat(zerovec(v), x));
196 : }
197 :
198 : static GEN
199 3523090 : fromdigitsu_dac(GEN x, GEN vB, long i, long l)
200 : {
201 : GEN a, b;
202 3523090 : long m = l>>1;
203 3523090 : if (l==1) return utoi(uel(x,i));
204 2916183 : if (l==2) return addui(uel(x,i), mului(uel(x,i+1), gel(vB, m)));
205 1325651 : a = fromdigitsu_dac(x, vB, i, m);
206 1325651 : b = fromdigitsu_dac(x, vB, i+m, l-m);
207 1325651 : return addii(a, mulii(b, gel(vB, m)));
208 : }
209 :
210 : static GEN
211 871788 : fromdigitsu_i(GEN x, GEN B)
212 : {
213 871788 : long n = lg(x)-1;
214 : GEN vB;
215 871788 : if (n == 0) return gen_0;
216 871788 : vB = get_vB(B, n, NULL, &Z_ring);
217 871788 : return fromdigitsu_dac(x, vB, 1, n);
218 : }
219 : GEN
220 871774 : fromdigitsu(GEN x, GEN B)
221 871774 : { pari_sp av = avma; return gerepileuptoint(av, fromdigitsu_i(x, B)); }
222 :
223 : static int
224 28 : ZV_in_range(GEN v, GEN B)
225 : {
226 28 : long i, l = lg(v);
227 1708 : for (i = 1; i < l; i++)
228 : {
229 1694 : GEN vi = gel(v, i);
230 1694 : if (signe(vi) < 0 || cmpii(vi, B) >= 0) return 0;
231 : }
232 14 : return 1;
233 : }
234 : static int
235 21 : zv_nonnegative(GEN v)
236 : {
237 21 : long i, l = lg(v);
238 763 : for (i = 1; i < l; i++) if (v[i] < 0) return 0;
239 14 : return 1;
240 : }
241 :
242 : GEN
243 217 : fromdigits(GEN x, GEN B)
244 : {
245 217 : pari_sp av = avma;
246 217 : long tx = typ(x);
247 217 : if (tx == t_VECSMALL)
248 : {
249 28 : if (lg(x)==1) return gen_0;
250 21 : if (zv_nonnegative(x))
251 : {
252 14 : B = check_basis(B); x = vecsmall_reverse(x);
253 14 : return gerepileuptoint(av, fromdigitsu_i(x, B));
254 : }
255 7 : x = zv_to_ZV(x);
256 : }
257 189 : else if (!is_vec_t(tx) || !RgV_is_ZV(x)) pari_err_TYPE("fromdigits",x);
258 196 : if (lg(x) == 1) return gen_0;
259 189 : B = check_basis(B);
260 189 : if (Z_ispow2(B) && ZV_in_range(x, B)) return fromdigits_2k(x, expi(B));
261 175 : x = vecreverse(x);
262 175 : return gerepileuptoint(av, gen_fromdigits(x, B, NULL, &Z_ring));
263 : }
264 :
265 : static const ulong digsum[] ={
266 : 0,1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9,10,2,3,4,5,6,7,8,9,10,11,3,4,5,6,7,8,
267 : 9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,9,10,11,
268 : 12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,
269 : 12,13,14,15,16,17,18,1,2,3,4,5,6,7,8,9,10,2,3,4,5,6,7,8,9,10,11,3,4,5,6,7,8,
270 : 9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,9,10,11,
271 : 12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,
272 : 12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,2,3,4,5,6,7,8,9,10,11,3,
273 : 4,5,6,7,8,9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,
274 : 9,10,11,12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,
275 : 9,10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,13,14,15,
276 : 16,17,18,19,20,3,4,5,6,7,8,9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,
277 : 11,12,13,14,6,7,8,9,10,11,12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,
278 : 12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,
279 : 19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,19,20,21,4,5,6,7,8,9,
280 : 10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,9,10,11,12,13,14,15,7,8,9,10,11,
281 : 12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,
282 : 11,12,13,14,15,16,17,18,19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,
283 : 18,19,20,21,13,14,15,16,17,18,19,20,21,22,5,6,7,8,9,10,11,12,13,14,6,7,8,9,
284 : 10,11,12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,
285 : 10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,13,14,15,16,
286 : 17,18,19,20,12,13,14,15,16,17,18,19,20,21,13,14,15,16,17,18,19,20,21,22,14,
287 : 15,16,17,18,19,20,21,22,23,6,7,8,9,10,11,12,13,14,15,7,8,9,10,11,12,13,14,
288 : 15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,11,12,13,
289 : 14,15,16,17,18,19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,19,20,
290 : 21,13,14,15,16,17,18,19,20,21,22,14,15,16,17,18,19,20,21,22,23,15,16,17,18,
291 : 19,20,21,22,23,24,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,
292 : 10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,13,14,15,16,
293 : 17,18,19,20,12,13,14,15,16,17,18,19,20,21,13,14,15,16,17,18,19,20,21,22,14,
294 : 15,16,17,18,19,20,21,22,23,15,16,17,18,19,20,21,22,23,24,16,17,18,19,20,21,
295 : 22,23,24,25,8,9,10,11,12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,11,
296 : 12,13,14,15,16,17,18,19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,
297 : 19,20,21,13,14,15,16,17,18,19,20,21,22,14,15,16,17,18,19,20,21,22,23,15,16,
298 : 17,18,19,20,21,22,23,24,16,17,18,19,20,21,22,23,24,25,17,18,19,20,21,22,23,
299 : 24,25,26,9,10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,
300 : 13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,19,20,21,13,14,15,16,17,18,19,
301 : 20,21,22,14,15,16,17,18,19,20,21,22,23,15,16,17,18,19,20,21,22,23,24,16,17,
302 : 18,19,20,21,22,23,24,25,17,18,19,20,21,22,23,24,25,26,18,19,20,21,22,23,24,
303 : 25,26,27
304 : };
305 :
306 : ulong
307 355152 : sumdigitsu(ulong n)
308 : {
309 355152 : ulong s = 0;
310 1361843 : while (n) { s += digsum[n % 1000]; n /= 1000; }
311 355152 : return s;
312 : }
313 :
314 : /* res=array of 9-digits integers, return sum_{0 <= i < l} sumdigits(res[i]) */
315 : static ulong
316 14 : sumdigits_block(ulong *res, long l)
317 : {
318 14 : ulong s = sumdigitsu(*--res);
319 355138 : while (--l > 0) s += sumdigitsu(*--res);
320 14 : return s;
321 : }
322 :
323 : GEN
324 35 : sumdigits(GEN n)
325 : {
326 35 : const long L = (long)(ULONG_MAX / 81);
327 35 : pari_sp av = avma;
328 : ulong *res;
329 : long l;
330 :
331 35 : if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
332 35 : switch(lgefint(n))
333 : {
334 7 : case 2: return gen_0;
335 14 : case 3: return utoipos(sumdigitsu(n[2]));
336 : }
337 14 : res = convi(n, &l);
338 14 : if (l < L)
339 : {
340 14 : ulong s = sumdigits_block(res, l);
341 14 : return gc_utoipos(av, s);
342 : }
343 : else /* Huge. Overflows ulong */
344 : {
345 0 : GEN S = gen_0;
346 0 : while (l > L)
347 : {
348 0 : S = addiu(S, sumdigits_block(res, L));
349 0 : res += L; l -= L;
350 : }
351 0 : if (l)
352 0 : S = addiu(S, sumdigits_block(res, l));
353 0 : return gerepileuptoint(av, S);
354 : }
355 : }
356 :
357 : GEN
358 126 : sumdigits0(GEN x, GEN B)
359 : {
360 126 : pari_sp av = avma;
361 : GEN z;
362 : long lz;
363 :
364 126 : if (!B) return sumdigits(x);
365 98 : if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
366 98 : B = check_basis(B);
367 98 : if (Z_ispow2(B))
368 : {
369 35 : long k = expi(B);
370 35 : if (k == 1) return gc_utoi(av,hammingweight(x));
371 21 : if (k < BITS_IN_LONG)
372 : {
373 14 : GEN z = binary_2k_nv(x, k);
374 14 : if (lg(z)-1 > 1L<<(BITS_IN_LONG-k)) /* may overflow */
375 0 : return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
376 14 : return gc_utoi(av,zv_sum(z));
377 : }
378 7 : return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
379 : }
380 63 : if (!signe(x)) { set_avma(av); return gen_0; }
381 63 : if (abscmpii(x,B)<0) { set_avma(av); return absi(x); }
382 56 : if (absequaliu(B,10)) { set_avma(av); return sumdigits(x); }
383 49 : x = absi_shallow(x); lz = logint(x,B) + 1;
384 49 : z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
385 49 : return gerepileuptoint(av, ZV_sum(z));
386 : }
|
