So, generating $q each time and comparing to $diff, also comparing totals $pp and $p.
Code:
#!/usr/bin/perl -w
$a = 1500;
$b = 1000;
$m = 400;
$diff = $b;
$p = $a;
$pp = $a;
$ddiff = 1;
for ($i = 1; $i <= 40; ++$i) {
$ddiff += $i;
$diff += $m*$ddiff;
$p += $diff;
# $diff = sum(1 to n) (caterer(i))
# i.e. sum(1 to n) i^2 / 2 + sum(1 to n) i/2 + sum(1 to n ) 2
# sum(1 to n) i^2 / 2 = (n * (n + 1) * (2n + 1)) / 6 / 2
# sum(1 to n) i / 2 = (n*(n+1))/2/2.
# sum(1 to n) 2/2 = n.
$q1 = ($i * ($i+1) * (2*$i + 1)) / 6; # sum (n^2)
$q2 = ($i * ($i+1)) / 2; # sum (n)
# must add before dividing by 2, to ensure even
$q = (($q1 + $q2) / 2 ) + $i; # /2 + sum 1
$pp += $b + $q*$m;
printf("%4d%12d%12d%12d%12d%12d\n", $i, $p, $diff, $ddiff, $b+$q*$m, $pp);
}
So now we need to take the sum of each component of the total sum to arrive at a formula for each iteration of the loop. I think it involves sum of cubes.
All-in-all, the formula is probably more clearly expressed as an algorithm.