Clever Programmers: I need to improve the efficiency of this program by 8 orders of magnitude :(

[Shalmanese]

Okay, we have this programming lab assignment where we have to write a recursive function to find out the first 20 numbers that are both square and triangular and it has to do it in less than 60 seconds o a P3 550.

So I go on the web and do a google and come up with this recursive formula to find out all square triangular numbers:

a(n) = 34 * a(n - 1) - a (n - 2) + 2

so I just whip up a quick code and everythings working... whatever.

However, today, in our lecture, our lecturer basically said that we cant use any formulas that would not be apparent to somebody without any deep understanding in number theory and a friend of mine who had a Lab on Wednesday who used that formula did not get the mark. His lab teacher said that what he was doing was "black magic".

So, I have come up with a non-black magic solution where it says:

starting at n = 1, m = 1

if the nth square is larger than the mth triangular number, add one to m and repeat recursivly

if the nth square is less than the mth triangular number, add one to n and repeat recursivly

if m = n, print it out, add one to both m and n and keep going until you have printed out 20 numbers.

The problem is, the 20th square triangular number is 361786555939836 ^2 and my algorithm is doing roughly 400,000 numbers per minute so it would take roughly 700 years for me to get up to the 20th number this way.

Can anybody think of a faster way to run this task so that it could get under 60 seconds?

 

[rgwalt]

Could you clue me in on what a triangular number is?

Edit: NM, found it for myself. Let me think about this.

Ryan

 

[rgwalt]

Could you just find the various triangular numbers and test to see if they were perfect squares?

Edit: I tried this. Fairly clever, but I overflowed my sqrt function after finding the 7th number. There has to be a better method to do this. Sorry I couldn't help. Good luck!

Ryan

 

[Mark R]

To continue with rgwalt's algorithm - this seems like the way to do it.

However, you don't actually need to evaluate the triangular number to assess it for squareness.

Try using a loop which finds prime factors of non-prime n and n+1 where 2T(n) = n.n+1

Note though that this is not the sort of problem where I can see a recursive loop being appropriate.

 

[m0ti]

Here's something helpful:

1/2*((n+1)*(n+2) - (n)*(n+1)) = 1/2*(n+1)*(n+2-n) = 1/2*(n+1)*2 = n+1.

This is actually sort of obvious, if one thinks about that we're talking about the difference between adding 1...n, and 1...n+1. obviously, the latter sum is larger by n+1.

Let's assume that we're at the m'th square number and it's larger than our current n'th triangular number. don't bother increasing n by 1 and then checking...you know that it'll be n+1 bigger than the previous number.

So what we do is: let's say that the difference between the m'th square and the n'th triangle is x.

let's assume that we have to add y on to n so that the triangle works out. that means that the (n+1) + (n+2) + ... + (n+y) = x. or that y*(y+(n+1)) = 2x. That can be quickly solved. If there is no solution, then set n to the newest lowest possible guess (based on the solutions to the equation).

Now then, we have to find a solution to incrementing m, in the other direction. if the n'th triangle is greater than the m'th square, how can we increment m?

again assuming x is the difference and we want to increment m by y: (m+y)^2 -m^2 = 2my + y^2 = x. again, easily solvable. if there isn't a whole number solution then increment m appropriately.

 

[Sahakiel]

I can think of a few ways.

1. Use assembly and pit your programming forte against the compiler. Might be advisable in any case, since you could then specialize the code to get the most out of the processor.
2. Use a different algorithm. Recursive functions require lots of memory, but I'm not sure if this results in your program spilling over into main memory or still being stored in the caches. Basically, make sure your program stays in cache, L1 if possible.
3. If you insist on scanning numbers one at a time, use two algorithms. One for N<2^32 and one for N>2^32. I'd suggest finding a better alternative to scanning.
4. I'm not big on number theory, but you may want to look into other derivation methods.
5. Bribe your prof. If your prof is male and you're female, you're in luck. If your prof is female and you're male, grin and bear it. If you're both female and sexy, invite me. If you're both male, stay as far away from me as possible. Otherwise, you're on your own.

 

[blahblah99]

Triangular numbers are given by the formula f(x) = x*(x+1)/2 - test this out and convince yourself that this is correct.

If you want a number that's triangular AND square, then f(x)^2 = x*(x+1)/2.... Using the quadratic formula, solve for x and you get [-1 + sqrt( 1 + 8*f(x)^2)]/2. Your only taking the positive solution since you only want postive triangular, square numbers. In order for x to be a positive integer, the quantity inside sqrt MUST be a square value.

And your recursive function would just increment a counter and test if that value is a square number. If it is, plug it into the equation for x and spit out the value. I'm not sure if your professor would consider this as requiring "deep understanding in number theory" as it is just basic high school math. But I dont see why he would want you to run a recursive function when a simple for() loop would suffice.

If the assignment is to test your knowledge on reducing algorithm speed, then for sure he is not going to let you use what I described above - I'll have to think about it a little more

 

[Turkey]

Remember, when you're calling a recursive function you don't have to analytically calculate a square or triangular number, just the next square or triangular number.

So the square numbers are given by the recursive formula

sq(1) = 1
sq(n+1) = sq(n) + 2n + 1

Triangular numbers are given by the recursive formula

tri(1) = 1
tri(m + 1) = tri(m) + m + 1

I think if you used your original program with these formulas instead of the analytical formulas you would have no problems. You might even be able to use reference parameters instead of value parameters if you coded it right.

EDIT: just tried it, I run out of stack space after the 7th number using g++ and HP-UX 10. I'd guess you'd have to use some better algorithm for incrementing N and M than just adding one, like the one m0ti suggested.

 

[m0ti]

I think that he's got a bit of a problem with most of the solutions.

If the 20th square triangle is that big (361,786,555,939,836 ^2), there is no way he can find it in 60 seconds on that machine if he goes through the squares one at a time. Let's assume that he only had to summations, and he only had to do a single sum per step. Let's even say that he can do 1 Trillion summations per second (I wish I had a computer that could do that many). It would still take him 6 minutes!

Therefore, there has to be some shortcut math involved so to not have to test all of them.

 

[blahblah99]

Like I said, solve the equation I gave you in the earlier post and you should be able to find 20 triangular square numbers fast.

 

[m0ti]

my statement is still correct. you're still checking every single triangular number to check if it's value squared times 8, plus 1 is a square. Therefore, you're still taking the numbers one at a time (f(1), f(2), f(3)...), and the calculation will take much too long. There has to be some process whereby you can eliminate large amounts of numbers based on previous answers, in order to do the calculation fast enough.

 

[Agent004]

If it's recursive, usually(I do mean usually, don't kill me on this ) all you need is the first two terms, as all the terms after will be given in terms of those.

