111 lines
3.1 KiB
ObjectPascal
111 lines
3.1 KiB
ObjectPascal
program CountAlkanes;
|
||
|
||
{$mode objfpc}{$H+}
|
||
|
||
uses SysUtils; // only for output
|
||
|
||
type TArrayUint64 = array of uint64;
|
||
{
|
||
Function to count alkanes, based on: Shinsaku Fujita,
|
||
"Numbers of Alkanes and Monosubstituted Alkanes.
|
||
A Long-Standing Interdisciplinary Problem over 130 Years",
|
||
Bull. Chem. Soc. Jpn. Vol. 83, No. 1, 1–18 (2010)
|
||
}
|
||
function CountAlkanes() : TArrayUint64;
|
||
const
|
||
MAX_RESULT_INDEX = 49; // as far as this code can get without multi-precision
|
||
MAX_R_INDEX = MAX_RESULT_INDEX div 2;
|
||
var
|
||
R : array [0..MAX_R_INDEX] of uint64;
|
||
nrCentUnb : uint64; // number of centroidal unbalanced alkanes
|
||
temp : uint64;
|
||
m, n, h, i, j, k : integer;
|
||
begin
|
||
SetLength( result, MAX_RESULT_INDEX + 1); // zero-based
|
||
{
|
||
Calculate enough of the coefficients R[], where the generating function
|
||
r(x) = R[0] + R[1]x + R[2]x^2 + R[3]x^3 + ... satifies
|
||
r(x) = 1 + (x/6)[r(x)^3 + 2r(x^3) + 3r(x)r(x^2)] (Fujita, equation 4)
|
||
}
|
||
R[0] := 1;
|
||
n := 0;
|
||
repeat
|
||
if (n mod 3 = 0) then temp := 2*R[n div 3]
|
||
else temp := 0;
|
||
for j := 0 to (n div 2) do begin
|
||
inc( temp, 3 * R[j] * R[n - 2*j]);
|
||
end;
|
||
for j := 0 to n do begin
|
||
for k := 0 to (n - j) do begin
|
||
inc(temp, R[j] * R[k] * R[n - j - k]);
|
||
end;
|
||
end;
|
||
Assert( temp mod 6 = 0); // keep an eye on it
|
||
inc(n);
|
||
R[n] := temp div 6;
|
||
until (n = MAX_R_INDEX);
|
||
{
|
||
Now use the generating function
|
||
(x/24)[r(x)^4 + 3r(x^2)^2 + *r(x)r(x^3) + 6r(x)^2r(x^2) + 6r(x^4)]
|
||
where inserting r(x) up to the term in x^m will give the number of alkanes
|
||
of orders 2m+1 and 2m+2, as the coefficients of x^(2m+1) and x^(2m+2).
|
||
|
||
Note: In Fujita's paper, equation 23, the factor is 1/24 not x/24,
|
||
but x/24 seems to be needed to give correct results.
|
||
}
|
||
result[0] := 1; // conventional
|
||
for n := 1 to MAX_RESULT_INDEX do begin
|
||
m := (n - 1) div 2; // so n = 2*m + 1 or 2*m + 2
|
||
temp := 0;
|
||
|
||
// These loops are written for clarity not efficiency
|
||
for k := 0 to m do begin
|
||
for j := 0 to m do begin
|
||
for i := 0 to m do begin
|
||
h := n - 1 - i - j - k;
|
||
if (h >= 0) and (h <= m) then inc( temp, R[h]*R[i]*R[j]*R[k]);
|
||
end;
|
||
end;
|
||
end;
|
||
|
||
if Odd(n) then begin
|
||
for k := 0 to m do begin
|
||
inc( temp, 3*R[k]*R[m - k]);
|
||
end;
|
||
end;
|
||
|
||
for k := 0 to (n - 1) div 3 do begin
|
||
j := n - 1 - 3*k;
|
||
if (j <= m) then inc( temp, 8*R[j]*R[k]);
|
||
end;
|
||
|
||
for k := 0 to m do begin
|
||
for j := 0 to m do begin
|
||
i := n - 1 - 2*k - j;
|
||
if (i >= 0) and (i <= m) then inc( temp, 6*R[i]*R[j]*R[k]);
|
||
end;
|
||
end;
|
||
|
||
if (n mod 4 = 1) then inc( temp, 6*R[(n - 1) div 4]);
|
||
|
||
Assert( temp mod 24 = 0); // keep an eye on it
|
||
nrCentUnb := temp div 24;
|
||
if Odd(n) then
|
||
result[n] := nrCentUnb
|
||
else begin
|
||
temp := R[n div 2];
|
||
result[n] := nrCentUnb + (temp*(temp + 1) div 2);
|
||
end;
|
||
end;
|
||
end;
|
||
|
||
// Call function and display the results
|
||
var
|
||
nrAlkanes : TArrayUint64;
|
||
k : integer;
|
||
begin
|
||
nrAlkanes := CountAlkanes();
|
||
for k := 0 to Length( nrAlkanes) - 1 do
|
||
WriteLn( SysUtils.Format( '%2d %d', [k, nrAlkanes[k]]));
|
||
end.
|