1:- module(fpseries, [ 2 open_array/1 3 , open_array/2 4 , close_array/1 5 , poly_new/1]). 6 7:- use_module(zdd('zdd-array')). 8:- use_module(pac(basic)). 9:- use_module(pac('expand-pac')). 10:- use_module(pac(op)). 11 12term_expansion --> pac:expand_pac. 13 14% ?- open_array(F), close_array(F). 15 16open_array(S):- open_array(S, 1). 17% 18open_array(a(Vector), VsizeExp):- 19 Vsize is VsizeExp, 20 functor(Vector, #, Vsize). 21% 22close_array(S):- setarg(1, S, []).
?- poly_new(P), writeln(P), index0(3, P, V), writeln(P).
?- poly_new(S), writeln(S), index(5, S, V), writeln(S).
31index0(I, S, V):- I0 is I + 1, index(I0, S, V).37% ?- poly_new(X), writeln(X). 38% ?- poly_new([X,Y]), writeln(X), writeln(Y). 39poly_new(X):- var(X), !, poly_new([X]). 40poly_new([]):-!. 41poly_new([P|Ps]):- open_array(P), poly_new(Ps). 42 43% ?- poly_init([],P). 44% ?- poly_init([1,2,3],P), writeln(P). 45poly_init(X, P):- open_array(P), 46 poly_unify_args(X, 1, P). 47% 48poly_unify_args([], _, _). 49poly_unify_args([A|As], I, P):- 50 index(I, P, A), 51 J is I+1, 52 poly_unify_args(As, J, P). 53 54% ?- poly_new(S), poly_factorial(5, F, S), writeln(S). 55poly_factorial(0, 1, _):-!. 56poly_factorial(N, F, S):- index(N, S, F), 57 ( nonvar(F) -> true 58 ; N0 is N - 1, 59 poly_factorial(N0, F0, S), 60 F is N*F0 61 ). 62 63% ?- open_array(S, 2), fibo(3, F, S). 64% ?- open_array(S, 2), fibo(1000, F, S). 65% ?- time((open_array(S, 2), fibo(10000, V, S))). 66% ?- time((poly_new(S), fibo(10000, V, S))). 67% ?- time((open_array(S, 2), fibo(10000, V, S), fibo(10000, V0, S))). 68%@ % 110,052 inferences, 0.026 CPU in 0.033 seconds (79% CPU, 4191659 Lips) 69%@ S = a(#(..)), 70% ?- time((open_array(S, 2), fibo(10000, F, S), index(10000, S, U))). 71%@ % 110,051 inferences, 0.016 CPU in 0.023 seconds (70% CPU, 6894130 Lips) 72% ?- time((poly_new(S), fibo(10000, F, S))). 73%@ % 130,043 inferences, 0.016 CPU in 0.022 seconds (74% CPU, 8082225 Lips) 74%@ S = a(#(..)), 75% ?- open_array(F), fibo(3, V, F), writeln(F). 76fibo(0, 0, _). 77fibo(1, 1, _). 78fibo(N, F, S):- N>1, index0(N, S, F), 79 ( nonvar(F) -> true 80 ; N0 is N - 1, 81 N1 is N - 1, 82 fibo(N0, F0, S), 83 fibo(N1, F1, S), 84 F is F0 + F1 85 ). 86 87% const(2), for example, is 2 + 2x+ 2x^2 + 2x^3 + 2x^4 + ... 88const(K, _, K). 89const(K, _, _, K). 90 91% nat for 1+2x+3x^2 + 4x^3 + 5x^4 + ... 92nat(I, J):- J is I+1. 93 94% p0 for 1-x 95p0(0, 1). 96p0(1, -1). 97p0(J, 0):- J>1. 98 99% 100unit(0, 1). 101unit(I, 0):- I>0. 102 103% ?- poly_new(P), poly(pn, 0, P, V), writeln(P). 104% ?- poly_new(P), poly(pn, 1, P, V), writeln(P). 105% ?- poly_new(P), poly(pn, 100, P, V). 106% ?- poly_new(P), poly(pn, 1000, P, V). 107% ?- poly_new(P), poly(pn, 10000, P, V). 108% ?- time((poly_new(P), poly(pn, 100000, P, V))). 109% ?- poly_new(P), poly(#(const(3)), 5, P, V), writeln(P). 110% ?- poly_new(P), poly(#(3), 5, P, V), writeln(P). 111% ?- poly_new(P), poly(3, 5, P, V), writeln(P). 112 113poly(#(X), I, _, V):-!, 114 ( number(X) -> V = X 115 ; call(X, I, V) 116 ). 117poly(X, I, P, V):- 118 ( number(X) -> F = const(X) 119 ; F = X 120 ), 121 poly_rec(F, I, P, V). 122% 123poly_rec(F, I, P, V):- index0(I, P, V), 124 ( nonvar(V) -> true 125 ; ( I = 0 -> true 126 ; I0 is I - 1, 127 poly_rec(F, I0, P, _) 128 ), 129 call(F, I, P, V) 130 ). 131 132% ?- poly_new([H, A, B]), poly_add(100, poly(#1)-A, poly(1)-B, U, H), 133% writeln(A), writeln(B), writeln(H). 134%@ a(#(_1656)) 135% ?- poly_new([H, A, B]), poly_add(100, poly(#1)-A, poly(#1)-B, U, H), 136% writeln(A), writeln(B), writeln(H). 137 138poly_add(N, F-S, G-T, W, H):- index0(N, H, W), 139 ( nonvar(W) -> true 140 ; ( N = 0 -> true 141 ; N0 is N-1, 142 poly_add(N0, F-S, G-T, _, H) 143 ), 144 call(F, N, S, U), 145 call(G, N, T, V), 146 W is U+V 147 ). 148 149% ?- poly_new([H, A, B]), poly_subtr(10, poly(#1)-A, poly(#1)-B, U, H). 150poly_subtr(N, F-S, G-T, W, H):- index0(N, H, W), 151 ( nonvar(W) -> true 152 ; ( N=0 -> true 153 ; N0 is N -1, 154 poly_subtr(N0, F-S, G-T, _, H) 155 ), 156 call(F, N, S, U), 157 call(G, N, T, V), 158 W is U-V 159 ). 160 161% Multiply a scalar for a polynomial. 162% ?- poly_new([H, A]), poly_scalar(10, 3, poly(#1)-A, U, H), writeln(H). 163% ?- poly_new([H, A, B]), poly_scalar(10, 3, poly(1)-A, U, H), 164% poly_minus(10, index0-H, V, B), maplist(writeln, [A, H, B]). 165poly_scalar(N, A, F-S, V, H):- index0(N, H, V), 166 ( nonvar(V) -> true 167 ; ( N = 0 -> true 168 ; N0 is N- 1, 169 poly_scalar(N0, A, F-S, _, H) 170 ), 171 call(F, N, S, U), 172 V is A*U 173 ). 174 175% ?- poly_new([H, A]), poly_minus(10, poly(#1)-A, U, H), writeln(H). 176% ?- poly_new([H, A, B]), poly_minus(10, poly(1)-A, U, H), 177% poly_minus(10, index0-H, V, B), maplist(writeln, [A, H, B]). 178poly_minus(N, FS, V, H):- poly_scalar(N, -1, FS, V, H). 179 180 181% ?- poly_new([H, A, B]), poly_mul(0, poly(1)-A, poly(1)-B, V, H), writeln(H). 182% ?- poly_new([H, A, B]), poly_mul(1, poly(1)-A, poly(1)-B, V, H), writeln(H). 183% ?- poly_new([H, A, B]), poly_mul(100, poly(1)-A, poly(1)-B, V, H), writeln(H). 184% ?- poly_new([H, A, B]), poly_mul(100, poly(p0)-A, poly(p0)-B, V, H), writeln(H). 185% ?- poly_new([H, A, B]), poly_mul(100, poly(p0)-A, poly(1)-B, V, H), writeln(H). 186% ?- poly_new([H]), poly_mul(100, poly(#1)-_, poly(#1)-_, V, H), writeln(H). 187 188poly_mul(N, F-P, G-Q, V, H):- index0(N, H, V), 189 ( nonvar(V) -> true 190 ; ( N = 0 -> true 191 ; N0 is N - 1, 192 poly_mul(N0, F-P, G-Q, _, H) 193 ), 194 poly_convolute(N, F-P, G-Q, V) 195 ). 196 197% ?- poly_new([H, A, B]), poly_convolute(2, poly(1)-A, poly(1)-B, V). 198% ?- poly_new([H, A, B]), poly_convolute(100, poly(1)-A, poly(1)-B, V). 199poly_convolute(N, FP, GQ, V):- poly_convolute(0, N, FP, GQ, 0, V). 200 201% 202poly_convolute(I, N, _, _, V, V):- I>N, !. 203poly_convolute(I, N, F-P, G-Q, U, V):- 204 call(F, I, P, A), 205 J is N-I, 206 call(G, J, Q, B), 207 U0 is U + A*B, 208 K is I + 1, 209 poly_convolute(K, N, F-P, G-Q, U0, V). 210 211% 212poly_inverse_aux(I, N, _, _, V, V):- I>N, !. 213poly_inverse_aux(I, N, F-P, Q, U, V):- 214 call(F, I, P, A), 215 J is N-I, 216 index0(J, Q, B), 217 U0 is U+A*B, 218 I0 is I + 1, 219 poly_inverse_aux(I0, N, F-P, Q, U0, V). 220 221% ?- poly_new([P, Q]), poly_inverse(0, poly(1)-P, Q, V). 222% ?- poly_new([P, Q]), poly_inverse(1, poly(1)-P, Q, V). 223% ?- poly_new([P, Q]), poly_inverse(10, poly(#p0)-P, Q, V), writeln(P), writeln(Q). 224% ?- poly_new([P, Q]), poly_inverse(100, poly(#p0)-P, Q, V), writeln(P), writeln(Q). 225% ?- poly_new([P, Q]), poly_inverse(1, poly(1)-P, Q, V), writeln(P), writeln(Q). 226% ?- poly_new([P, Q]), poly_inverse(2, poly(1)-P, Q, V), writeln(P), writeln(Q). 227% ?- poly_new([P, Q]), poly_inverse(1, poly(#nat)-P, Q, V). 228% ?- poly_new([P, Q]), poly_inverse(10, poly(#nat)-P, Q, V), writeln(Q). 229% ?- poly_new([P, Q]), poly_inverse(10, fibo_plus_one-P, Q, V), writeln(Q). 230% ?- poly_new([P, Q]), poly_inverse(100, fibo_plus_one-P, Q, V), writeln(Q). 231 232fibo_plus_one(X, S, Y):-fibo(X, Y0, S), Y is Y0 + 1. 233 234% ?- poly_new([P, Q]), poly_inverse(1, poly(#p0)-P, Q, V), writeln(Q). 235% ?- poly_new([P, Q]), poly_inverse(2, poly(#p0)-P, Q, V), writeln(Q). 236% ?- poly_new([P, Q]), poly_inverse(10, poly(#p0)-P, Q, V), writeln(Q). 237% ?- poly_new([P, Q]), poly_inverse(100, pn-P, Q, V), writeln(P), writeln(Q). 238 239poly_inverse(I, F-P, Q, V):- index0(I, Q, V), 240 ( nonvar(V) -> true 241 ; I = 0 242 -> call(F, 0, P, A0), 243 V is rdiv(1, A0) 244 ; I0 is I-1, 245 poly_inverse(I0, F-P, Q, _), 246 poly_inverse_aux(1, I, F-P, Q, 0, V0), 247 call(F, 0, P, A0), 248 V is -rdiv(V0, A0) 249 ). 250 251 /************************** 252 * partition_number * 253 **************************/ 254 255% ?- time(pn(4)). 256% ?- time(pn(10)). 257% ?- time(pn(1000)). 258% ?- time(pn(100000)). 259 260pn(N):- poly_new(S), 261 index0(0, S, 1), 262 index0(1, S, 1), 263 partition_number(N, _, S), !, 264 ( between(0, N, I), 265 index0(I, S, X), 266 var(X), 267 writeln(I), 268 fail 269 ; writeln("done.") 270 ). 271 272% ?- time(pn(100, X)). 273% ?- time(pn(10000, X)). 274% ?- time(pn(100000, X)). 275 276pn(N, V):- partition_number(N, V). 277% 278pn(N, S, P):- partition_number(N, P, S). 279 280% ?- check_inverse_pn(10). 281% ?- check_inverse_pn(100). 282% ?- time(check_inverse_pn(1000)). 283 284check_inverse_pn(N):- poly_new([S, D, A]), 285 index0(0, S, 1), 286 index0(1, S, 1), 287 partition_number(N, _, S), 288 M is N//2, 289 poly_inverse(M, index0-S, D, _), 290 K is M//2, 291 poly_mul(K, index0-S, index0-D, _, A), 292 format("inverse check done.\n"), 293 writeln(A). 294 295% ?- parity(0, X). 296parity(J, S):- % S is power of -1 by J: S = (-1)^ J. 297 ( J mod 2 =:= 0 -> S = 1 298 ; S = -1 299 ). 300 301% ?- jm(1,X). 302% ?- jp(1,X). 303jm(J, X):- X is (J*(3*J-1)) //2. 304jp(J, X):- X is (J*(3*J+1)) //2. 305 306% ?- time(zdd((partition_number(100000, P)))). 307%@ % 456,034,166 inferences, 39.034 CPU in 39.279 seconds (99% CPU, 11682862 Lips) 308% ?- time(partition_number(100000, P)). 309%@ % 456,033,880 inferences, 39.124 CPU in 39.368 seconds (99% CPU, 11656207 Lips) 310 311partition_number(N, P):- poly_new(S), 312 index0(0, S, 1), 313 index0(1, S, 1), 314 partition_number(N, P, S). 315 316partition_number(N, 0, _):- N < 0, !. 317partition_number(0, 1, S):-!, index0(0, S, 1). 318partition_number(1, 1, S):-!, index0(1, S, 1). 319partition_number(N, P, S):- index0(N, S, P), 320 ( nonvar(P) -> true 321 ; partition_number_convolute(N, 1, 0, P0, S), 322 P is -P0 323 ). 324% 325partition_number_convolute(N, J, A, B, S):- 326 jm(J, Jm), 327 jp(J, Jp), 328 Nm is N-Jm, 329 Np is N-Jp, 330 ( Nm < 0 -> B = A 331 ; parity(J, P), 332 partition_number(Nm, U0, S), 333 partition_number(Np, V0, S), 334 A0 is A + P*(U0+V0), 335 J0 is J + 1, 336 partition_number_convolute(N, J0, A0, B, S) 337 )