I think both Turkey and m0ti are on the right track.

 

[DuffMan]

You are probably already doing this, but make sure you use full optimization with your compiler.

 

[blahblah99]

Ok, in case you havn't solved this problem yet... I worked it out a little more.

And you guys are correct, you can't just compute the square of that formula. I tried compiling a simple C code using gcc and it gave me an error at output. SOOOOOOOOOOOO
with a little research and linear algebra, it turns out the solution to this problem is a classic example of PELL'S EQUATION.

You want to solve for 1+8*y^2 = f(y)^2. Rewriting that gives you f(y)^2 - 8*y^2 = 1. Solve that equation by finding the continued fraction of 8 (which is the coefficient in y^2)

Let f(y) = z
Doing that leads to this difference equation:

1. z[k] = 6*z[k-1] - z[k-2]. Rewriting this gives z[k+2] - 6*z[k+1] + z[k] = 0

The solution to 1. has the form A*c1^k + B*c2^k

Characteristic equation is s^2 - 6s + 1

Using quadratic formula, you get s = 3+ 2sqrt(2), 3-2sqrt(2) and those are the c1, and c2. Now going back to the equation f(y)^2 = 1+8*y^2, you can find the initial values. Let y = 0 then f(0) = 1. Let y = 1 then f(1) = 3. Plug that back into A*c1^k + B*c2^k and solve for A and B and you get A=B=1/2.

z[k] = 1/2 * (3+2sqrt(2))^k + 1/2 * (3-2sqrt(2))^k, which gives n[k] = (z[k] - 1) / 2, where n[k] is the the value of (-1+sqrt(1+8*y^2)) / 2.

For example, if

k=0 z[0] = 1, n[0] = 0
k=1 z[1] = 3, n[1] = 1
k=2 z[2] = 17, n[2] = 8
k=3 z[3] = 99, n[3] = 49

