PROGRAM gridmethod
! produce tiling of the plane using the so-called grid
! method of projecting the five dimensional hypercube onto two dimensions.
DIM xp(0 to 1,0 to 1),yp(0 to 1,0 to 1),kk(0 to 4),gamma(0 to 4)
DIM cj(0 to 4),sj(0 to 4)
DECLARE DEF delta, ceiling
ASK PIXELS px,py
LET ar = px/py
SET WINDOW -20*ar,20*ar,-20,20
INPUT prompt "Enter a small integer k = ? ": kmax
CLEAR
FOR j = 1 to 4
    LET gamma(j) = (rnd -.5)
    LET sum = sum + gamma(j)
NEXT j
LET gamma(0) = -sum               ! set sum of gamma(j) = 0
FOR j = 0 to 4
    LET cj(j) = cos(2*pi*j/5)     ! calculate sines and cosines
    LET sj(j) = sin(2*pi*j/5)
NEXT j
FOR  k1 = -kmax to kmax
    FOR k2 = -kmax to kmax
        FOR j1 = 0 to 3
            LET cj1 = cj(j1)
            LET sj1 = sj(j1)
            FOR j2 = j1 + 1 to 4
                LET cj2 = cj(j2)
                LET sj2 = sj(j2)
                ! find intersection of two grid lines
                LET y = ((k1-gamma(j1))*cj2
                let y = y  - (k2-gamma(j2))*cj1)/(sj1*cj2-sj2*cj1)
                LET x = (k2 - y*sj2 - gamma(j2))/cj2
                ! Find Kj(x,y)
                FOR j = 0 to 4
                    IF (j <> j1) and (j <> j2) then
                       LET kk(j) = ceiling(x*cj(j) + y*sj(j) + gamma(j))
                    ELSE IF (j = j1) then
                       LET kk(j) = k1
                    ELSE
                       LET kk(j) = k2
                    END IF
                NEXT j
                FOR e1 = 0 to 1
                    FOR e2 = 0 to 1
                        LET xp(e1,e2) = 0
                        LET yp(e1,e2) = 0
                        FOR j = 0 to 4
                            LET n = kk(j) + e1*delta(j,j1) + e2*delta(j,j2)
                            LET xp(e1,e2) = xp(e1,e2) + n*cj(j)
                            LET yp(e1,e2) = yp(e1,e2) + n*sj(j)
                        NEXT j
                    NEXT e2
                NEXT e1
                PLOT LINES: xp(0,0),yp(0,0);xp(0,1),yp(0,1);xp(1,1),yp(1,1)
                PLOT LINES: xp(1,1),yp(1,1);xp(1,0),yp(1,0);xp(0,0),yp(0,0)
            NEXT j2
        NEXT j1
    NEXT k2
NEXT k1
END
 
FUNCTION delta(i,j)
    IF i = j then
       LET delta = 1
    ELSE
       LET delta = 0
    END IF
END DEF
 
FUNCTION ceiling(x)
    LET ix = truncate(x,0)
    IF x > 0 then LET ix = ix + 1
    LET ceiling = ix
END DEF