\\ Pari/GP script to compute terms for OEIS sequence A061343,
\\ program by Joerg Arndt, 2013-May-09.
\\ See
\\   R. M. Thrall: A combinatorial problem,
\\   Michigan Mathematical Journal, vol.1, no.1, pp.81-88, (1952).
\\   http://projecteuclid.org/euclid.mmj/1028989731

\\default(echo,1);

dm(v, m)=
{  \\ return product of differences of distinct pairs
    my(p=1);
    for (j=1, m, for (i=j+1, m, p *= ( v[i] - v[j] ); ); );
    return(p);
} /* ----- */

dp(v, m)=
{  \\ return product of sums of distinct pairs
    my(p=1);
    for (j=1, m, for (i=j+1, m, p *= ( v[i] + v[j] ); ); );
    return(p);
} /* ----- */

pfact(v, m) = prod(j=1,m,v[j]!);  \\ product of factorials

\\ Thrall's formula (4) at page 83:
ntab(v, m, nf)= ( nf * dm(v,m) ) / ( pfact(v,m) * dp(v, m) );
\\ A061343(n) is the sum of this function over
\\ all partitions into distinct parts.

pdist(n)=
{
\\ Iterate over all partitions of n into distinct parts
\\ (represented as non-decreasing compositions).
\\ Algorithm by Joerg Arndt.

    my( tct=0, nf=n! );
    my( m, a );  \\ number of parts, vector with parts
    my( s, k );  \\ aux
    \\ init:
    m = floor((sqrt(1+8*n)-1)/2);
    if (m==0, m=1);  \\ handle n==0 as one part ==zero
    a = vector(m);
    s = n;
    for (j=1, m-1, a[j]=j; s-=j);
    a[m] = s;

    while ( 1,
        \\ visit partition:
\\        print(vector(m, j, a[j]));
        tct += ntab(a, m, nf);

        if ( m<=1,  break() );  \\ current partition is last

        s = a[m] + a[m-1];
        k = a[m-1] + 1;
        m -= 1;

        while ( s >= ( k + (k+1) ),
            a[m] = k;
            s -= k;
            k += 1;
            m += 1;
        );
        a[m] = s;

    );
    return(tct);
} /* ----- */

a(n)=pdist(n);

vector(32, n, a(n))
quit;

\\ format for b-file:
for (n=1, 101, print(n, " ", a(n)) );
quit;
\\ end of file