Extendable array for formal power series as array of virtually infinite size

I wrote codes for extendable array with two interface predicates poly_new/1 and index0/3, which are just simple wrappers of existing predicates in zdd library library(pac(prolog/zdd)). All sample queries in this post works after pack_installing pac and loading attached codes below.

Example of formal power series (Wikipedia) are

  1. 1
  2. 1 - x
  3. 1+x+x^2+x^3+\cdots

\Sigma_{n\geq0}a_nx^n= a_0 + a_1x+a_2x^2+\cdots is an FPS (formal power series) if a_i are numbers. That is, a formal power series is just a function which assigns a number to an integer n\geq 0. Basic operations on FPS are straightforward. With
p(x)=\Sigma_{n\geq 0}a_nx^n=a_0 + a_1x+\cdots+a_nx^n+\cdots, and
q(x)=\Sigma_{n\geq 0}b_nx^n,

  1. (add) p(x)+q(x)=\Sigma_{n\geq 0} (a_n+b_n) x^n.
  2. (scalar) For a number c, cp(x)=\Sigma_{n\geq 0}ca_nx^n.
  3. (multiply) p(x)q(x)=\Sigma_{n\geq 0}c_nx^n, where c_n= \Sigma_{i,j\geq0,i+j=n}a_ib_j (n\geq0)

Let F[[x]] be the set of FPSs over a certain fixed number field F.
For p(x)\in F[[x]], if q(x)\in F[[x]] satisfies p(x)q(x)=1 then q(x) is called a (multiplicative) inverse of p(x). Clearly an inverse of p(x) exists if and only if the constant term p(0)\ne 0. Also it is clear that the inverse exists uniquely, written p(x)^{-1}.
p(x)p(x)^{-1}= p(x)^{-1}p(x)=1.

It is easy to see that (1-x)^{-1} = 1+x+x^2+x^3+\cdots., usually written
\frac 1{1-x}= 1+x+x^2+x^3+\cdots.
As second example in this post, with p_n being the partition number of n\geq0, and P(x)= p_0 + p_1x+p_2x^2+\cdots= 1+1x+2x^2+3x^3+5x^4+7x^5+11x^6+15x^7+22x^8+30x^9+42x^{10}+56x^{11}+\cdots
Then what is P(x)^{-1} upto degree d ? (d=100 for example)

