-- File IterativeFormulaOSequences.CoCoA5
-- Last update: July 30, 2025

-- by Francesca Cioffi

-- AIM: 
-- Computation of the number of all O-sequences with a given multiplicity by an iterative formula

-- Algorithm: 
-- see Corollary 4.4 in the paper "Counting finite O-sequences of a given multiplicity" by F. Cioffi and M. Guida
 

define FOSequenze(P,N,K,D)
  Print "P=", P , " N=", N, " K=", K, " D=", D, NewLine();

  StartTime := CpuTime();

  if P<1 or K<0 or D<1 Or K>=D or K>N then
    Out:=0;
  else
    OS:=[[],[]];
    V:=[ [ [ 0 | Di In 1..D] | Ki In 1..(N+1)  ] | Ni In 1..(N+1) ];
    OS[1]:=V;
    Ind:=1;
    Ind_1:=2;

-- Assignment of data for the case Pi=1 (base induction)

    for Di:=1 to D do
      for Ni:=Di-1 to N do
        OS[1,Ni+1,Di,Di]:=1;
      endfor;
    endfor;

-- Recursion on P, D, N, K

-- Case P>1:

    for Pi:=2 To P do
      Print "Pi=", Pi, NewLine();

      App:=Ind;
      Ind:=Ind_1;
      Ind_1:=App;
      OS[Ind]:= V;

      for Di:=1 to D do
        For Ni:=0 To N Do
          OS[Ind,Ni+1,1,Di]:= Sum([OS[Ind_1,Ni+1,H+1,Di] | H In 0..Ni]); 
        endfor;
      endfor;

      for Ki:=1 to N do
        for Di:=Ki+1 to D do
          for Ni:=Ki to N do
            OS[Ind,Ni+1,Ki+1,Di]:=
 Sum([ Sum([ OS[Ind_1,Ni+1,I+1,Di-J]*OS[Ind,I,Ki,J] | I In Ki..N ]) | J In 1..(Di-1)]);
          endfor; 
        endfor;
      endfor;
    endfor;
    Out:=OS[Ind,N+1,K+1,D];
  endif;
  
  EndTime := CpuTime();
  TimeTaken := EndTime - StartTime;
  PrintLn "Computation time: ", DecimalStr(TimeTaken);

  Return Out;
enddefine;