# Create image of Ulam spiral

SIZE = 201    # image size (must be odd)

P = [0] * (SIZE * SIZE + 1)

for C in range(2, len(P)):
    if P[C] == 0:
        #print(C)
        for x in range(2, len(P)):
            z = x * C
            if z >= len(P):
                break
            P[z] = 1

cdir = 0
dirs = [(1,0), (0,-1), (-1,0), (0,1)]
out = [["0"] * SIZE for y in range(SIZE)]
px, py = SIZE//2, SIZE//2
nstep = []
cstep = 1
for x in range(SIZE * SIZE):
    nstep.append(cstep)
    nstep.append(cstep)
    cstep += 1

cstep = nstep.pop(0)
for p in range(SIZE * SIZE):
    if p > 0 and P[p+1] == 0:
        out[py][px] = "1"
    px += dirs[cdir % 4][0]
    py += dirs[cdir % 4][1]
    cstep -= 1
    if cstep == 0:    
        cdir += 1
        cstep = nstep.pop(0)

# Save image in NetPBM bitmap format
f = open("ulam.pbm", "w")
f.write(f"P1\n{SIZE} {SIZE}\n")
for y in range(SIZE):
    f.write(''.join(out[y]) + "\n")
f.close()

# Convert PBM to PNG with ImageMagick (and increase resolution 10x):
#  mogrify -format png -interpolate Nearest -filter point -resize 1000% ulam.pbm