Simple path finding algorithm in assembler

[sn3j]

So I've been working on a path finding algorithm (in Z80 assembler) in the last days, here's the results.
The algorithm takes a 32x40 tiled terrain map as input and computes the shortest path between two positions.
Performance is about 6 frames per 1K visited tiles.
The average workload - to find the path for two positions halfway across a map - is 500 visited tiles, so 3 frames.
You could use this algorithm to build games like Ultima.

Ram requirements:

• 1280 bytes of terrain data
• 1280+512 bytes workspace
• 450 bytes code

If anyone's interested I can pm the code. Here's how it works:

• there are 4 queues: L, L+1, L+2, L+3, initially all empty
• start position goes into the L-queue and so becomes the first item to work on
• foreach item in the L-queue do this:
  - all diagonal unvisited neighbor tiles (NE,SE,SW,NW) are being added to the L+3-queue
  - all straight ones (N,E,S,W) to the L+2-queue
  - any neighbor added to a queue is also marked as 'visited'
• when the L-queue is done, all queues shift: L+1-queue is now L-queue, L+2-queue is now L+1 queue, L+3 gets L+2, and L+3 is set to empty
• repeat working on the new L-queue, unless the goal position is found, or the new queues contain no more items.

[Einar Saukas]

There's something similar called PATHFINDER:

The PATHFINDER library is a small efficient Assembly implementation of the "simplified Dijkstra" algorithm for finding shortest paths in small mazes, that you can freely use in your own programs. It uses 512 bytes as temporary queue and only 89 bytes of code

The incremental PATHFINDER works exactly like the original full version, except it will only solve a limited number of maze steps each time it's called. This way, an action game that calls the incremental version only once every 3 frames would be able to recalculate all paths in a full screen maze about once per second. This should work just fine for ghosts chasing a constantly moving player in a Pacman game, for instance.

[ParadigmShifter]

And if you want fast instead you can't go much wrong with A*. (Can use Chebyshev distance d = max(|x1-x2|, |y1-y2|) as heuristic or the taxicab metric d = |x1-x2| + |y1-y2|). Where (x1, y1) is current position, (x2, y2) is goal position.

|a| is absolute value of a.

If you do A* you'll want to implement a binary heap though for an efficient priority queue.

I had a heap implementation for a heap sort but it had a bug if there were more than 64 entries I think (which wasn't going to be a problem since I'd likely have less than that amount of things to process... can post the code if you like, bug is probably 8 bit maths instead of 16, can't remember).

Implementing a heap (data structure, nothing to do with the heap in C) is not hard. A heap is a binary tree which does not need pointers to each node since they can be calculated based on where you are in an array. A min heap or a max heap is a way of putting values into the heap so that they aren't guaranteed to be fully sorted but an element with the minimum value (for a min heap, maximum value for a max heap) is always at the root, so ideal for a priority queue.

Pacman doesn't really use a search it just tries to move towards a goal position each time it reaches an intersection (and can't turn around although scatter mode makes all ghosts reverse direction immediately IIRC). So if the goal position is further away in the x direction than the y direction it will move horizontally towards it if it can... otherwise move vertically towards goal. Can't remember what it does if the x and y differences are the same... keep going straight ahead I think? Goal positions are different for each ghost which is what makes the AI more fun than always homing in (which would just result in all the ghosts being on top of each other pretty soon)... see https://gameinternals.com/understanding ... t-behavior for a full explanation. Targets don't even need to be inside the maze path either.

[sn3j]

There's something similar called PATHFINDER:

Ah, I've looked up the DB for Dijkstra and A-Star before trying myself, but the search didn't turn up results.
I like the idea of storing the path lengths directy in the maze. I am using 8-way movement (instead of the 4-way in your algorithm) but this could be worked out by adding either 2 or 3 to the path (instead of just 1), keeping the 1:1,5 ratio.

BTW this

Code: Select all

        ld      de, WIDTH
        scf
        sbc     hl, de

could be replaced with

Code: Select all

        ld      de, -WIDTH
        add     hl, de

if it's not too obfuscated.

[sn3j]