To compute coefficients upto degree 100 of inverse of P(x)
% ?- poly_new([P, Q]),  poly_inverse(100, pn-P, Q, V), writeln(Q).
 a(#(1,-1,-1,0,0,1,0,1,0,0,0,0,-1,0,0,-1,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,-1,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,_286,_288,_290,_292,_294,_296,_298,_300,_302,_304,_306,_308,_310,_312,_314,_316,_318,_320,_322,_324,_326,_328,_330,_332,_334,_336,_338))

Once I wrote prolog codes for FPS using lists (I was not aware of extendable array at that time). At least for the moment, I like the current version. One of the reason is that it is easy by take array as table to compute only once for the same input (index).
Stuff of this post may be a kind of folklore for this forum, but I would like to hear any comment. Thanks in advance.

Source code with sample queries
:- module(fpseries, [
			open_array/1
		,	open_array/2
		,	close_poly/1
		,	poly_new/1]).

:- use_module(zdd('zdd-array')).
:- use_module(pac(basic)).
:- use_module(pac('expand-pac')).
:- use_module(pac(op)).

term_expansion --> pac:expand_pac.

% ?- open_array(F), close_poly(F).
open_array(S):- open_array(S, 1).
%
open_array(a(Vector), VsizeExp):-
	Vsize is VsizeExp,
	functor(Vector, #, Vsize).
%
close_poly(S):- setarg(1, S, []).

%% index0(I, S, V)
%	S is an array of virtually infinete lengthl. V is a vlue at Ith element
%	of S. (I>=0.)
%
% ?-  poly_new(P), writeln(P), index0(3, P, V), writeln(P).
% ?-  poly_new(S), writeln(S), index(5, S, V), writeln(S).

index0(I, S, V):-  I0 is I + 1, index(I0, S, V).

%% poly_new(X)
%	open a new set of virtual formal power series for  variables
%	specified in X.

% ?- poly_new(X), writeln(X).
% ?- poly_new([X,Y]), writeln(X), writeln(Y).
poly_new(X):- var(X), !, poly_new([X]).
poly_new([]):-!.
poly_new([P|Ps]):- open_array(P), poly_new(Ps).

% ?- poly_init([],P).
% ?- poly_init([1,2,3],P), writeln(P).
poly_init(X, P):- open_array(P),
	poly_unify_args(X, 1, P).
%
poly_unify_args([], _, _).
poly_unify_args([A|As], I, P):-
	index(I, P, A),
	J is I+1,
	poly_unify_args(As, J, P).

% ?- poly_new(S), poly_factorial(5, F, S), writeln(S).
poly_factorial(0, 1, _):-!.
poly_factorial(N, F, S):- index(N, S, F),
	(	nonvar(F) -> true
	;	N0 is N - 1,
		poly_factorial(N0, F0, S),
		F is N*F0
	).

% ?-  open_array(S, 2), fibo(3, F, S).
% ?-  open_array(S, 2), fibo(1000, F, S).
% ?-  time((open_array(S, 2), fibo(10000, V, S))).
% ?-  time((poly_new(S), fibo(10000, V, S))).
% ?-  time((open_array(S, 2), fibo(10000, V, S), fibo(10000, V0, S))).
%@ % 110,052 inferences, 0.026 CPU in 0.033 seconds (79% CPU, 4191659 Lips)
%@ S = a(#(..)),
% ?-  time((open_array(S, 2), fibo(10000, F, S), index(10000, S, U))).
%@ % 110,051 inferences, 0.016 CPU in 0.023 seconds (70% CPU, 6894130 Lips)
% ?-  time((poly_new(S), fibo(10000, F, S))).
%@ % 130,043 inferences, 0.016 CPU in 0.022 seconds (74% CPU, 8082225 Lips)
%@ S = a(#(..)),
% ?- open_array(F), fibo(3, V, F), writeln(F).
fibo(0, 0, _).
fibo(1, 1, _).
fibo(N, F, S):- N>1, index0(N, S, F),
	(	nonvar(F) -> true
	;	N0 is N - 1,
		N1 is N - 1,
		fibo(N0, F0, S),
		fibo(N1, F1, S),
		F is F0 + F1
	).

% const(2), for example, is  2 + 2x+ 2x^2 + 2x^3 + 2x^4 + ...
const(K, _, K).
const(K, _, _, K).

% nat for 1+2x+3x^2 + 4x^3 + 5x^4 + ...
nat(I, J):- J is I+1.

%  p0 for 1-x
p0(0, 1).
p0(1, -1).
p0(J, 0):- J>1.

%
unit(0, 1).
unit(I, 0):- I>0.

% ?- poly_new(P), poly(pn, 0, P, V), writeln(P).
% ?- poly_new(P), poly(pn, 1, P, V), writeln(P).
% ?- poly_new(P), poly(pn, 100, P, V).
% ?- poly_new(P), poly(pn, 1000, P, V).
% ?- poly_new(P), poly(pn, 10000, P, V).
% ?- time((poly_new(P), poly(pn, 100000, P, V))).
% ?- poly_new(P), poly(#(const(3)), 5, P, V), writeln(P).
% ?- poly_new(P), poly(#(3), 5, P, V), writeln(P).
% ?- poly_new(P), poly(3, 5, P, V), writeln(P).

poly(#(X), I, _, V):-!,
	(	number(X) -> V = X
	;	call(X, I, V)
	).
poly(X, I, P, V):-
	(	number(X)  -> F = const(X)
	;	F = X
	),
	poly_rec(F, I, P, V).
%
poly_rec(F, I, P, V):- index0(I, P, V),
	(	nonvar(V) -> true
	;	(	I = 0 -> true
		;   I0 is I - 1,
			poly_rec(F, I0, P, _)
		),
		call(F, I, P, V)
	).

% ?- poly_new([H, A, B]), poly_add(100, poly(#1)-A, poly(1)-B, U, H),
% writeln(A), writeln(B), writeln(H).
%@ a(#(_1656))
% ?- poly_new([H, A, B]), poly_add(100, poly(#1)-A, poly(#1)-B, U, H),
% writeln(A), writeln(B), writeln(H).

poly_add(N, F-S, G-T, W, H):- index0(N, H, W),
	(	nonvar(W) -> true
	; 	(	N = 0 -> true
		;	N0 is N-1,
			poly_add(N0, F-S, G-T, _, H)
		),
		call(F, N, S, U),
		call(G, N, T, V),
		W is U+V
	).

% ?- poly_new([H, A, B]),  poly_subtr(10, poly(#1)-A, poly(#1)-B, U, H).
poly_subtr(N, F-S, G-T, W, H):- index0(N, H, W),
	(	nonvar(W) -> true
	;	(	N=0 -> true
		;	N0 is N -1,
			poly_subtr(N0, F-S, G-T, _, H)
		),
		call(F, N, S, U),
		call(G, N, T, V),
		W is U-V
	).

% Multiply a scalar for a polynomial.
% ?- poly_new([H, A]), poly_scalar(10, 3,  poly(#1)-A, U, H), writeln(H).
% ?- poly_new([H, A, B]), poly_scalar(10, 3, poly(1)-A, U, H),
%  poly_minus(10, index0-H, V, B), maplist(writeln, [A, H, B]).
poly_scalar(N, A, F-S, V, H):- index0(N, H, V),
	(	nonvar(V) -> true
	;	(	N = 0 -> true
		;	N0 is N- 1,
			poly_scalar(N0, A, F-S, _, H)
		),
		call(F, N, S, U),
		V is A*U
	).

% ?- poly_new([H, A]), poly_minus(10, poly(#1)-A, U, H), writeln(H).
% ?- poly_new([H, A, B]), poly_minus(10, poly(1)-A, U, H),
%  poly_minus(10, index0-H, V, B), maplist(writeln, [A, H, B]).
poly_minus(N, FS, V, H):- poly_scalar(N, -1, FS, V, H).


% ?- poly_new([H, A, B]), poly_mul(0, poly(1)-A, poly(1)-B, V, H), writeln(H).
% ?- poly_new([H, A, B]), poly_mul(1, poly(1)-A, poly(1)-B, V, H), writeln(H).
% ?- poly_new([H, A, B]), poly_mul(100, poly(1)-A, poly(1)-B, V, H), writeln(H).
% ?- poly_new([H, A, B]), poly_mul(100, poly(p0)-A, poly(p0)-B, V, H), writeln(H).
% ?- poly_new([H, A, B]), poly_mul(100, poly(p0)-A, poly(1)-B, V, H), writeln(H).
% ?- poly_new([H]), poly_mul(100, poly(#1)-_, poly(#1)-_, V, H), writeln(H).

poly_mul(N, F-P, G-Q, V, H):- index0(N, H, V),
	(	nonvar(V) -> true
	;	(	N = 0 	->	true
		;   N0 is N - 1,
			poly_mul(N0, F-P, G-Q, _, H)
		),
		poly_convolute(N, F-P, G-Q, V)
	).

% ?- poly_new([H, A, B]), poly_convolute(2, poly(1)-A,  poly(1)-B, V).
% ?- poly_new([H, A, B]), poly_convolute(100, poly(1)-A, poly(1)-B, V).
poly_convolute(N, FP, GQ, V):-	poly_convolute(0, N, FP, GQ, 0, V).

%
poly_convolute(I, N, _, _, V, V):- I>N, !.
poly_convolute(I, N, F-P, G-Q, U, V):-
	call(F, I, P, A),
	J is N-I,
	call(G, J, Q, B),
	U0 is U + A*B,
	K is I + 1,
	poly_convolute(K, N, F-P, G-Q, U0, V).

%
poly_inverse_aux(I, N, _, _, V, V):- I>N, !.
poly_inverse_aux(I, N, F-P, Q, U, V):-
	call(F, I, P, A),
	J is N-I,
	index0(J, Q, B),
	U0 is U+A*B,
	I0 is I + 1,
	poly_inverse_aux(I0, N, F-P, Q, U0, V).

% ?- poly_new([P, Q]),  poly_inverse(0, poly(1)-P, Q, V).
% ?- poly_new([P, Q]),  poly_inverse(1, poly(1)-P, Q, V).
% ?- poly_new([P, Q]),  poly_inverse(10, poly(#p0)-P, Q, V), writeln(P), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(100, poly(#p0)-P, Q, V), writeln(P), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(1, poly(1)-P, Q, V), writeln(P), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(2, poly(1)-P, Q, V), writeln(P), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(1, poly(#nat)-P, Q, V).
% ?- poly_new([P, Q]),  poly_inverse(10, poly(#nat)-P, Q, V), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(10, fibo_plus_one-P, Q, V), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(100, fibo_plus_one-P, Q, V), writeln(Q).

fibo_plus_one(X, S, Y):-fibo(X, Y0, S), Y is Y0 + 1.

% ?- poly_new([P, Q]),  poly_inverse(1, poly(#p0)-P, Q, V), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(2, poly(#p0)-P, Q, V), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(10, poly(#p0)-P, Q, V), writeln(Q).
% ?- poly_new([P, Q]),  poly_inverse(100, pn-P, Q, V), writeln(P), writeln(Q).

poly_inverse(I, F-P, Q, V):- index0(I, Q, V),
	(	nonvar(V) ->  true
	;	I = 0
	->  call(F, 0, P, A0),
		V is rdiv(1, A0)
	;	I0 is I-1,
		poly_inverse(I0, F-P, Q, _),
		poly_inverse_aux(1, I, F-P, Q, 0, V0),
		call(F, 0, P, A0),
		V is -rdiv(V0, A0)
	).

		/**************************
		*    partition_number     *
		**************************/

% ?- time(pn(4)).
% ?- time(pn(10)).
% ?- time(pn(1000)).
% ?- time(pn(100000)).

pn(N):- poly_new(S),
	index0(0, S, 1),
	index0(1, S, 1),
	partition_number(N, _, S), !,
	(	between(0, N, I),
		index0(I, S, X),
		var(X),
		writeln(I),
		fail
	;	writeln("done.")
	).

% ?- time(pn(100, X)).
% ?- time(pn(10000, X)).
% ?- time(pn(100000, X)).

pn(N, V):- partition_number(N, V).
%
pn(N, S, P):-  partition_number(N, P, S).

% ?- check_inverse_pn(10).
% ?- check_inverse_pn(100).
% ?- time(check_inverse_pn(1000)).

check_inverse_pn(N):- poly_new([S, D, A]),
	index0(0, S, 1),
	index0(1, S, 1),
	partition_number(N, _, S),
	M is N//2,
	poly_inverse(M, index0-S, D, _),
	K is M//2,
	poly_mul(K, index0-S, index0-D, _, A),
	format("inverse check done.\n"),
	writeln(A).

% ?- parity(0, X).
parity(J, S):-  % S is power of -1 by J:  S = (-1)^ J.
	(	J mod 2 =:= 0  -> S = 1
	;	S = -1
	).

% ?- jm(1,X).
% ?- jp(1,X).
jm(J, X):- X is (J*(3*J-1)) //2.
jp(J, X):- X is (J*(3*J+1)) //2.

% ?- time(zdd((partition_number(100000, P)))).
%@ % 456,034,166 inferences, 39.034 CPU in 39.279 seconds (99% CPU, 11682862 Lips)
% ?- time(partition_number(100000, P)).
%@ % 456,033,880 inferences, 39.124 CPU in 39.368 seconds (99% CPU, 11656207 Lips)

partition_number(N, P):- poly_new(S),
	index0(0, S, 1),
	index0(1, S, 1),
	partition_number(N, P, S).

partition_number(N, 0, _):- N < 0, !.
partition_number(0, 1, S):-!, index0(0, S, 1).
partition_number(1, 1, S):-!, index0(1, S, 1).
partition_number(N, P, S):- index0(N, S, P),
	(	nonvar(P) -> true
	;	partition_number_convolute(N, 1, 0, P0, S),
		P is -P0
	).
%
partition_number_convolute(N, J, A, B, S):-
	jm(J, Jm),
	jp(J, Jp),
	Nm is N-Jm,
	Np is N-Jp,
	(	Nm < 0 -> B = A
	;	parity(J, P),
		partition_number(Nm, U0, S),
	    partition_number(Np, V0, S),
		A0 is A + P*(U0+V0),
		J0 is J + 1,
		partition_number_convolute(N, J0, A0, B, S)
	).

Kuniaki Mukai

2 Likes