Arjen Markus (13 january 2006) Suppose you want to examine the properties of points on a two-dimensional grid (or more-dimensional for that matter). To give an example: to find all the points of a curve f(x,y) = 0 with integer x and y (that problem is known as solving Diophantine equations). You can do that with a double loop:
for { set j 0 } { $j < $jmax } { incr j } {
for { set i 0 } { $i < $imax } { incr i } {
... check the properties ...
}
}The problem is that you need to specify the upper limits and that you will spend a lot of time in regions far from the origin, whereas in many cases, you will be interested mostly in points near the origin.
This becomes worse with increasing dimensions.
Is there an alternative? Yes, you could enumerate the points like this:
| .
| .
(5)--(8)--
| . .
| . .
(2)--(4)--(7)--
| . . .
| . . .
(0)--(1)--(3)--(6)--So that:
This scheme allows you to split a single integer into two integer coordinates. If you use the same procedure on either coordinate, you can enumerate the points on a three-dimensional grid and so on.
The advantages are that you only need a single loop, that you stick close to the origin all the time and that if you want more points, you simply expand that one loop. You can also easily select a slice away from the origin.
The script below illustrates this technique by enumerating the points with positive integer coordinates in the ball |x| < 5 in 4D space.
namespace eval ::IntegerGrid {
variable current 0
}
proc ::IntegerGrid::enum2d {} {
variable current
incr current
split2d $current
}
proc ::IntegerGrid::enum3d {} {
foreach {x yz} [enum2d] {break}
foreach {y z} [split2d $yz] {break}
return [list $x $y $z]
}
proc ::IntegerGrid::enum4d {} {
#
# To avoid double points in the origin
#
set xy 0
set yw 0
while { $xy == 0 || $zw == 0 } {
foreach {xy zw} [enum2d] {break}
}
foreach {x y} [split2d $xy] {break}
foreach {z w} [split2d $zw] {break}
return [list $x $y $z $w]
}
proc ::IntegerGrid::split2d {number} {
set ncount 1
while { $number > $ncount } {
set number [expr {$number-$ncount}]
incr ncount
}
set x [expr {$ncount-$number}]
set y [expr {$ncount-$x-1}]
list $x $y
}
# main --
# Enumerate the points in a ball |x| < 5 in 4D
# Estimating the space that we need to examine
# is trifle difficult, but take it wider than
# necessary to be sure
#
#
set rad 0
set count 0
while { $rad < 10 } {
foreach {x y z w} [::IntegerGrid::enum4d] {break}
set rad [expr {sqrt($x*$x+$y*$y+$z*$z+$w*$w)}]
if { $rad < 5 } {
puts "$x $y $z $w -- $rad"
incr count
}
}
puts "$count points found"
#
# Using the straightforward approach ...
#
set count 0
for { set x 0 } { $x < 5 } { incr x } {
for { set y 0 } { $y < 5 } { incr y } {
for { set z 0 } { $z < 5 } { incr z } {
for { set w 0 } { $w < 5 } { incr w } {
set rad [expr {sqrt($x*$x+$y*$y+$z*$z+$w*$w)}]
if { $rad < 5 } {
puts "$x $y $z $w -- $rad"
incr count
}
}
}
}
}puts "$count points found"
Lars H: On the issue of choosing a termination radius in the 4D example, the limit 10 is actually optimal (at least with respect to rational points in the ball; the integral points may allow a stricter bound, but it's not realistic to do such analysis in preparation for a simple enumeration such as this one). To see this, consider the point {2.5 2.5 2.5 2.5}. This is at distance 5 from the origin, so it lies on the boundary of the ball. The sum of its coordinates is 10, so no hyperplane $x+$y+$z+$w=$n with $n<10 falls outside the ball. Since [enum4d] walks through each such hyperplane before it proceeds to the next one, it is necessary to allow processing of $n=9 and thus also e.g. the point {9 0 0 0}. The least integral radius admitting this point is indeed 10.