Line data Source code
1 : /* Copyright (C) 2002 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 : /* Original code contributed by: Ralf Stephan (ralf@ark.in-berlin.de).
16 : * Updated by Bill Allombert (2014) to use Selberg formula for L
17 : * following http://dx.doi.org/10.1112/S1461157012001088
18 : *
19 : * This program is basically the implementation of the script
20 : *
21 : * Psi(n, q) = my(a=sqrt(2/3)*Pi/q, b=n-1/24, c=sqrt(b));
22 : * (sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c))
23 : * L(n,q)=sqrt(k/3)*sum(l=0,2*k-1,
24 : if(((3*l^2+l)/2+n)%k==0,(-1)^l*cos((6*l+1)/(6*k)*Pi)))
25 : * part(n) = round(sum(q=1,5 + 0.24*sqrt(n),L(n,q)*Psi(n,q)))
26 : *
27 : * only faster.
28 : *
29 : * ------------------------------------------------------------------
30 : * The first restriction depends on Pari's maximum precision of floating
31 : * point reals, which is 268435454 bits in 2.2.4, since the algorithm needs
32 : * high precision exponentials. For that engine, the maximum possible argument
33 : * would be in [5*10^15,10^16], the computation of which would need days on
34 : * a ~1-GHz computer. */
35 :
36 : #include "pari.h"
37 : #include "paripriv.h"
38 :
39 : /****************************************************************/
40 :
41 : /* Given c = sqrt(2/3)*Pi*sqrt(N-1/24)
42 : * Psi(N, q) = my(a = c/q); sqrt(q) * (a*cosh(a) - sinh(a)) */
43 : static GEN
44 1109212 : psi(GEN c, ulong q, long prec)
45 : {
46 1109212 : GEN a = divru(c, q), ea = mpexp(a), invea = invr(ea);
47 1109212 : GEN cha = shiftr(addrr(ea, invea), -1); /* ch(a) */
48 1109212 : GEN sha = shiftr(subrr(ea, invea), -1); /* sh(a) */
49 1109212 : return mulrr(sqrtr(utor(q,prec)), subrr(mulrr(a,cha), sha));
50 : }
51 :
52 : /* L(n,q)=sqrt(k/3)*sum(l=0,2*k-1,
53 : if(((3*l^2+l)/2+n)%k==0,(-1)^l*cos((6*l+1)/(6*k)*Pi)))
54 : * Never called with q < 3, so ignore this case */
55 : static GEN
56 1295452 : L(GEN n, ulong k, long bitprec)
57 : {
58 : ulong r, l, m;
59 1295452 : long pr = nbits2prec(bitprec / k + k);
60 1295452 : GEN s = utor(0,pr), pi = mppi(pr);
61 1295452 : pari_sp av = avma;
62 :
63 1295452 : r = 2; m = umodiu(n,k);
64 34034238 : for (l = 0; l < 2*k; l++)
65 : {
66 32738786 : if (m == 0)
67 : {
68 2590332 : GEN c = mpcos(divru(mulru(pi, 6*l+1), 6*k));
69 2590332 : if (odd(l)) subrrz(s, c, s); else addrrz(s, c, s);
70 2590332 : set_avma(av);
71 : }
72 32738786 : m += r; if (m >= k) m -= k;
73 32738786 : r += 3; if (r >= k) r -= k;
74 : }
75 : /* multiply by sqrt(k/3) */
76 1295452 : return mulrr(s, sqrtr((k % 3)? rdivss(k,3,pr): utor(k/3,pr)));
77 : }
78 :
79 : /* Return a low precision estimate of log p(n). */
80 : static GEN
81 70000 : estim(GEN n)
82 : {
83 70000 : pari_sp av = avma;
84 70000 : GEN p1, pi = mppi (DEFAULTPREC);
85 :
86 70000 : p1 = divru( itor(shifti(n,1), DEFAULTPREC), 3 );
87 70000 : p1 = mpexp( mulrr(pi, sqrtr(p1)) ); /* exp(Pi * sqrt(2N/3)) */
88 70000 : p1 = divri (shiftr(p1,-2), n);
89 70000 : p1 = divrr(p1, sqrtr( utor(3,DEFAULTPREC) ));
90 70000 : return gerepileupto(av, mplog(p1));
91 : }
92 :
93 : /* c = sqrt(2/3)*Pi*sqrt(n-1/24); d = 1 / ((2*b)^(3/2) * Pi); */
94 : static void
95 70000 : pinit(GEN n, GEN *c, GEN *d, ulong prec)
96 : {
97 70000 : GEN b = divru( itor( subiu(muliu(n,24), 1), prec ), 24 ); /* n - 1/24 */
98 70000 : GEN sqrtb = sqrtr(b), Pi = mppi(prec), pi2sqrt2, pisqrt2d3;
99 :
100 70000 : pisqrt2d3 = mulrr(Pi, sqrtr( divru(utor(2, prec), 3) ));
101 70000 : pi2sqrt2 = mulrr(Pi, sqrtr( utor(8, prec) ));
102 70000 : *c = mulrr(pisqrt2d3, sqrtb);
103 70000 : *d = invr( mulrr(pi2sqrt2, mulrr(b,sqrtb)) );
104 70000 : }
105 :
106 : /* part(n) = round(sum(q=1,5 + 0.24*sqrt(n), L(n,q)*Psi(n,q))) */
107 : GEN
108 70021 : numbpart(GEN n)
109 : {
110 70021 : pari_sp ltop = avma, av;
111 : GEN sum, est, C, D, p1, p2;
112 : long prec, bitprec;
113 : ulong q;
114 :
115 70021 : if (typ(n) != t_INT) pari_err_TYPE("partition function",n);
116 70021 : if (signe(n) < 0) return gen_0;
117 70021 : if (abscmpiu(n, 2) < 0) return gen_1;
118 70007 : if (cmpii(n, uu32toi(0x38d7e, 0xa4c68000)) >= 0)
119 7 : pari_err_OVERFLOW("numbpart [n < 10^15]");
120 70000 : est = estim(n);
121 70000 : bitprec = (long)(rtodbl(est)/M_LN2) + 32;
122 70000 : prec = nbits2prec(bitprec);
123 70000 : pinit(n, &C, &D, prec);
124 70000 : sum = utor(0, prec);
125 : /* Because N < 10^16 and q < sqrt(N), q fits into a long
126 : * In fact q < 2 LONG_MAX / 3 */
127 70000 : av = avma; togglesign(est);
128 1365452 : for (q = (ulong)(sqrt(gtodouble(n))*0.24 + 5); q >= 3; q--, set_avma(av))
129 : {
130 1295452 : GEN t = L(n, q, bitprec);
131 1295452 : if (abscmprr(t, mpexp(divru(est,q))) < 0) continue;
132 :
133 969212 : t = mulrr(t, psi(gprec_w(C, nbits2prec(bitprec / q + 32)), q, prec));
134 969212 : affrr(addrr(sum, t), sum);
135 : }
136 70000 : p1 = addrr(sum, psi(C, 1, prec));
137 70000 : p2 = psi(C, 2, prec);
138 70000 : affrr(mod2(n)? subrr(p1,p2): addrr(p1,p2), sum);
139 70000 : return gerepileuptoint (ltop, roundr(mulrr(D,sum)));
140 : }
141 :
142 : /* for loop over partitions of integer k.
143 : * nbounds can restrict partitions to have length between nmin and nmax
144 : * (the length is the number of non zero entries) and
145 : * abounds restrict to integers between amin and amax.
146 : *
147 : * Start from central partition.
148 : * By default, remove zero entries on the left.
149 : *
150 : * Algorithm:
151 : *
152 : * A partition of k is an increasing sequence v1,... vn with sum(vi)=k
153 : * The starting point is the minimal n-partition of k: a,...a,a+1,.. a+1
154 : * (a+1 is repeated r times with k = a * n + r).
155 : *
156 : * The procedure to obtain the next partition:
157 : * - find the last index i<n such that v{i-1} != v{i} (that is vi is the start
158 : * of the last constant range excluding vn).
159 : * - decrease vi by one, and set v{i+1},... v{n} to be a minimal partition (of
160 : * the right sum).
161 : *
162 : * Examples: we highlight the index i
163 : * 1 1 2 2 3
164 : * ^
165 : * 1 1 1 3 3
166 : * ^
167 : * 1 1 1 2 4
168 : * ^
169 : * 1 1 1 1 5
170 : * ^
171 : * 0 2 2 2 3
172 : * ^
173 : * This is recursive in nature. Restrictions on upper bounds of the vi or on
174 : * the length of the partitions are straightforward to take into account. */
175 :
176 : static void
177 217 : parse_interval(GEN a, long *amin, long *amax)
178 : {
179 217 : switch (typ(a))
180 : {
181 49 : case t_INT:
182 49 : *amax = itos(a);
183 49 : break;
184 168 : case t_VEC:
185 168 : if (lg(a) != 3)
186 0 : pari_err_TYPE("forpart [expect vector of type [amin,amax]]",a);
187 168 : *amin = gtos(gel(a,1));
188 168 : *amax = gtos(gel(a,2));
189 168 : if (*amin>*amax || *amin<0 || *amax<=0)
190 0 : pari_err_TYPE("forpart [expect 0<=min<=max, 0<max]",a);
191 168 : break;
192 0 : default:
193 0 : pari_err_TYPE("forpart",a);
194 : }
195 217 : }
196 :
197 : void
198 238 : forpart_init(forpart_t *T, long k, GEN abound, GEN nbound)
199 : {
200 :
201 : /* bound on coefficients */
202 238 : T->amin=1;
203 238 : if (abound) parse_interval(abound,&T->amin,&T->amax);
204 126 : else T->amax = k;
205 : /* strip leading zeros ? */
206 238 : T->strip = (T->amin > 0) ? 1 : 0;
207 : /* bound on number of nonzero coefficients */
208 238 : T->nmin=0;
209 238 : if (nbound) parse_interval(nbound,&T->nmin,&T->nmax);
210 133 : else T->nmax = k;
211 :
212 : /* non empty if nmin*amin <= k <= amax*nmax */
213 238 : if ( T->amin*T->nmin > k || k > T->amax * T->nmax )
214 : {
215 0 : T->nmin = T->nmax = 0;
216 : }
217 : else
218 : {
219 : /* to reach nmin one must have k <= nmin*amax, otherwise increase nmin */
220 238 : if ( T->nmin * T->amax < k )
221 154 : T->nmin = 1 + (k - 1) / T->amax; /* ceil( k/tmax ) */
222 : /* decrease nmax (if strip): k <= amin*nmax */
223 238 : if (T->strip && T->nmax > k/T->amin)
224 14 : T->nmax = k / T->amin; /* strip implies amin>0 */ /* fixme: take ceil( ) */
225 : /* no need to change amin */
226 : /* decrease amax if amax + (nmin-1)*amin > k */
227 238 : if ( T->amax + (T->nmin-1)* T->amin > k )
228 56 : T->amax = k - (T->nmin-1)* T->amin;
229 : }
230 :
231 238 : if ( T->amax < T->amin )
232 21 : T->nmin = T->nmax = 0;
233 :
234 238 : T->v = zero_zv(T->nmax); /* partitions will be of length <= nmax */
235 238 : T->k = k;
236 238 : }
237 :
238 : GEN
239 3161753 : forpart_next(forpart_t *T)
240 : {
241 3161753 : GEN v = T->v;
242 3161753 : long n = lg(v)-1;
243 : long i, s, a, k, vi, vn;
244 :
245 3161753 : if (n>0 && v[n])
246 : {
247 : /* find index to increase: i s.t. v[i+1],...v[n] is central a,..a,a+1,..a+1
248 : keep s = v[i] + v[i+1] + ... + v[n] */
249 3161508 : s = a = v[n];
250 4638872 : for(i = n-1; i>0 && v[i]+1 >= a; s += v[i--]);
251 3161508 : if (i == 0) {
252 : /* v is central [ a, a, .. a, a+1, .. a+1 ] */
253 1120 : if ((n+1) * T->amin > s || n == T->nmax) return NULL;
254 980 : i = 1; n++;
255 980 : setlg(v, n+1);
256 980 : vi = T->amin;
257 : } else {
258 3160388 : s += v[i];
259 3160388 : vi = v[i]+1;
260 : }
261 : } else {
262 : /* init v */
263 245 : s = T->k;
264 245 : if (T->amin == 0) T->amin = 1;
265 245 : if (T->strip) { n = T->nmin; setlg(T->v, n+1); }
266 245 : if (s==0)
267 : {
268 21 : if (n==0 && T->nmin==0) {T->nmin++; return v;}
269 7 : return NULL;
270 : }
271 224 : if (n==0) return NULL;
272 217 : vi = T->amin;
273 217 : i = T->strip ? 1 : n + 1 - T->nmin; /* first nonzero index */
274 217 : if (s <= (n-i)*vi) return NULL;
275 : }
276 : /* now fill [ v[i],... v[n] ] with s, start at vi */
277 3161571 : vn = s - (n-i)*vi; /* expected value for v[n] */
278 3161571 : if (T->amax && vn > T->amax)
279 140 : {
280 : /* do not exceed amax */
281 : long ai, q, r;
282 140 : vn -= vi;
283 140 : ai = T->amax - vi;
284 140 : q = vn / ai; /* number of nmax */
285 140 : r = vn % ai; /* value before nmax */
286 : /* fill [ v[i],... v[n] ] as [ vi,... vi, vi+r, amax,... amax ] */
287 469 : while ( q-- ) v[n--] = T->amax;
288 140 : if ( n >= i ) v[n--] = vi + r;
289 371 : while ( n >= i ) v[n--] = vi;
290 : } else {
291 : /* fill as [ v[i], ... v[i], vn ] */
292 7798840 : for ( k=i; k<n; v[k++] = vi );
293 3161431 : v[n] = vn;
294 : }
295 3161571 : return v;
296 : }
297 :
298 : GEN
299 0 : forpart_prev(forpart_t *T)
300 : {
301 0 : GEN v = T->v;
302 0 : long n = lg(v)-1;
303 : long j, ni, q, r;
304 : long i, s;
305 0 : if (n>0 && v[n])
306 : {
307 : /* find index to decrease: start of last constant sequence, excluding v[n] */
308 0 : i = n-1; s = v[n];
309 0 : while (i>1 && (v[i-1]==v[i] || v[i+1]==T->amax))
310 0 : s+= v[i--];
311 0 : if (!i) return NULL;
312 : /* amax condition: cannot decrease i if maximal on the right */
313 0 : if ( v[i+1] == T->amax ) return NULL;
314 : /* amin condition: stop if below except if strip & try to remove */
315 0 : if (v[i] == T->amin) {
316 0 : if (!T->strip) return NULL;
317 0 : s += v[i]; v[i] = 0;
318 : } else {
319 0 : v[i]--; s++;
320 : }
321 : /* zero case... */
322 0 : if (v[i] == 0)
323 : {
324 0 : if (T->nmin > n-i) return NULL; /* need too many non zero coeffs */
325 : /* reduce size of v ? */
326 0 : if (T->strip) {
327 0 : i = 0; n--;
328 0 : setlg(v, n+1);
329 : }
330 : }
331 : } else
332 : {
333 0 : s = T->k;
334 0 : i = 0;
335 0 : if (s==0)
336 : {
337 0 : if (n==0 && T->nmin==0) {T->nmin++; return v;}
338 0 : return NULL;
339 : }
340 0 : if (n*T->amax < s || s < T->nmin*T->amin) return NULL;
341 : }
342 : /* set minimal partition of sum s starting from index i+1 */
343 0 : ni = n-i;
344 0 : q = s / ni;
345 0 : r = s % ni;
346 0 : for(j=i+1; j<=n-r; j++) v[j]=q;
347 0 : for(j=n-r+1; j<=n; j++) v[j]=q + 1;
348 0 : return v;
349 : }
350 :
351 : static long
352 98 : countpart(long k, GEN abound, GEN nbound)
353 : {
354 98 : pari_sp av = avma;
355 : long n;
356 : forpart_t T;
357 98 : if (k<0) return 0;
358 91 : forpart_init(&T, k, abound, nbound);
359 819 : for (n=0; forpart_next(&T); n++)
360 728 : set_avma(av);
361 91 : return n;
362 : }
363 :
364 : GEN
365 98 : partitions(long k, GEN abound, GEN nbound)
366 : {
367 : GEN v;
368 : forpart_t T;
369 98 : long i, n = countpart(k,abound,nbound);
370 98 : if (n==0) return cgetg(1, t_VEC);
371 70 : forpart_init(&T, k, abound, nbound);
372 70 : v = cgetg(n+1, t_VEC);
373 798 : for (i=1; i<=n; i++)
374 728 : gel(v,i)=zv_copy(forpart_next(&T));
375 70 : return v;
376 : }
377 :
378 : void
379 77 : forpart(void *E, long call(void*, GEN), long k, GEN abound, GEN nbound)
380 : {
381 77 : pari_sp av = avma;
382 : GEN v;
383 : forpart_t T;
384 77 : forpart_init(&T, k, abound, nbound);
385 3160206 : while ((v=forpart_next(&T)))
386 3160129 : if (call(E, v)) break;
387 77 : set_avma(av);
388 77 : }
389 :
390 : void
391 84 : forpart0(GEN k, GEN code, GEN abound, GEN nbound)
392 : {
393 84 : pari_sp av = avma;
394 84 : if (typ(k) != t_INT) pari_err_TYPE("forpart",k);
395 84 : if (signe(k)<0) return;
396 77 : push_lex(gen_0, code);
397 77 : forpart((void*)code, &gp_evalvoid, itos(k), abound, nbound);
398 77 : pop_lex(1);
399 77 : set_avma(av);
400 : }
|
