The previous post stated a formula for f(n), the nth square triangular number (i.e. the nth triangular number that is also a square number):
((17 + 12√2)n + (17 – 12√2)n – 2)/32
Now 17 – 12√2 is 0.029… and so the term (17 – 12√2)n approaches zero very quickly as n increases. So the f(n) is very nearly
((17 + 12√2)n – 2)/32
The error in the approximation is less than 0.001 for all n, so the approximation is exact when rounded to the nearest integer. So the nth square triangular number is
⌊((17 + 12√2)n +14)/32⌋
where ⌊x⌋ is the greatest integer less than x.
Here’s how you might implement this in Python:
from math import sqrt, floor
def f(n):
x = 17 + 12*sqrt(2)
return floor((x**n + 14)/32.0)
Unfortunately this formula isn’t that useful in ordinary floating point computation. When n = 11 or larger, the result is needs more bits than are available in the significand of a floating point number. The result is accurate to around 15 digits, but the result is longer than 15 digits.
The result can also be computed with a recurrence relation:
f(n) = 34 f(n − 1) − f(n − 2) + 2
for n > 2. (Or n > 1 if you define f(0) to be 0, a reasonable thing to do).
This only requires integer arithmetic, so it’s only limitation is the size of the integers you can compute with. Since Python has arbitrarily large integers, the following works for very large integers.
def f(n):
if n < 2:
return n
return f(n-1)*34 - f(n-2) + 2
This is a direct implementation, though it’s inefficient because of the redundant function evaluations. It could be made more efficient by, for example, using memoization.
That’s like the extremely inefficient naive way of computing the fibonacci sequence recursively. A better way would be
def f(): a, b = 0, 1 yield a yield b while True: a, b = b, 34*b - a + 2 yield b(based on a function from http://stackoverflow.com/a/499245/161801). If you want the nth term, you can use the nth() recipe from https://docs.python.org/3/library/itertools.html.
If n < 2, wouldn't "return n" be more succint than "return [0, 1][n]"?
Peter: You’re absolutely right. I updated the code per your suggestion.
I initially wrote the code starting at 1 and returned [1, 36][n] for n < 3. I changed the code blindly to handle 0 as an argument without stopping to think how that would let me simplify the code.
My version of Aaron’s approach, in Perl 6:
my @st := 0, 1, -> $a, $b { $b * 34 – $a + 2 } … *;
I went ahead and wrote the mathematical formula that you gave in Perl 6:
sub circumfix:« ⌊ ⌋ » ( Numeric() $_ ) { .floor } sub prefix:« √ » ( Numeric() $_ ) is equiv(&[**]) { .sqrt } # a cheat so that we don't need a * between 12 and √ sub infix:« √ » ( Numeric() $a, Numeric() $b ) is equiv(&[**]) { $a * $b.sqrt } sub square-triangle ( Int $n ) { # another cheat to rid it of ** sub postfix:« ⁿ » ( $a ) { $a ** $n } ⌊((17 + 12√2)ⁿ +14)/32⌋ }Since
sqrtforces it to a floating point number it has the same limitations you stated about the Python version. Also the use of a postfix sub for exponentiation by n (which makes it so that the calculation can look exactly like the formula you provided) will have some performance drawbacks until and unless spesh figures out that it can cheat and inline the exponentiation. ( It may not be able to cheat because the sub is at least theoretically different each time the outer function is called. )I probably would have written the previous Perl 6 comment using placeholder variables like this:
my @st ::= 0, 1, { $^b * 34 - $^a + 2 } ... *;Using a Whatever lambda instead makes it … well lets just say cryptic:
A more idiom for idiom translation from the Python example would be:
( Do not even attempt to use on a pre-GLR Rakudo, and absolutely make sure you have the
lazyprefix! I may or may not have found this out the hard way. )sub f (){ lazy gather { my ($a, $b) = 0, 1; take $a; take $b; loop { ($a, $b) = $b, 34*$b - $a + 2; take $b; } } }or the more Perlish:
sub f (){ lazy gather loop { state ($a,$b) = take(0), take(1); ($a, $b) = $b, 34*$b - $a + 2; take $b; } }You could also assign the lazy gather loop directly to an array variable if you want. “
my @st = lazy gather loop { ... }”“
loop” is the rough equivalent of the “for” loop in C, I just left off the unneeded “( ; ; )” part.( Note to self: don’t use <pre><code></code></pre> on this site )
@Brad: Thanks. Perl 6 is intriguing. I’d like to look into it for fun, even though the bulk of my work is in other languages.
I edited your comment to change most of your <code> tags to <pre> tags. <code> is an inline element and <pre> is a block element. So <code> is what you want inside a paragraph, and I left those alone. But <pre> is what you want around a block of code.
On StackOverflow putting a
<code>inside of a<pre>is how you enable syntax highlighting if you are using HTML instead of Markdown. I think there is another site that I post on that works similarly so you can see why I would do that. ( I generally prefer Markdown, so much so, that I format my plain-text emails that way as well )I think you should at least provide a link to a page that discusses how to format responses.
The most up-to-date info on Perl 6 is on http://docs.perl6.org/ other information is linked on http://perl6.org. The main discussion is at #perl6 on freenode.net ( I think you would have the most in common with DrForr there )
Interesting; I thought this was a case where the intermediate pieces got too big, but that the end product wasn’t that big. So I figured replacing x**n with exp(n*log(x)) would get around the need for bignums. Turns out the 500th number in this sequence is just too big for this trick to help. Thanks for inducing me to write my first python program in at least 10 years!
Well if you care about efficiency, then you want sublinear time. This is O(log n) time to calculate f(n), in Haskell:
g 0 = (1, 0, 1, 0)
g 1 = (17, 12, 17, -12)
g n
| even n = let (a1, b1, a2, b2) = g (n `div` 2)
in (a1*a1 + 2*b1*b1, 2*a1*b1, a2*a2 + 2*b2*b2, 2*a2*b2)
| otherwise = let (a1, b1, a2, b2) = g (n-1)
in (17*a1 + 24*b1, 12*a1 + 17*b1, 17*a2 – 24*b2, -12*a2 + 17*b2)
f n = let (a1,b1,a2,b2) = g n in (a1 + a2 – 2) `div` 32
I wrote a blog post on the same subject a few years ago; you might find it interesting. I took a somewhat different approach, although of course I arrived at the same answer.
http://blog.plover.com/math/square-triangular.html
Mark, Nice post.