I had a heap implementation for a heap sort but it had a bug if there were more than 64 entries I think (which wasn't going to be a problem since I'd likely have less than that amount of things to process... can post the code if you like, bug is probably 8 bit maths instead of 16, can't remember).

I'd be interested to see how you implemented it. My version works with max 128 items and uses 512 bytes storage.

And if you want fast instead you can't go much wrong with A*. (Can use Chebyshev distance d = max(|x1-x2|, |y1-y2|) as heuristic or the taxicab metric d = |x1-x2| + |y1-y2|). Where (x1, y1) is current position, (x2, y2) is goal position.

For my A* implementation I used a different heuristic, I named it the tile-distance. There's probably a proper name for it but I didn't want to spend too much time studying this. It's a compromise between Chebyshev and Taxicab, because I'm doing 8-way movements only I cannot use Taxicab, and Chebyshev has flaws, too.
That is, both Taxicab and Chebyshev do not incur costs for zig-zagged paths under some conditions, resulting in unnatural paths.

[Evil Genius]

Idle curiosity: can this be adapted for multiple terrain types with differing movement costs, even at a limited range of say, 6 or 8 hexes?

[Einar Saukas]

And if you want fast instead you can't go much wrong with A*.

A* "optimality" makes it faster for large mazes. However Spectrum games don't have large enough maps to compensate it.

For smaller mazes, it's probably faster to just solve the entire maze instead of spending time with heuristics. As mentioned in PATHFINDER documentation:

Code: Select all

For such small maps, the overhead from using a fairly complex heuristics-based
pathfinding algorithm (in terms of data structure management time and memory) is 
not worth it. Instead, these maps can be easily solved by a "simplified Dijkstra"
algorithm.

Another advantage of PATHFINDER is supporting multiple sources and multiple targets (unlike A*). Thus you only need a single execution for all ghosts to find the player.

Pacman doesn't really use a search it just tries to move towards a goal position each time it reaches an intersection

That's because the original Pacman is not really a "maze" you need to solve. It's basically an open area filled with small obstacles to walk around. Therefore it doesn't really need a maze solver, it's enough to walk in the general direction towards the goal. That wouldn't work in other Pacman variants with harder mazes.

I mentioned Pacman in the documentation just because it made easier to understand the idea of several enemies chasing a common target in a maze. But a pathfinder is certainly more useful in more complex mazes.

Goal positions are different for each ghost which is what makes the AI more fun than always homing in (which would just result in all the ghosts being on top of each other pretty soon)

They all move in the general direction of the player, with small differences. Probably the main reason ghosts don't get on top of each other is, from time to time, they stop chasing the player and each ghost move towards a different corner of screen ("scattering"). This way, next time they start chasing the player again, each one will be coming from a different direction.

[sn3j]

Idle curiosity: can this be adapted for multiple terrain types with differing movement costs, even at a limited range of say, 6 or 8 hexes?

If the movement costs do not differ substantially, maybe. The algorithm could be expanded to work with, say, 8 queues instead of 4. Now, when a neighbor tile is to be added, it must still fit into one of the L+1 .. L+7 queues. So with costs in the range of 1..7 you would need 8 queues, in the range of 1..15 you already need 16 queues. And this is going to hurt, implementation wise. Somewhere, you'll end up spending more time in managing the queues than you had done in maintaining a proper heap / priority queue. It's a sort of struggle between tradeoffs and gains.

[Einar Saukas]

Idle curiosity: can this be adapted for multiple terrain types with differing movement costs, even at a limited range of say, 6 or 8 hexes?

In case of PATHFINDER, support for weighted terrain would require using one priority queue (implemented as binary heap) instead of a simple queue.

[Einar Saukas]

Ah, I've looked up the DB for Dijkstra and A-Star before trying myself, but the search didn't turn up results.
I like the idea of storing the path lengths directy in the maze. I am using 8-way movement (instead of the 4-way in your algorithm) but this could be worked out by adding either 2 or 3 to the path (instead of just 1), keeping the 1:1,5 ratio.

Exactly! This way, a weighted variant of PATHFINDER could support diagonal movement too.

BTW this

Code: Select all

        ld      de, WIDTH
        scf
        sbc     hl, de

could be replaced with

Code: Select all

        ld      de, -WIDTH
        add     hl, de

if it's not too obfuscated.

The PATHFINDER documentation explains how to change the WIDTH by POKEing certain addresses. I didn't want to make this process any harder by having to explain how to POKE a negative 16-bits number from BASIC, just to save a single byte.

[ParadigmShifter]

I'd be interested to see how you implemented it. My version works with max 128 items and uses 512 bytes storage.



For my A* implementation I used a different heuristic, I named it the tile-distance. There's probably a proper name for it but I didn't want to spend too much time studying this. It's a compromise between Chebyshev and Taxicab, because I'm doing 8-way movements only I cannot use Taxicab, and Chebyshev has flaws, too.
That is, both Taxicab and Chebyshev do not incur costs for zig-zagged paths under some conditions, resulting in unnatural paths.

I only have a heapsort implemented (and there is a bug when it has > 64 items so it may not be that much use. I didn't have more than 64 objects to sort though). That's because it was simple to write compared to other O(n log n) sorts, and bubblesort (easiest to write) is O(n^2).

https://en.wikipedia.org/wiki/Heapsort

If you want to try a heap based priority queue you probably want to have a look at this

https://www.geeksforgeeks.org/priority- ... nary-heap/

Heapsort basically puts all items in a priority queue and then removes them all, they come out in order, item to remove is always the one at index 0. Adding and removal are both O(n log n) O(log n). (Peek at the top of the queue is just the item at index 0 though so is O(1)).

[Einar Saukas]

Adding and removal are both O(n log n).

Actually O(log n).

[ParadigmShifter]

I'd be interested to see how you implemented it. My version works with max 128 items and uses 512 bytes storage.



For my A* implementation I used a different heuristic, I named it the tile-distance. There's probably a proper name for it but I didn't want to spend too much time studying this. It's a compromise between Chebyshev and Taxicab, because I'm doing 8-way movements only I cannot use Taxicab, and Chebyshev has flaws, too.
That is, both Taxicab and Chebyshev do not incur costs for zig-zagged paths under some conditions, resulting in unnatural paths.

Do you mean diagonal movement costs more than orthogonal movement then? (Most people use moving diagonally 1.5 times more expensive for speed and multiply distances by 2)

So distances using that would be (using 2 for orthogonal movement and 3 for diagonal)

Code: Select all

6 7 8 9
4 5 6 8
2 3 5 7
0 2 4 6

Where you are at position 0 and nodes have the cost to get there otherwise? (Other quadrants are just mirrored obvs so I didn't show them)

Distances for Chebyshev (using cost 1 for all movement)

Code: Select all

3 3 3 3
2 2 2 3
1 1 2 3
0 1 2 3

Distances for Taxicab/Manhattan distance

Code: Select all

3 4 5 6
2 3 4 5
1 2 3 4
0 1 2 3

Chebyshev moves on diagonals if it can. Taxicab encourages orthogonal moves.

You'll notice that the one where diagonal moves cost 3 and orthogonal moves cost 2 is just Chebyshev distance + Taxicab distance i.e. it's

d = max(|x1-x2|, |y1-y2|) + |x1-x2| + |y1-y2|


EDIT: Anyway for A* the distance is just used as a heuristic it does not need to be 100% accurate. Once you visit a node you update the value with the actual cost. It just uses the heuristic to determine which way to start looking.

Djikstra's algorithm is equivalent to A* with a heuristic of 0 for all nodes. EDIT2: But obviously it's more efficient to implement Djikstra without actually doing an A* with a heuristic of 0, they are just functionally equivalent.

[ParadigmShifter]

Actually O(log n).

Yup I will correct my post. It's the sort (simply n insertions and then n removals) which is O(n log n).

[sn3j]

Do you mean diagonal movement costs more than orthogonal movement then? (Most people use moving diagonally 1.5 times more expensive for speed and multiply distances by 2)

So distances using that would be (using 2 for orthogonal movement and 3 for diagonal)

Code: Select all

6 7 8 9
4 5 6 8
2 3 5 7
0 2 4 6

Where you are at position 0 and nodes have the cost to get there otherwise? (Other quadrants are just mirrored obvs so I didn't show them)

Distances for Chebyshev (using cost 1 for all movement)

Code: Select all

3 3 3 3
2 2 2 3
1 1 2 3
0 1 2 3

Distances for Taxicab/Manhattan distance

Code: Select all

3 4 5 6
2 3 4 5
1 2 3 4
0 1 2 3

Chebyshev moves on diagonals if it can. Taxicab encourages orthogonal moves.

You'll notice that the one where diagonal moves cost 3 and orthogonal moves cost 2 is just Chebyshev distance + Taxicab distance i.e. it's

d = max(|x1-x2|, |y1-y2|) + |x1-x2| + |y1-y2|

Yes, this is what I came up with, orthogonal cost of 2 and diagonal 3.
The problem I had with Taxicab was that 'Manhattan' paths where equally prioritized in contrast to diagonally straight paths, and vice-versa with Chebyshev:

(Taxicab metric, all paths have the same cost, but Red would be preferred.)

(Chebyshev metric, same thing)

So I did a combination of the two, it is 2*max(|dx|, |dy|) + min(|dx|, |dy|) which I think equals the equation you gave above.

[Einar Saukas]

If anyone's interested I can pm the code.

Or you could release your demo, source code and documentation as a new entry in the archive!

[ParadigmShifter]

Yes, this is what I came up with, orthogonal cost of 2 and diagonal 3.
The problem I had with Taxicab was that 'Manhattan' paths where equally prioritized in contrast to diagonally straight paths, and vice-versa with Chebyshev:

(Taxicab metric, all paths have the same cost, but Red would be preferred.)

(Chebyshev metric, same thing)

So I did a combination of the two, it is 2*max(|dx|, |dy|) + min(|dx|, |dy|) which I think equals the equation you gave above.

Yes that is the same equation. I'll use dx and dy instead of |dx| and |dy| to save annoying myself typing it


Proof:

Case 1: dx <= dy

max(dx, dy) = dy, min(dx, dy) = dx

so
max(dx, dy) + dx + dy = dy + dx + dy = (2*dy) + dx = 2*max(dx, dy) + min(dx, dy)

Case 2: dx > dy

max(dx, dy) = dx, min(dx, dy) = dy

so
max(dx, dy) + dx + dy = dx + dx + dy = (2*dx) + dy = 2*max(dx, dy) + min(dx, dy)

Depends on your implementation of max really which one is better. I'm assuming you don't also evaluate the min. If max returns the max in a register and you have to evaluate the min as well though the first one is better.

Cost value is at most (board maximum extent-1)*3 with that metric so you should be able to use 8 bits for the cost with the board size you are using.

I'm not sure if that metric has a name... Taxicab and Chebyshev are part of an infinite family of metrics

https://en.wikipedia.org/wiki/Lp_space

with p>=1

Which is why Taxicab metric is probably easier to refer to as the L-1 norm and Chebyshev the L-infinity norm. The standard Euclidean metric is the L-2 norm.

EDIT: Whoops 0 != 1 except for small values of 1

[ParadigmShifter]

Anyway it sounds like you have basically implemented A* when the heuristic always gives the exact cost and not an estimate of the cost (which simplifies A* anyway that assumption).

Sounds like Einar's code starts at the goal and flood fills the cost value to each node which as he said is good if multiple agents have the same goal. Also small code and memory overhead. I don't know how he is resolving ties though. A* and Djikstra just return the first shortest path they find.

Since you are basically doing an A* with perfect heuristic anyway you probably want to look into implementing a minheap (which is very easy, only uses shifts to move to descendent left/right or parent nodes) although you have to do quite a lot of swaps, O(log n) in the worst case. Since your algorithm sounds like you always get the exact cost you never have to put a node back into the open list IIRC so you don't need to do the more complicated case in A* where you find your actual cost is lower than the estimated cost (you need to check if the node is in the open list then I think, which is not efficient with a binary heap, but if you have a spare bit in the tile map you can use a 1 bit flag for "is in the open list" I guess but that means you have to clear that bit when you do another search).

I was looking for the most efficient 16 bit swap routine when writing my heapsort, this is what I came up with, probably a much better way. (And I wasn't that great at Z80 when I wrote it, was several years ago).

It swaps (DE) with (HL) and preserves DE, HL. Trashes A and A'
If what DE and HL are pointing at are aligned to a 2 byte boundary can use INC L instead of INC HL etc. for a speedup.
Can be made faster by using B or C or maybe both may be required, ages since I wrote it, if you don't mind trashing those. But you can only do LD (HL), B or C and LD (HL), B or C, you can't use B or C for (DE), have to use A.

Code: Select all

	; DE address of word1
	; HL address of word2
	; Trashes A, A'
	; preserves DE, HL
	; no alignment assumed for DE or HL
	MACRO SWAP_WORD
	ld a, (hl) ; 7
	ex af, af' ; 4
	ld a, (de) ; 7
	ld (hl), a ; 7
	ex af, af' ; 4
	ld (de), a ; 7
	inc hl ; 6
	inc de ; 6
	ld a, (hl) ; 7
	ex af, af' ; 4
	ld a, (de) ; 7
	ld (hl), a ; 7
	ex af, af' ; 4
	ld (de), a ; 7
	dec hl ; 6
	dec de ; 6 = 96
	ENDM

I mainly used that to speed up collision detection between player and static objects (can do a binary search on position then) and I did the sort once at level startup after generating the random obstacles. I did plan on using it to sort the sprites on their Y coordinate later on though. Was one of my earliest attempts at drawing directly to the screen racing the raster, got a lot better at doing that recently though.

[sn3j]

I'm not sure if that metric has a name... Taxicab and Chebyshev are part of an infinite family of metrics

https://en.wikipedia.org/wiki/Lp_space

with p>=1

Which is why Taxicab metric is probably easier to refer to as the L-1 norm and Chebyshev the L-infinity norm. The standard Euclidean metric is the L-2 norm.

So I was using Chebyshev+Taxicab all the time without knowing it. Interesting.
Though it's still not precise enough. The 1.5 is just a bit too far from the sqrt(2).

The gold standard would be Euclidean distance but you'd ruin the performance in any attempt to get a square root.
I was thinking along of using a 5:7 ratio, that is 5*Taxicab+2*Chebyshev, I think.
I'll give the Lp-space article a good read, though I'm afraid it's a bit over my head, but thank you anyway.

[sn3j]

Here's an example where the Chebyshev+Taxicab metric is incorrect.

(Chebyshev+Taxicab: Red gets 13.5, losing against blue; Euclidean cost: 12.7)

The cost for the blue paths is 13, but the red one at (9,9) gets 13.5 which means it loses preference over the blue ones.
With Euclidean distance it would cost 12.7 (the circle segment shows the radius).
This can result in paths slightly deviating from the best path on large maps - only if the best path is roughly diagonal of course.
With small maps the error is not easily noticeable I think, but it's going to have impact on larger scales.

[ParadigmShifter]

That's only going to make a difference when there is a choice of ways to get to the goal though and they are forced to visit either the red or blue node and it picks the blue one incorrectly.

You could always use a lookup table for the distance if you are really worried about it but you'll need to use more than 8 bits for the cost then. You only need 1 quadrant so similar to the diagrams I drew earlier like this, probably best to flip it though like so, so the origin is in the top left corner at (0, 0)

Code: Select all

0 2 4 6
2 3 5 7
4 5 6 8
6 7 8 9

You'll need a lookup table size of ((max map size/2)+1)^2 * precision bits for that (but probs want to align each row to a power of 2 to easier convert row, column to index in table?).

EDIT: Actually you only need a bit more than half of the table (one octant plus a bit more for one diagonal and one full horizontal row) since the table is symmetric i.e d(x,y) = d(y,x) so it is the same as it's transpose (flip along leading diagonal)... but that would make lookup more complicated.

i.e. you could store this

Code: Select all

0 2 4 6
  3 5 7
    6 8
      9

which can pack as

0 2 4 6 3 5 7 6 8 9

I could probably work out what the accessor function for that is if I could be bothered There's n(n+1)/2 entries anyway (and you won't need to store 0 either since d(x,x) = 0, so n(n+1)/2-1 entries to pack a lookup range of -n to n in both axes).

Alternatively (and you'll need a multiply for this of course) you could work out dx^2 + dy^2 (or you can use a table of squares for that too so you don't need a multiply) and then have a square root lookup table so

0 -> 0
1 -> 1
2 -> sqrt(2)
(you can't make 3 with a sum of 2 squares)
4 -> 2
5 -> sqrt(5)
etc.

You can either leave gaps in the table (will probably be quite large then) or else do a binary search looking for your value since the left hand side sequence is strictly increasing). EDIT: If you have gaps in the table you might save more space overall though since you don't have to store the key values then (sequence on the left above)

To find out what entries you would not need in the lookup table you can use this theorem to check you haven't made a mistake

https://en.wikipedia.org/wiki/Sum_of_tw ... es_theorem

I don't really like how they describe the result there though

To find out if an integer N can be expressed as a sum of 2 squares:

1) Factorise N. Get wolfram alpha to do that for you
2) If any of the odd prime factors (i.e. ignore factors of 2) have remainder 3 when divided by 4, they must occur to an even power

Example:

1234567890 factorises as

2 x 3 x 3 x 5 x 3607 x 3803 = 2 * 3^2 * 5 * 3607 * 3803

factor of 3 occurs to an even power, that's ok

3607 / 4 = 901 remainder 3. That's an odd power -> 1234567890 is NOT a sum of 2 squares

Less silly examples:

42 is not a sum of 2 squares, 42 = 2 * 3 * 7 and both 3 and 7 mod 4 = 3, so have odd powers in the prime factorisation

18 = 2 * 3^2 which although it has 3 as a factor and 3 mod 4 = 3, that occurs to an even power. So 18 is a sum of 2 squares (3^2 + 3^2).

I have a remarkable short proof of the result (well it's in a book I have), but this post is too small for me to fit it in

The case for whether a prime is a sum of 2 squares was proven by Fermat. A prime p is a sum of 2 squares if and only if p mod 4 != 3 (so 3 is not a sum of 2 squares, 5 is (2^2 + 1^2), 7 isn't, 11 isn't, 13 is (3^2 + 2^2), etc.

[ParadigmShifter]

Why work out the accessor function for a symmetric matrix when you can look one up

Index (M (i, j)) = (i * (i – 1) / 2) + (j-1)

But you'd need a multiply function to implement that (and can subtract 1 from the index if you don't store 0).

I feel some BASIC coding about to start lol

[sn3j]

But if we restrict the player (or NPCs) to 8-way movements only, the cost would always be m + n*sqrt(2) with m,n being natural numbers, right?
I am meanwhile beginning to think that Euclidean distance may not be the ideal metric for tiled maps.
And how about hexagonals like mentioned by @Evil Genius. There shouldn't be any square roots required at all.

Regarding lookup tables, the sums might get into the 1000s with bigger maps. You'd have to spend a lot of time to do the factorisation and looking up the results.

[ParadigmShifter]

ZX Basic coding complete.

Symmetric matrix stored as a triangular array implementation. Watch out for the arrays being 1-based thing in basic and the accessor FN a(i,j) only works properly if i <= j, have to use FN a(j,i) in that case.



Output



The accessor function in assembler needs either a multiply function or you can just repeatedly add increasing numbers in a loop:

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

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

About a zillion different ways to prove that result

Note that 1 + n = 2 + (n-1) = 3 + (n-2) etc. though and there are n/2 pairs of those. Formula works for odd and even n though.

EDIT: That output is bobbins though there's a bug in my code, I hate basic lol Bug is in FN a(i,j) I believe, array contents look fine. It works for asize <= 3

[ParadigmShifter]

But if we restrict the player (or NPCs) to 8-way movements only, the cost would always be m + n*sqrt(2) with m,n being natural numbers, right?
I am meanwhile beginning to think that Euclidean distance may not be the ideal metric for tiled maps.
And how about hexagonals like mentioned by @Evil Genius. There shouldn't be any square roots required at all.

Regarding lookup tables, the sums might get into the 1000s with bigger maps. You'd have to spend a lot of time to do the factorisation and looking up the results.

You don't need to do any factorisation. That just explains why there are gaps in the table (why you can't make 3 as a sum of 2 squares so won't be in the table but 5 = 1^2 + 2^2 so is in the table).

You'll need a lot of space though yeah for that, especially if you leave gaps (but that might need less space than having a key/value pair and doing a bsearch for the answer).

values being m + n * sqrt(2) is not true. if dx = 1 and dx = 2 then dx^2 + dy^2 = 5 and sqrt(5) is not a multiple of sqrt(2) (sqrt(a) = k * sqrt(2) if and only if 2 divides a, again due to prime factorisation and sqrt(a*b) = sqrt(a)*sqrt(b)). Individual steps are always either distance 1 or distance sqrt(2) though of course.

Using the triangular array thing is probs best for space... my basic program is storing i^2 + j^2 in the array but you want to store sqrt(i^2+j^2) in the array instead. Then you won't need a lookup for the sqrt.