1:- module(int_part, []). 2 3:- use_module(library(clpfd)). 4:- use_module(util(misc)). 5:- use_module(pac(basic)). 6% :- expects_dialect(pac). 7term_expansion --> pac:expand_pac. 8:- use_module(pac(op)). 9 10 11 12% ?- prolog_stack_property(X, Y). 13% ?- set_prolog_stack(global, limit(2*10**9)). 14 15% ?- fold_poly_zip_bounded([], [], 1, [], X). 16% ?- fold_poly_zip_bounded([poly(1,[1,1])], [poly(1,[1,1])], 3, [], X). 17% ?- fold_poly_zip_bounded([poly(1,[1,1])], [poly(1,[1,1])], 5, [], X). 18% ?- fold_poly_zip_bounded([poly(1,[1,1])], [poly(1,[1,1])], 6, [], X). 19 20fold_poly_zip_bounded(Ps, Qs, N, X, Y):- 21 fold_poly_zip(pred(N, ([K, K, P, Q, R] 22 :- math:poly_mul_bd_poly(P, Q, N, R))), 23 Ps, Qs, 0, X, Y). 24 25% ?- scm(int_part). 26partition_number(N, P, Opt):- 27 ( memberchk(formula, Opt) -> partition_number(N, P) 28 ; memberchk(rec, Opt) -> integer_partition(N, Ps), 29 length(Ps, P) 30 ; memberchk(history, Opt) -> course_of_partition(N, [A|_]), 31 length(A, P) 32 ; memberchk(history2, Opt) -> partition_history_tr(N, [A|_]), 33 length(A, P) 34 ; memberchk(clpfd, Opt) -> setof(Vs, partition_clpfd(N, Vs), All), 35 length(All, P) 36 ; partition_number(N, P) 37 ). 38 39% ?- between(1, 52, N), partition_number(N, X, [history]), writeln(#('Y'(N)) = X), fail; true. 40 41% ?- partition_number(500, P, [formula]). 42% ?- partition_number(500, P, []). 43% ?- partition_number(40, P, []). 44% ?- partition_number(40, P, [formula]). 45% ?- partition_number(40, P, [rec]). 46% ?- partition_number(40, P, [history]). 47% ?- partition_number(40, P, [history2]). 48% ?- partition_number(20, P, [history]). 49% ?- partition_number(20, P, [clpfd]). 50% ?- time(partition_number(52, X, [formula])). 51% ?- time(partition_number(52, X, [history])). 52%@ % 47,461,267 inferences, 5.356 CPU in 5.489 seconds (98% CPU, 8861529 Lips) 53%@ X = 281589 . 54% ?- time(partition_number(52, X, [history2])). 55%@ % 47,461,222 inferences, 5.722 CPU in 5.833 seconds (98% CPU, 8295086 Lips) 56%@ X = 281589 . 57% ?- time(partition_number(53, X, [history])). 58%@ % 56,645,575 inferences, 6.636 CPU in 6.749 seconds (98% CPU, 8536582 Lips) 59%@ X = 329931 . 60% ?- time(partition_number(53, X, [history2])). 61%@ % 56,645,538 inferences, 6.478 CPU in 6.598 seconds (98% CPU, 8743664 Lips) 62% ?- time(partition_number(54, X, [history2])). 63%@ % 67,495,266 inferences, 7.389 CPU in 7.528 seconds (98% CPU, 9135086 Lips) 64%@ X = 386155 . 65% ?- time(partition_number(55, X, [history2])). 66%@ % 80,293,565 inferences, 9.192 CPU in 9.345 seconds (98% CPU, 8734916 Lips) 67%@ X = 451276 . 68% ?- time(partition_number(60, X, [history2])). 69%@ % 187,006,232 inferences, 20.617 CPU in 20.900 seconds (99% CPU, 9070386 Lips) 70%@ X = 966467 . 71 72 73% ?- time(partition_number(52, X, [rec])). % Out of global stack 74% ?- time(partition_number(53, X, [])). 75% ?- time(partition_number(53, X, [history])). % max 76% ?- time(partition_number(53, X, [history2])). % Out of Global Stack 77 78% Conclusion: induction with `history' is a little bit 79% better than the simple recursion on the integer. 80% 81% ?- time(partition_number(53, X, [rec])). % out of stack. 82% ?- time(partition_number(54, X, [history])). % out of stack. 83 84 85 /****************************************************** 86 * ?- set_prolog_stack(global, limit(2*10**9)). * 87 ******************************************************/ 88 89% ?- time(partition_number(65, X, [history])). 90%@ X = 2012558 . 91% ?- time(partition_number(65, X, [history2])). 92%@ X = 2012558 . 93 94 /******************************************************************** 95 * (F) Fast partion number based on the regression formula: * 96 * * 97 * $\sum_{j\geq0} p(n-\frac{j(3j\pm 1)}2)(-1)^j = 0$ * 98 ********************************************************************/ 99 100% ?- partition_number(1, X, [formula]). 101% ?- partition_number(2, X, [formula]). 102% ?- partition_number(3, X, [formula]). 103% ?- partition_number(4, X, [formula]). 104% ?- partition_number(6, X, [formula]). 105% ?- partition_number(7, X, [formula]). 106% ?- partition_number(20, X, [formula]). 107% ?- time(partition_number(10000, X, [formula])). 108%@ % 9,075,313 inferences, 1.560 CPU in 1.578 seconds (99% CPU, 5816379 Lips) 109%@ X = 36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144. 110 111% 112partition_number(N, P):- succ(N, N1), 113 fast_partition(N1, V), 114 arg(N1, V, P). 115 116% ?- fast_partition(1, X). 117% ?- fast_partition(2, X). 118% ?- fast_partition(3, X). 119% ?- fast_partition(4, X). 120% ?- fast_partition(5, X). 121% ?- fast_partition(6, X). 122% ?- fast_partition(7, X). 123% ?- fast_partition(8, X). 124% ?- fast_partition(9, X). 125% ?- fast_partition(10, X). 126% ?- fast_partition(11, X). 127% ?- fast_partition(21, X). 128 129% 130fast_partition(N, Vector):- functor(Vector, p, N), 131 arg(1, Vector, 1), 132 fast_partition(2, N, Vector). 133 134% 135fast_partition(K, N, _):- K> N, !. 136fast_partition(K, N, V):- fast_partition(1, K, V, 0, S), 137 arg(K, V, S), 138 succ(K, K1), 139 fast_partition(K1, N, V). 140% 141fast_partition(J, K, V, S, S0):- 142 sign_index(J, Sign, X, Y), 143 K1 is K-X, 144 K2 is K-Y, 145 ( K1=<0, K2=<0 -> S0 = S 146 ; ( K1 > 0 -> arg(K1, V, A1); A1 = 0 ), 147 ( K2 > 0 -> arg(K2, V, A2); A2 = 0 ), 148 S1 is S + Sign * (A1 + A2), 149 succ(J, J1), 150 fast_partition(J1, K, V, S1, S0) 151 ). 152 153% 154sign_index(J, Sign, X, Y):- 155 ( J mod 2 =:= 0 -> Sign = -1; Sign = 1 ), 156 A is 3 * J * J, 157 X is (A - J) // 2, 158 Y is (A + J) // 2. 159 160% ?- time(partition_number_list(10000, X)). 161%@ % 1,462,095,308 inferences, 96.689 CPU in 96.820 seconds (100% CPU, 15121631 Lips) 162%@ X = 36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144. 163 164partition_number_list(N, P):- succ(N, N1), 165 fast_partition_list(N1, V), 166 nth1(N1, V, P). 167 168% ?- fast_partition_list(10, X). 169fast_partition_list(N, Vector):- length(Vector, N), 170 nth1(1, Vector, 1), 171 fast_partition_list(2, N, Vector). 172 173% 174fast_partition_list(K, N, _):- K> N, !. 175fast_partition_list(K, N, V):- fast_partition_list(1, K, V, 0, S), 176 nth1(K, V, S), 177 succ(K, K1), 178 fast_partition_list(K1, N, V). 179% 180fast_partition_list(J, K, V, S, S0):- 181 sign_index(J, Sign, X, Y), 182 K1 is K-X, 183 K2 is K-Y, 184 ( K1=<0, K2=<0 -> S0 = S 185 ; ( K1 > 0 -> nth1(K1, V, A1); A1 = 0 ), 186 ( K2 > 0 -> nth1(K2, V, A2); A2 = 0 ), 187 S1 is S + Sign * (A1 + A2), 188 succ(J, J1), 189 fast_partition_list(J1, K, V, S1, S0) 190 ). 191 192 /**************************************** 193 * integer partition by recursion * 194 ****************************************/ 195 196% ?- scm(int_part). 197% ?- partition_number(1, X, [rec]). 198% ?- partition_number(2, X, [rec]). 199% ?- partition_number(3, X, [rec]). 200% ?- partition_number(10, X, [rec]). 201 202integer_partition(1, [[1-1]]). 203integer_partition(N, Ps):- N>1, 204 succ(N0, N), 205 integer_partition(N0, Qs), 206 bump_shift_list(Qs, P0s, []), 207 sort(P0s, Ps). 208 209% 210reverse([], X, X). 211reverse([A|R], X, Y):- reverse(R, [A|X], Y). 212 213% 214bump_shift_list([], X, X). 215bump_shift_list([A|B], X, Y):- 216 bump_shift(A, X, X0), 217 bump_shift_list(B, X0, Y). 218 219% ?- bump_shift([1-2], X, Y). 220% ?- bump_shift([2-1, 3-1], X, Y). 221bump_shift([], X, X). 222bump_shift([1-K|U], [[1-K0|U]|X], Y):- !, succ(K, K0), 223 bump_shift([], [1-K|U], X, Y). 224bump_shift(U, [[1-1|U]|X], Y):- bump_shift([], U, X, Y). 225 226% ?- bump_shift([],[1-2], X, Y). 227%@ X = [[1-1, 2-1]|Y] . 228bump_shift(_Rev, [], X, X). 229bump_shift(Rev, [I-J|V], X, Y):- succ(I, I0), 230 succ(J0, J), 231 bump_shift(Rev, I, J0, I0, V, X, X0), 232 bump_shift([I-J|Rev], V, X0, Y). 233% 234bump_shift(Rev, _I, 0, I0, V, [U|X], X):- !, bump_shift_(I0, V, U0), 235 reverse(Rev, U0, U). 236bump_shift(Rev, I, J, I0, V, [U|X], X):- bump_shift_(I0, V, U0), 237 reverse(Rev, [I-J|U0], U). 238 239bump_shift_(I, [I-K|V], [I-K0|V]):- !, succ(K, K0). 240bump_shift_(I, V, [I-1|V]). 241 242 /************************************** 243 * Inversion of Polynomials. * 244 **************************************/ 245 246% For a given polynimial p(x) = ax + ... (a is non zero rational), 247% there is a unique polynomial q(x) of a specified degree 248% such that p(q(x) + r(x)) = x as a formal power series 249% for some power series r(x) whose lowest degree is 250% larger than the degree of p(x). q(x) is considered 251% as a polynomial approximation to the inverse of 252% of the polynomial function p(x) around zero. 253 254% This program is written as an application of the integer partition 255% just for fun. 256 257% A list [1, 0, 1, 0, 0, 0], for example, means 258% a polynimial x + x^3. 259 260% ?- scm(int_part). 261% ?- inverse_poly([3, 1], R, [verify]). 262% ?- inverse_poly([2, 1], R, [verify]). 263% ?- inverse_poly([1, 0, 0], R, [verify]). 264% ?- inverse_poly([1, 1, 0, 0, 0, 0], R,[verify]). 265%@ [1,0,0,0,0,0,-132,165,-308,616,-1176,1764,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] 266%@ R = [1, -1, 2, -5, 14, -42] . 267% ?- inverse_poly([1, 0, 1], R,[verify]). 268% ?- inverse_poly([1, 1, 0, 0, 0, 0], R,[verify]). 269% ?- inverse_poly([3, 0, 1, 0, 0, 0], R,[verify]). 270% ?- inverse_poly([1, 0, 1, 0, 0, 0], R,[verify]). 271% ?- inverse_poly([1, 1, 0, 0, 0, 0, 0, 0, 0, 0], R,[verify]). 272% ?- inverse_poly([1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], R,[verify]). 273% ?- inverse_poly([1,1], R, [degree(20), verify]). 274% ?- inverse_poly([1,1], R, [degree(30), verify]). 275% ?- inverse_poly([1,1], R, [degree(35), verify]). 276% ?- inverse_poly([1,1], R, [degree(35), verify]). 277% ?- inverse_poly([1,1], R, [degree(36), verify]). 278% ?- inverse_poly([1,1], R, [degree(37), history, verify]). 279% ?- inverse_poly([1,1], R, [degree(37), rec, verify]). 280%@ ERROR: Out of local stack 281%@ ERROR: Out of local stack 282 283% ?- inverse_poly([1,1], R, [degree(37), power_table, verify]). 284% ?- scm(int_part). 285% ?- inverse_poly([-2], R, []). 286% ?- inverse_poly([-2], R, []). 287% ?- inverse_poly([-2], R). 288% ?- inverse_poly([1, 0], R). 289% ?- inverse_poly([1, 1], R, []). 290% ?- inverse_poly([1, 1, 0], R, [verify]). 291% ?- inverse_poly_by_power_increment(120, [1, 1], R). 292% ?- inverse_poly_by_power_increment(121, [1, 1], R). 293% ?- inverse_poly_by_power_increment(120, [1, 1], R). 294% ?- inverse_poly([1, 1], R, [degree(50), power_table]). 295% ?- inverse_poly([1, 1], R, [degree(50), rec]). 296% ?- inverse_poly([1, 1], R, [degree(50)]). 297% ?- inverse_poly([1, 1], R, [degree(120), power_table]). 298%@ R = [1, -1, 2, -5, 14, -42, 132, -429, 1430|...] . 299 300% ?- inverse_poly([1, 1], R, [degree(121), power_table]). 301% ?- inverse_poly([1, 1], R, [degree(122), power_table]). 302%@ R = [1, -1, 2, -5, 14, -42, 132, -429, 1430|...] . 303 304% ?- inverse_poly([1, 1], R, [degree(123), power_table]). 305%@ ERROR: Out of local stack 306% ?- inverse_poly([1, 1], R, [degree(122), power_table, verify]). 307 308% ?- inverse_poly([1, 1], R, [degree(11), rec]). 309%@ R = [1, -1, 2, -5, 14, -42, 132, -429, 1430|...] . 310% 311inverse_poly(Poly, R, Opt0):- 312 ( memberchk(degree(Deg), Opt0) 313 -> Opt = Opt0 314 ; length(Poly, Deg), 315 Opt = [degree(Deg)|Opt0] 316 ), 317 length(Poly, Lp), 318 N is max(Deg, Lp), 319 length(Poly0, N), 320 append(Poly, Q, Poly0), 321 maplist(=(0), Q), 322 Vec=..[u|Poly0], 323 length(R, N), 324 IVec=..[v|R], 325 Poly=[C1|_], 326 D1 is rdiv(1, C1), 327 R = [D1|_], 328 ( memberchk(rec, Opt) 329 -> integer_partition(2, Ps), 330 inverse_poly(2, N, Vec, IVec, Ps) 331 ; course_of_partition(2, Hs), 332 inverse_poly_history(2, N, Vec, IVec, Hs) 333 ), 334 ( memberchk(verify, Opt0) 335 -> math:poly(##([0|Poly], [0|R]), [_|U]), 336 % polynimial composition. 337 writeln(U) 338 ; true 339 ). 340 341% ?- poly_to_list(poly(3, [1]), R). 342%@ R = [0, 0, 0, 1]. 343% ?- poly_to_list(poly(3, [0]), R). 344poly_to_list(poly(0, L), L):-!. 345poly_to_list(poly(E, L), M):- !, succ(E0, E), 346 poly_to_list(poly(E0, [0|L]), M). 347poly_to_list(L, L). 348 349% 350inverse_poly(I, N, _, _, _):- I > N, !. 351inverse_poly(I, N, Vec, IVec, Ps):- 352 nth_coefficient(I, Vec, IVec, C, Ps), 353 arg(I, IVec, B), 354 arg(1, Vec, A1), 355 B is - (C rdiv A1), 356 succ(I, I0), 357 bump_shift_list(Ps, Q0s, []), 358 sort(Q0s, Qs), 359 inverse_poly(I0, N, Vec, IVec, Qs). 360 361% 362inverse_poly_history(I, N, _, _, _):- I > N, !. 363inverse_poly_history(I, N, Vec, IVec, Hs):- 364 Hs=[Ps|_], 365 nth_coefficient(I, Vec, IVec, C, Ps), 366 arg(I, IVec, B), 367 arg(1, Vec, A1), 368 B is - (C rdiv A1), 369 succ(I, I0), 370 extend_partition_list(Hs, I0, Qs, []), 371 sort(Qs, Q0s), 372 inverse_poly_history(I0, N, Vec, IVec, [Q0s|Hs]). 373 374% ?- fold_partition([1-1, 2-2], 6, v(5, 6, 1), V). 375%@ V = 540 376% ?- fold_partition([3-1, 2-2], 120, v(5, 6, 1), V). 377%@ V = 2160 . 378% ?- X is (120*(5**3)*(6**2)) // (6*2). 379%@ X = 45000. 380 381fold_partition(P, Jfac, Cs, D):- 382 foldl(pred(Cs, ([N-K, U, V] 383 :- arg(N, Cs, B), 384 math:scalar_exp(B, K, P0), 385 V is P0 * U)), 386 P, 1, D0), 387 foldl(pred(([_-K, U, V] 388 :- factorial(K, Fk), 389 V is Fk*U)), 390 P, 1, F), 391 D is (Jfac rdiv F)*D0. 392 393% ?- fold_partition_list([[1-1, 2-2],[1-1, 2-2]], 3, v(5, 6), V). 394%@ V = 1080 . 395 396fold_partition_list(Ps, J, Cs, E):- 397 factorial(J, Jfac), 398 foldl(pred([Cs, Jfac], ( [P, U, V] 399 :- fold_partition(P, Jfac, Cs, A), 400 V is A+U)), 401 Ps, 0, E). 402 403% 404nth_coefficient(I, U, V, C, Ps):- 405 nth_coefficient(I, 2, Ps, U, V, 0, C). 406% 407nth_coefficient(I, J, _, _, _, C, C):- J > I, !. 408nth_coefficient(I, J, Ps, U, V, C, C0):- 409 collect_by_partition_size(Ps, J, Qs), 410 fold_partition_list(Qs, J, V, A), 411 arg(J, U, M), 412 C1 is C + M*A, 413 succ(J, J0), 414 nth_coefficient(I, J0, Ps, U, V, C1, C0). 415 416% 417collect_by_partition_size([A|R], K, [A0|R0]):- 418 ( A = p(_, A0); A =A0 ), 419 partition_size(A0, K), 420 !, 421 collect_by_partition_size(R, K, R0). 422collect_by_partition_size([_|R], K, R0):- !, 423 collect_by_partition_size(R, K, R0). 424collect_by_partition_size([], _, []). 425 426% 427partition_size([], 0). 428partition_size([_-I|R], K):- K >= 0, K0 is K -I, 429 partition_size(R, K0). 430 431% 432factorial(N, F):- factorial(N, 1, F). 433% 434factorial(1, N, N). 435factorial(K, N, N0):- N1 is K*N, succ(K0, K), 436 factorial(K0, N1, N0). 437% 438mul_list(L, M):- foldl(pred([X,Y,Z]:- Z is X*Y), L, 1, M). 439 440% ?- num_of_bags(4, [1,2,3], V). 441 442num_of_bags(N, L, V):- factorial(N, F), 443 mul_list(L, M), 444 V is F//M. 445 446 447 /******************************************************** 448 * inverse a polynomial by incrementing exponents * 449 ********************************************************/ 450 451% [2016/04/20] Drastic improvent without integer partitions !! 452% A drastic improvent for inversing polynomials 453% by decreasing multiplications on polynomials. 454% It is based on a familiar and well-known binomial formula: 455% 456% (a + b)^n = Sum of nCr a^(n-r)b^r 457% 458% with i running from 0 to n. 459% 460% For examples: 461% (a + b)^2 = a^2 +2ab + b^2 462% (a + b)^3 = a^3 +3a^2b + 3ab^2 + b^3 463 464% combinations 465% ?- combination(4, 0, X). 466% ?- combination(4, 1, X). 467% ?- combination(4, 2, X). 468% ?- combination(4, 3, X). 469% ?- combination(4, 4, X). 470 471% ?- scm(int_part). 472% ?- poly_coeff(1, poly(1, [2]), V). 473%@ V = 2. 474% ?- poly_coeff(1, poly(1, []), V). 475%@ V = 0. 476% ?- poly_coeff(1, [], V). 477%@ V = 0. 478% ?- poly_coeff(1, [1,2], V). 479%@ V = 2. 480 481poly_coeff(I, poly(J, L), V):- !, 482 plus(K, J, I), 483 ( nth0(K, L, V) 484 -> true 485 ; V=0 486 ). 487poly_coeff(I, L, V):- 488 ( nth0(I, L, V) 489 -> true 490 ; V=0 491 ). 492 493% ?- scm(int_part). 494% ?- fold_poly_coeff(1,[3,5], [poly(1, [1]), poly(1, [2])], X). 495%@ X = 13 . 496% ?- fold_poly_coeff(1,[3], [poly(1, [1]), poly(1, [2])], X). 497%@ X = 3 . 498% ?- fold_poly_coeff(1,[3,5], [poly(1, [1])], X). 499%@ X = 3. 500 501fold_poly_coeff(I, Ws, Ps, X):- fold_poly_coeff(I, Ws, Ps, 0, X). 502% 503fold_poly_coeff(I, [W|Ws], [P|Ps], X, Y):-!, 504 poly_coeff(I, P, C), 505 X0 is W*C+X, 506 fold_poly_coeff(I, Ws, Ps, X0, Y). 507fold_poly_coeff(_, [], _, X, X). 508fold_poly_coeff(_, _, [], X, X). 509 510% ?- scm(int_part). 511% ?- inverse_poly_by_power_increment(3, [1], R). 512% ?- inverse_poly([1], R, [degree(3), power_table]). 513% ?- inverse_poly([1, 1], R, [degree(10), power_table]). 514% ?- inverse_poly([1, 1], R, [degree(20), power_table]). 515% ?- inverse_poly([1, 1], R, [degree(100), power_table]). 516% ?- inverse_poly([1, 1], R, [degree(120), power_table]). 517% ?- inverse_poly([1, 1], R, [degree(10), power_table,verify]). 518% ?- inverse_poly([1, 1], R, [degree(20), power_table]). 519% ?- inverse_poly([1, 1], R, [degree(100), power_table]). 520% ?- inverse_poly([1, 1], R, [degree(120), power_table]). 521% ?- inverse_poly([1, 1], R, [degree(121), power_table]). 522%@ R = [1, -1, 2, -5, 14, -42, 132, -429, 1430|...] . 523 524% ?- set_prolog_stack(global, limit(2*10**9)). 525%@ true. 526 527inverse_poly_by_power_increment(N, [A1|As], R):- 528 B1 is 1 rdiv A1, 529 update_power(1, poly(1, [B1]), [[]], N, Ps), 530 repeat_update_power(2, N, A1, As, Ps, Qs), 531 last(Qs, R0), 532 once(trim_constant_term(R0, R)). 533 534% For compatibility with other methods 535trim_constant_term(poly(0, [_|R]), R). 536trim_constant_term(poly(1, R), R). 537trim_constant_term(poly(J, R), S):- succ(J0, J), 538 trim_constant_term(poly(J0, [0|R]), S). 539trim_constant_term([_|R], R). 540 541% 542repeat_update_power(I, N, _, _, Ps, Ps):- I > N, !. 543repeat_update_power(I, N, A, Ws, Ps, Qs):- 544 reverse(Ps, [_|Us]), 545 fold_poly_coeff(I, Ws, Us, B), 546 C is - (B rdiv A), 547 succ(I, I0), 548 update_power(I, poly(I, [C]), Ps, N, Rs), 549 repeat_update_power(I0, N, A, Ws, Rs, Qs). 550 551% ?- update_power(1, poly(1,[1]), [[]], 3, X). 552%@ X = [poly(0, [0, 0, 1]), poly(0, [0, 1])] . 553% ?- update_power(2, poly(3,[1]), [poly(2,[1]), poly(1,[1])], 10, P). 554%@ P = [poly(0, [0, 0, 0, 1, 0, 3, 0, 3, 0, 1]), poly(0, [0, 0, 1, 0, 2, 0, 1]), poly(0, [0, 1, 0, 1])] . 555 556update_power(N, A, Ps, Max, [QH|Qs]):- 557 update_power(N, A, Ps, Max, Qs, []), 558 last(Qs, Q), 559 Qs = [H|_], 560 math:poly_mul_bd_poly(Q, H, Max, QH). 561% 562update_power(0, _, _, _, X, X):-!. 563update_power(I, A, Ps, Max, [P0|X], Y):- 564 fold_poly_power_combi(I, A, Ps, Max, P0), 565 succ(I0, I), 566 Ps=[_|Qs], 567 update_power(I0, A, Qs, Max, X, Y). 568 569% ?- fold_poly_power_combi(2, poly(3,[1]), [poly(2,[1]), poly(1,[1])], 10, P). 570fold_poly_power_combi(N, A, Ps, Max, P):- 571 fold_poly_power_combi(Ps, 0, N, A, poly(0,[1]), Max, [], P). 572 573% 574fold_poly_power_combi([], _, _, _, V, _, P, Q):- !, 575 math:poly_add(V, P, Q). 576fold_poly_power_combi([U|Us], I, N, A, V, Max, P, Q):- 577 combination(N, I, C), 578 math:poly_mul_bd_poly(V, U, Max, VU), 579 math:poly_scalar(C, VU, CVU), 580 math:poly_add(CVU, P, P0), 581 succ(I, I0), 582 math:poly_mul_bd_poly(A, V, Max, AV), 583 fold_poly_power_combi(Us, I0, N, A, AV, Max, P0, Q). 584 585% combi(N, R, C) <=> 586% Let X, R be a set s.t. |X| = N, a number R, respectively, 587% then C = |{Y in powerset(X) | |Y| = R}| 588 589combination(_, 0, 1):- !. 590combination(N, N, 1):- !. 591combination(N, R, C):- N >= R, 592 NR is N - R, 593 ( R < NR 594 -> combination(N, R, NR, C) 595 ; combination(N, NR, R,C) 596 ). 597 598% 599combination(N, R, S, C):- succ(S, S1), 600 nat_prod(S1, N, 1, M), 601 factorial(R, FR), 602 C is M rdiv FR. 603 604% 605nat_prod(K, N, X, X):- K>N, !. 606nat_prod(K, N, X, Y):- X0 is K*X, 607 succ(K, K1), 608 nat_prod(K1, N, X0, Y). 609 610% A bench mark of two methods (A) and (B) on integer partitions. 611 612 /************************************************************ 613 * (A) integer partition by course-of-value induction * 614 ************************************************************/ 615 616% ?- partition_number(40, X, [history]). 617%@ X = 37338 . 618% ?- partition(10, X), maplist(writeln, X). 619 620partition(N, Ps):- course_of_partition(N, [Ps|_]). 621 622% ?- course_of_partition(4, X). 623 624incremental_partition(Ps, N, [Qs|Ps]):- 625 extend_partition_list(Ps, N, Q0s, []), 626 sort(Q0s, Qs). 627 628% ?- scm(int_part). 629% ?- course_of_partition(2, X). 630% ?- course_of_partition(3, X). 631% ?- course_of_partition(4, X). 632% ?- course_of_partition(1,X), incremental_partition(X,2,Y). 633% ?- course_of_partition(1,X), incremental_partition(X,2,Y), incremental_partition(Y,3,Z). 634 635course_of_partition(0,[[p(0, [])]]). 636course_of_partition(N, Ps) :- N > 0, 637 succ(N1, N), 638 course_of_partition(N1, Qs), 639 incremental_partition(Qs, N, Ps). 640 641% ?- scm(int_part). 642extend_partition_list([], _, P, P). 643extend_partition_list([A|R], N, P, Q):- 644 extend_partition_level(A, N, P, P0), 645 extend_partition_list(R, N, P0, Q). 646 647% 648extend_partition_level([], _, P, P). 649extend_partition_level([X|R], N, P, Q):- 650 extend_canonical(X, N, P, P0), 651 !, 652 extend_partition_level(R, N, P0, Q). 653 654% Increasing order on block size 655% (decreasing order on block size does not work. 656extend_canonical(p(_, []), N, [p(N, [N-1])|D], D). 657extend_canonical(p(I, [Block_length-K|R]), N, [P|D], D):- 658 N1 is N-I, 659 ( N1 > Block_length 660 -> P = p(N, [N1-1, Block_length-K|R]) 661 ; N1 == Block_length 662 -> succ(K, K1), 663 P = p(N, [Block_length-K1|R]) 664 ). 665extend_canonical(_, _, D, D). 666 667 668 /****************************** 669 * partition_history_tr * 670 ******************************/ 671% Tail recursion version. 672% Despite of tail-recursion, it shows weaker performance ! Why ? 673 674% ?- scm(int_part). 675%@ true. 676% ?- partition_history_tr(0, X). 677% ?- partition_history_tr(1, X). 678% ?- partition_history_tr(2, X). 679% ?- partition_history_tr(3, X). 680% ?- partition_history_tr(4, X). 681% ?- partition_history_tr(40, X). 682% ?- partition_history_tr(50, X). 683% ?- partition_history_tr(52, X). % stack overflow 684% ?- partition_history_tr(53, X). % stack overflow 685% ?- partition_history_tr(52, X), flatten(X, Y), length(Y, L). 686% ?- partition_history_tr(40, X), flatten(X, Y), length(Y, L). 687% ?- partition_number(40, N, [formula]). 688%@ N = 37338. 689 690partition_history_tr(0, [[p(0, [])]]). 691partition_history_tr(N, Ps) :- N > 0, 692 partition_history_tr(0, P0s), 693 partition_history_tr(1, N, P0s, Ps). 694 695% 696partition_history_tr(I, N, Ps, Ps):- I > N, !. 697partition_history_tr(I, N, Ps, Qs):- 698 extend_partition_list(Ps, I, A, []), 699 sort(A, B), 700 succ(I, J), 701 partition_history_tr(J, N, [B|Ps], Qs). 702 703 /************************************************* 704 * (B) integer partition in Finite Domain * 705 *************************************************/ 706 707% ?- time(partition_clpfd_count(20, X)). 708% ?- time(partition_clpfd_count(30, X)). 709% ?- time(partition_clpfd_count(40, X)). 710% ?- time(partition_clpfd_count(45, X)). 711% ?- time(partition_clpfd_count(50, X)). 712%@ % 2,800,717,455 inferences, 222.620 CPU in 223.016 seconds (100% CPU, 12580703 Lips) 713%@ X = 204226. 714% ?- partition_clpfd_count(51, X). 715%@ ERROR: Out of global stack 716 717% ?- partition_clpfd(10, X). 718partition_clpfd_count(N, N0):- setof(Vs, partition_clpfd(N, Vs), All), 719 length(All, N0). 720 721% 722partition_clpfd(N, Vs):- 723 length(Vs, N), 724 ins(Vs, ..(0, N)), 725 numlist(1, N, Ns), 726 inner_product(Vs, Ns, 0, N), 727 label(Vs). 728% 729inner_product([], [], A, A). 730inner_product([X|R], [I|S], A, B):- 731 #=(I*X + A, A0), 732 inner_product(R, S, A0, B). 733 734 /*************************************** 735 * (C) naive, but inefficient version * 736 ***************************************/ 737 738% ?- naive_young(2, Ps). 739% ?- naive_young(3, Ps). 740% ?- naive_young(4, Ps). 741% ?- naive_young(5, Ps). 742% ?- time(naive_young(6, Ps)). 743%@ % 3,098,266 inferences, 0.328 CPU in 0.367 seconds (89% CPU, 9441874 Lips) 744%@ Ps = [ ... ]. 745 746% ?- time(naive_young(7, Ps)). 747%@ % 42,369,834 inferences, 3.568 CPU in 3.833 seconds (93% CPU, 11874930 Lips) 748%@ ERROR: Out of local stack 749 750naive_young(N, Ps):- n_powered_n(N, Series), 751 collect(pred(N, [V] :- 752 fold_vector(V, N)), 753 Series, Ps). 754 755% ?- n_powered_n(8, P). 756% ?- n_powered_n(2, P). 757% ?- n_powered_n(4, P). 758 759% ?- util:choices([[1,2],[1,2]], R). 760n_powered_n(N, P):- numlist(0, N, L), 761 length(Q, N), 762 maplist(=(L), Q), 763 choices(Q, P). 764 765% ?- fold_vector([1,2,3], S). 766%@ S = 14. 767fold_vector(R, S):- fold_vector(R, 1, 0, S). 768 769% 770fold_vector([], _, A, A). 771fold_vector([X|R], J, A, B):- 772 C is J*X + A, 773 succ(J, J0), 774 fold_vector(R, J0, C, B)