If you continue with k, you'll notice that n[k] alternates between odd and even numbers, and this is shown in n[k] = 6*n[k-1] - n[k-2] + 2. Also, from the original equation,
f(n)^2=n(n+1)/2 ==> 2*f(n)^2 = n*(n+1) - which shows that x and x+1 are co-prime numbers so either n or n+1 must be an odd square and the other one HAS to be twice the square. As a result, k MUST BE ODD in order for n(n+1)/2 to be a square number.

Well, after all that confusing math ,

you get 2*n[k] = (3+2sqrt(2)^k + (3-2sqrt(2)^k, where k is an odd integer, and n[k] is the k-th value of n that you can plug into f(n)^2 = n*(n+1)/2 to get your numbers.

EDITED: DOH, A QUICK SEARCH ON GOOGLE YIELDED this

 

[skriefal]

<< Could you just find the various triangular numbers and test to see if they were perfect squares? >>



Yep. In fact, I just tried exactly that algorithm to solve the OP's problem. Using int's, it'll only calculcate the first 35 square triangles, and then the values become too large for an int to hold.

It's quite fast -- my function that calculates the first 35 square triangles can run 10000 times in a little over 1 second on an AMD Athlon XP 1600+.

Edit: Doh! I was thinking that by "square triangular" you meant a number for which both the cubed root and square root are whole integers. Oh well.

 

[Shalmanese]

<<
So I go on the web and do a google and come up with this recursive formula to find out all square triangular numbers:

a(n) = 34 * a(n - 1) - a (n - 2) + 2

so I just whip up a quick code and everythings working... whatever.
>>



thats not it, I need something that can be UNDERSTOOD by a person with understanding of basic number theory.

 

[DocBurner]

I think I have an useable solution for you, let me code it up and see if it works.

 

[bizmark]

I think that well before 20 square triangular numbers, you're going to hit a wall with the bit resolution of the processor. I don't think it's possible. My story:

I wrote a program trying to take advantage of the fact that the triangular numbers are very simply incremented, i.e. the first 5 triangular numbers are 1, 3, 6, 10, and 15, but 3=1+2, 6=3+3, 10=6+4, 15=10+5, etc. (see the number is just one plus the last increment...) and the same idea for squares 1, 4, 9, 16, 25, so the differences here are 3, 5, 7, 9 respectively, which is just 3+2(+2)(+2)(...). This was summarized by Turkey. I guess I just did the same thing that he did. Note that this was NOT a recursive algorithm.

Anyway, on a 64-bit Solaris system at my University (not sure the exact processor speed, but it sure ain't state-of-the-art) (using the "long long" integer implementation), I got the first 10 triangular squares (up to 63955431761796) within like 5 seconds. Then it seemed to hang; the next one (2172602007770041) came after about 30 more seconds. 73804512832419600 came at like 2 minutes into the program. And then after like 10 more minutes of waiting, I got the 13th square triangular number, 2507180834294496361. Then I killed the program. Total running time was ~13 minutes. I tried switching the type to double and I got the same answers, but once it got past 2.1726e+15, I started getting false positives left and right.

(Trying a "long double" implementation, I got impatient after a few minutes of letting it think --it took a lot more time to get any of the answers, and I quit before it even got as far as the normal double did.)

Trying all of this with "long" integers sucked. It hit the bit-wall very quickly, and where you'd expect it to: ~2^32.

So basically I don't know how this problem is really possible, given that you have a 32-bit processor on which to do your work, and we need to represent numbers that are in the range of 10^30=2^100 (i.e. 100-bit precision) if m0ti's statement about the 20th triangular number was correct.

Is there an implementation of integers or real numbers which overcomes the limits of the processor? I mean, you'd think there'd be something available, that could essentially split up 128 into 64 and 64, or 64 into 32 and 32....

Other ideas, just for general thought (I wrote all of these before I tried my program and realized the bits thing):

I think that it might be worth trying to eliminate some "tries" from your algorithm. Sure, this may slow it down a bit, but I think that the gain is worth the loss.

Like if you're going through computing triangular numbers, you can pre-test to see if they're candidates to be squares. No square will ever have a 2, 3, 7, or 8 in the ones space. That eliminates 2/5 of your tries already, assuming that these numbers come up regularly in triangular numbers (I don't know for sure whether they do, but 3 and 8 show up already in the first 9 triangular numbers). Of course, this depends on exactly what algorithm you're using to find these numbers; with my algorithm above, this type of checking doesn't make much sense.

So you're doing this little comparison, and yes, parsing the number may be a pain and take some time, but it would seem that computing the square root (or, alternatively, trying numbers to see if, when squared, they are this number) would be more computationally intensive than parsing the number and seeing whether it's a 2, 3, 7, or 8. Would the gains outweigh the losses? I dunno, and it may depend on the particular language/compiler/processor setup you're working with.

I know we can't see an order-of-magnitude gain from this type of optimization, but still it's something worth considering.

Another way is to do basically what blahblah99 said. So we really just need to look for m's, such that 1+8m^2 is also a perfect square. If it's not a perfect square, then there's no way in hell that n is an integer. If it is, though, we've got a pretty good candidate for a number that is also a triangular number.

Last addition to this post:
So that we can look for patterns to aid in designing a recursive algorithm, here are the first 13 triangular square numbers:

1 (n=1, m=1)
36 (n=6, m=8)
1225 (n=35, m=49)
41616 (n=204, m=288)
1413721 (n=1189, m=1681)
48024900 (n=6930, m=9800)
1631432881 (n=40,391, m=57,121)
55420693056 (n=235,416, m=332,928)
1882672131025 (my calculator won't get bigger than 12 digits, and I'm very tired right now!)
63955431761796
2172602007770041
73804512832419600
2507180834294496361

So looking at this, I notice the following:
the pattern for the ones digit is: 1, 6, 5, 6, 1, 0, 1, 6, 5, 1, 0, ... (so we can eliminate also 9 and 4 from consideration)
41616/36=1156
48024900/1225=39204
55420693056/41616=1331716.
However, 1631432881 is NOT evenly divisible by 1413721. It's very close but no cigar.
So, I thought I saw a pattern, but it's not there. And my calculator won't hold more than 12 digits. Hmmm well anybody have any bright ideas?

 

[blahblah99]

Given the requirement of finding the first 20 square triangular numbers under a minute on a P3, it is impossible to compute by brute force. So in order to solve this problem, you will need to do some math (see my post above).

Y = ( n[k]*n[k] + 1 ) / 2,

where n[k] = ( 3+2sqrt(2))^k / 2 + ( 3 -2sqrt(2))^k / 2, k = odd positive integer.

 

[zetter]

Here's a nice formula to work out the list of square and triangular numbers I found whilst travelling the 'net:

If L is a square-triangular number, then the next one is 1 + 17*L + 6*sqrt(L + 8*L^2) - Lekraj Beedassy (beedassylekraj@hotmail.com), Jun 27 2001

 

[Sukhoi]

I bet that function would work in less than 60 seconds.

 

[RobK]

Really, you're going to have to listen to these guys and use one of those two formulas. Ask the professor - never rely on the lab TA for a final answer!

 

[Shalmanese]

<< I think that well before 20 square triangular numbers, you're going to hit a wall with the bit resolution of the processor. >>



with the Language I am using (Haskell), it has a nice type signature called Integer which can store any integer within 256MB of address space. Resolution is not a problem.

 

[lowwall]

Not sure how much this will help, but looking at wbwither's list, I notice the following pattern:

for any triangular square > 1, the square root of the next triangular square is 6*the square root of the current triangular square - the square root of the previous triangular square

i.e.

6*6 - 1 = 35
6*35 - 6 = 204
6*204 - 35 = 1189
6*1189 - 204 = 6930
6*6930 - 1189 = 40391
...

Obviously the algorithm is trivial, if someone can just come up with a non magic reason why this is true, we are all set.

 

[bizmark]

<< for any triangular square > 1, the square root of the next triangular square is 6*the square root of the current triangular square - the square root of the previous triangular square >>



Good work. Reading back through, I see that that's what blahblah99 said earlier (although I didn't understand it at the time):



<< z[k] = 6*z[k-1] - z[k-2]. >>



So blahblah, you wanna explain the math here so we can judge whether it's voodoo or not?

I was up with you until the continued fraction of 8 thing.... good thing I'm a third-year f*cking math major at a top-10 university We don't learn anything at all practical. Or maybe I just learned it once and forgot