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 3831527 : set_vexp(GEN v, long l)
24 : {
25 : long m;
26 3831527 : if (v[l]) return;
27 1399702 : v[l] = 1; m = l>>1;
28 1399702 : set_vexp(v, m);
29 1399702 : set_vexp(v, l-m);
30 : }
31 :
32 : /* return all needed B^i for DAC algorithm, for lz digits */
33 : static GEN
34 1032123 : get_vB(GEN T, long lz, void *E, struct bb_ring *r)
35 : {
36 1032123 : GEN vB, vexp = const_vecsmall(lz, 0);
37 1032123 : long i, l = (lz+1) >> 1;
38 1032123 : vexp[1] = 1;
39 1032123 : vexp[2] = 1;
40 1032123 : set_vexp(vexp, lz);
41 1032123 : vB = zerovec(lz); /* unneeded entries remain = 0 */
42 1032123 : gel(vB, 1) = T;
43 2608071 : for (i = 2; i <= l; i++)
44 1575948 : if (vexp[i])
45 : {
46 1399702 : long j = i >> 1;
47 1399702 : GEN B2j = r->sqr(E, gel(vB,j));
48 1399702 : gel(vB,i) = odd(i)? r->mul(E, B2j, T): B2j;
49 : }
50 1032123 : return vB;
51 : }
52 :
53 : static void
54 327640 : 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 327640 : long m = l>>1;
59 327640 : if (l==1) { *z=x; return; }
60 156508 : q = div(E, x, gel(vB,m), &r);
61 156508 : gen_digits_dac(r, vB, m, z, E, div);
62 156508 : gen_digits_dac(q, vB, l-m, z+m, E, div);
63 : }
64 :
65 : static GEN
66 972580 : 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 972580 : long m = l>>1;
72 972580 : if (l==1) return gel(x,i);
73 441826 : a = gen_fromdigits_dac(x, vB, i, m, E, add, mul);
74 441826 : b = gen_fromdigits_dac(x, vB, i+m, l-m, E, add, mul);
75 441826 : return add(E, a, mul(E, b, gel(vB, m)));
76 : }
77 :
78 : static GEN
79 15240 : 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 15240 : if (n==1) retmkvec(gcopy(x));
84 14624 : vB = get_vB(B, n, E, r);
85 14624 : z = cgetg(n+1, t_VEC);
86 14624 : gen_digits_dac(x, vB, n, (GEN*)(z+1), E, div);
87 14624 : return z;
88 : }
89 :
90 : GEN
91 15183 : 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 15183 : pari_sp av = avma;
95 15183 : return gerepilecopy(av, gen_digits_i(x, B, n, E, r, div));
96 : }
97 :
98 : GEN
99 88928 : gen_fromdigits(GEN x, GEN B, void *E, struct bb_ring *r)
100 : {
101 88928 : pari_sp av = avma;
102 88928 : long n = lg(x)-1;
103 88928 : GEN vB = get_vB(B, n, E, r);
104 88928 : GEN z = gen_fromdigits_dac(x, vB, 1, n, E, r->add, r->mul);
105 88928 : return gerepilecopy(av, z);
106 : }
107 :
108 : static GEN
109 420595 : _addii(void *data /* ignored */, GEN x, GEN y)
110 420595 : { (void)data; return addii(x,y); }
111 : static GEN
112 1345389 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
113 : static GEN
114 827668 : _mulii(void *data /* ignored */, GEN x, GEN y)
115 827668 : { (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 224469 : check_basis(GEN B)
125 : {
126 224469 : if (!B) return utoipos(10);
127 224441 : if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
128 224441 : if (abscmpiu(B,2)<0) pari_err_DOMAIN("digits","abs(B)","<",gen_2,B);
129 224427 : 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 479894 : digits_dacsmall(GEN x, GEN vB, long l, ulong* z)
135 : {
136 479894 : pari_sp av = avma;
137 : GEN q,r;
138 : long m;
139 479894 : if (l==1) { *z=itou(x); return; }
140 212042 : m=l>>1;
141 212042 : q = dvmdii(x, gel(vB,m), &r);
142 212042 : digits_dacsmall(q,vB,l-m,z);
143 212042 : digits_dacsmall(r,vB,m,z+l-m);
144 212042 : set_avma(av);
145 : }
146 :
147 : /* x t_INT */
148 : static GEN
149 112091 : digits_pos(GEN x, GEN B)
150 : {
151 : long lz;
152 : GEN z;
153 112091 : if (abscmpii(x,B)<0) retmkvec(absi(x));
154 111748 : if (Z_ispow2(B))
155 : {
156 55930 : pari_sp av = avma;
157 55930 : long k = expi(B);
158 55930 : if (k == 1) return binaire(x);
159 27923 : if (k >= BITS_IN_LONG) return binary_2k(x, k);
160 27916 : (void)new_chunk(4*(expi(x) + 2)); /* HACK */
161 27916 : z = binary_2k_nv(x, k);
162 27916 : set_avma(av); return Flv_to_ZV(z);
163 : }
164 55818 : x = absi_shallow(x);
165 55818 : lz = logint(x,B) + 1;
166 55818 : if (lgefint(B) > 3)
167 : {
168 8 : z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
169 8 : vecreverse_inplace(z); return z;
170 : }
171 : else
172 : {
173 55810 : pari_sp av = avma;
174 55810 : GEN vB = get_vB(B, lz, NULL, &Z_ring);
175 55810 : (void)new_chunk(3*lz); /* HACK */
176 55810 : z = zero_zv(lz);
177 55810 : digits_dacsmall(x,vB,lz,(ulong*)(z+1));
178 55810 : set_avma(av); return Flv_to_ZV(z);
179 : }
180 : }
181 :
182 : static GEN
183 56028 : digits_neg(GEN n, GEN B)
184 : {
185 56028 : GEN V = digits_pos(n,B);
186 56028 : long i, l = lg(V), carry = 0, s = signe(n)==-1;
187 380478 : for (i = 1; i < l; i++)
188 : {
189 324450 : GEN u = odd(i+s) ? addsi(carry, gel(V, l-i)): subsi(carry,gel(V, l-i));
190 324450 : if (signe(u) < 0) { u = addii(u, B); carry=1; }
191 : else
192 203721 : if (cmpii(u, B) >= 0) { u = subii(u, B); carry=-1; }
193 : else
194 170457 : carry = 0;
195 324450 : gel(V,l-i) = u;
196 : }
197 56028 : if (carry > 0)
198 28014 : V = vec_prepend(V, stoi(carry));
199 28014 : else if (carry < 0)
200 8533 : V = shallowconcat(mkvec2(gen_1,addsi(carry, B)),V);
201 56028 : return V;
202 : }
203 :
204 : /* x t_INT */
205 : static GEN
206 112154 : digits_i(GEN x, GEN B)
207 : {
208 112154 : pari_sp av = avma;
209 112154 : B = check_basis(B);
210 112140 : if (!signe(x)) {set_avma(av); return cgetg(1,t_VEC); }
211 112077 : if (signe(B) > 0) return gerepilecopy(av, digits_pos(x, B));
212 56014 : return gerepilecopy(av, digits_neg(x, negi(B)));
213 : }
214 :
215 : GEN
216 112175 : digits(GEN x, GEN B)
217 : {
218 112175 : pari_sp av = avma;
219 112175 : long v = 0;
220 112175 : if (typ(x) == t_INT) return digits_i(x, B);
221 35 : if (typ(x) != t_PADIC || (v = valp(x)) < 0 || (B && !gequal(B, gel(x,2))))
222 7 : pari_err_TYPE("digits",x);
223 28 : if (!signe(gel(x, 4))) return cgetg(1, t_VEC);
224 14 : x = digits_i(gel(x, 4), gel(x, 2));
225 14 : vecreverse_inplace(x);
226 14 : if (!v) return x;
227 0 : return gerepileupto(av, concat(zerovec(v), x));
228 : }
229 :
230 : static GEN
231 3525149 : fromdigitsu_dac(GEN x, GEN vB, long i, long l)
232 : {
233 : GEN a, b;
234 3525149 : long m = l>>1;
235 3525149 : if (l==1) return utoi(uel(x,i));
236 2917797 : if (l==2) return addui(uel(x,i), mului(uel(x,i+1), gel(vB, m)));
237 1326194 : a = fromdigitsu_dac(x, vB, i, m);
238 1326194 : b = fromdigitsu_dac(x, vB, i+m, l-m);
239 1326194 : return addii(a, mulii(b, gel(vB, m)));
240 : }
241 :
242 : static GEN
243 872761 : fromdigitsu_i(GEN x, GEN B)
244 : {
245 872761 : long n = lg(x)-1;
246 : GEN vB;
247 872761 : if (n == 0) return gen_0;
248 872761 : vB = get_vB(B, n, NULL, &Z_ring);
249 872761 : return fromdigitsu_dac(x, vB, 1, n);
250 : }
251 : GEN
252 872747 : fromdigitsu(GEN x, GEN B)
253 872747 : { pari_sp av = avma; return gerepileuptoint(av, fromdigitsu_i(x, B)); }
254 :
255 : static int
256 28028 : ZV_in_range(GEN v, GEN B)
257 : {
258 28028 : long i, l = lg(v);
259 220822 : for (i = 1; i < l; i++)
260 : {
261 192808 : GEN vi = gel(v, i);
262 192808 : if (signe(vi) < 0 || cmpii(vi, B) >= 0) return 0;
263 : }
264 28014 : return 1;
265 : }
266 : static int
267 21 : zv_nonnegative(GEN v)
268 : {
269 21 : long i, l = lg(v);
270 763 : for (i = 1; i < l; i++) if (v[i] < 0) return 0;
271 14 : return 1;
272 : }
273 :
274 : GEN
275 112273 : fromdigits(GEN x, GEN B)
276 : {
277 112273 : pari_sp av = avma;
278 112273 : long tx = typ(x);
279 112273 : if (tx == t_VECSMALL)
280 : {
281 28 : if (lg(x)==1) return gen_0;
282 21 : if (zv_nonnegative(x))
283 : {
284 14 : B = check_basis(B); x = vecsmall_reverse(x);
285 14 : return gerepileuptoint(av, fromdigitsu_i(x, B));
286 : }
287 7 : x = zv_to_ZV(x);
288 : }
289 112245 : else if (!is_vec_t(tx) || !RgV_is_ZV(x)) pari_err_TYPE("fromdigits",x);
290 112252 : if (lg(x) == 1) return gen_0;
291 112189 : B = check_basis(B);
292 112189 : if (Z_ispow2(B) && ZV_in_range(x, B)) return fromdigits_2k(x, expi(B));
293 84175 : x = vecreverse(x);
294 84175 : return gerepileuptoint(av, gen_fromdigits(x, B, NULL, &Z_ring));
295 : }
296 :
297 : static const ulong digsum[] ={
298 : 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,
299 : 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,
300 : 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,
301 : 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,
302 : 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,
303 : 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,
304 : 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,
305 : 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,
306 : 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,
307 : 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,
308 : 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,
309 : 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,
310 : 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,
311 : 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,
312 : 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,
313 : 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,
314 : 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,
315 : 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,
316 : 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,
317 : 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,
318 : 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,
319 : 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,
320 : 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,
321 : 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,
322 : 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,
323 : 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,
324 : 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,
325 : 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,
326 : 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,
327 : 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,
328 : 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,
329 : 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,
330 : 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,
331 : 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,
332 : 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,
333 : 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,
334 : 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,
335 : 25,26,27
336 : };
337 :
338 : ulong
339 355152 : sumdigitsu(ulong n)
340 : {
341 355152 : ulong s = 0;
342 1361843 : while (n) { s += digsum[n % 1000]; n /= 1000; }
343 355152 : return s;
344 : }
345 :
346 : /* res=array of 9-digits integers, return sum_{0 <= i < l} sumdigits(res[i]) */
347 : static ulong
348 14 : sumdigits_block(ulong *res, long l)
349 : {
350 14 : ulong s = sumdigitsu(*--res);
351 355138 : while (--l > 0) s += sumdigitsu(*--res);
352 14 : return s;
353 : }
354 :
355 : GEN
356 35 : sumdigits(GEN n)
357 : {
358 35 : const long L = (long)(ULONG_MAX / 81);
359 35 : pari_sp av = avma;
360 : ulong *res;
361 : long l;
362 :
363 35 : if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
364 35 : switch(lgefint(n))
365 : {
366 7 : case 2: return gen_0;
367 14 : case 3: return utoipos(sumdigitsu(n[2]));
368 : }
369 14 : res = convi(n, &l);
370 14 : if (l < L)
371 : {
372 14 : ulong s = sumdigits_block(res, l);
373 14 : return gc_utoipos(av, s);
374 : }
375 : else /* Huge. Overflows ulong */
376 : {
377 0 : GEN S = gen_0;
378 0 : while (l > L)
379 : {
380 0 : S = addiu(S, sumdigits_block(res, L));
381 0 : res += L; l -= L;
382 : }
383 0 : if (l)
384 0 : S = addiu(S, sumdigits_block(res, l));
385 0 : return gerepileuptoint(av, S);
386 : }
387 : }
388 :
389 : GEN
390 140 : sumdigits0(GEN x, GEN B)
391 : {
392 140 : pari_sp av = avma;
393 : GEN z;
394 : long lz;
395 :
396 140 : if (!B) return sumdigits(x);
397 112 : if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
398 112 : B = check_basis(B);
399 112 : if (signe(B)<0) return gerepileuptoint(av, ZV_sum(digits_neg(x,negi(B))));
400 98 : if (Z_ispow2(B))
401 : {
402 35 : long k = expi(B);
403 35 : if (k == 1) return gc_utoi(av,hammingweight(x));
404 21 : if (k < BITS_IN_LONG)
405 : {
406 14 : GEN z = binary_2k_nv(x, k);
407 14 : if (lg(z)-1 > 1L<<(BITS_IN_LONG-k)) /* may overflow */
408 0 : return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
409 14 : return gc_utoi(av,zv_sum(z));
410 : }
411 7 : return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
412 : }
413 63 : if (!signe(x)) { set_avma(av); return gen_0; }
414 63 : if (abscmpii(x,B)<0) { set_avma(av); return absi(x); }
415 56 : if (absequaliu(B,10)) { set_avma(av); return sumdigits(x); }
416 49 : x = absi_shallow(x); lz = logint(x,B) + 1;
417 49 : z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
418 49 : return gerepileuptoint(av, ZV_sum(z));
419 : }
|
