Problem:
 Score : 800 points    Problem Statement Takahashi is doing a research on sets of points in a plane.
Takahashi thinks a set S of points in a coordinate plane is a good set when S satisfies both of the following conditions:   The distance between any two points in S is not \sqrt{D_1} .  The distance between any two points in S is not \sqrt{D_2} .   Here, D_1 and D_2 are positive integer constants that Takahashi specified.  Let X be a set of points (i,j) on a coordinate plane where i and j are integers and satisfy 0 ≤ i,j < 2N .  Takahashi has proved that, for any choice of D_1 and D_2 , there exists a way to choose N^2 points from X so that the chosen points form a good set.
However, he does not know the specific way to choose such points to form a good set.
Find a subset of X whose size is N^2 that forms a good set.      Constraints  1 ≤ N ≤ 300  1 ≤ D_1 ≤ 2×10^5  1 ≤ D_2 ≤ 2×10^5  All values in the input are integers.         Input Input is given from Standard Input in the following format:  N  D_1  D_2       Output Print N^2 distinct points that satisfy the condition in the following format:  x_1  y_1  x_2  y_2 : x_{N^2}  y_{N^2}   Here, (x_i,y_i) represents the i -th chosen point. 0 ≤ x_i,y_i < 2N must hold, and they must be integers.
The chosen points may be printed in any order.
In case there are multiple possible solutions, you can output any.        Sample Input 1 2 1 2      Sample Output 1 0 0
0 2
2 0
2 2  Among these points, the distance between 2 points is either 2 or 2\sqrt{2} , thus it satisfies the condition.       Sample Input 2 3 1 5      Sample Output 2 0 0
0 2
0 4
1 1
1 3
1 5
2 0
2 2
2 4

Solution (Python):
import sys
input = sys.stdin.readline

import numpy as np

N,*D = map(int,input().split())

grid = np.zeros((2*N,2*N), dtype=np.int8)

for i,d in enumerate(D):
    width = 1
    while d % 4 == 0:
        d //= 4
        width *= 2
    if d%4 == 1:
        x = np.arange(2*N)
        bl = (x%(2*width) >= width)
        grid += (bl[:,None]^bl[None,:]) * (1<<i)
    elif d%4 == 2:
        x = np.arange(2*N)
        bl = (x%(2*width) >= width)
        grid += bl * (1<<i)

for i in range(4):
    X,Y = np.where(grid == i)
    if len(X) >= N*N:
        break

X,Y = X[:N*N], Y[:N*N]
for x,y in zip(X,Y):
    print(x,y)

Solution (Python):
def check_odd(i, j):
    return (i + j) % 2 == 0


def check_even(i, j):
    return i % 2 == 0


def solve(n, d1, d2):
    s1, s2 = 0, 0
    while d1 % 4 == 0:
        d1 >>= 2
        s1 += 1
    while d2 % 4 == 0:
        d2 >>= 2
        s2 += 1

    f1 = check_odd if d1 % 2 else check_even
    f2 = check_odd if d2 % 2 else check_even
    lim = n ** 2

    buf = []
    cnt = 0
    for i in range(2 * n):
        _i1 = i >> s1
        _i2 = i >> s2
        for j in range(2 * n):
            if f1(_i1, j >> s1) and f2(_i2, j >> s2):
                buf.append('{} {}'.format(i, j))
                cnt += 1
                if cnt == lim:
                    return buf


n, d1, d2 = map(int, input().split())
print('\n'.join(solve(n, d1, d2)))

Solution (Python):
from itertools import product


def check_odd(i, j):
    return (i + j) % 2 == 0


def check_even(i, j):
    return i % 2 == 0


def solve(n, d1, d2):
    s1, s2 = 0, 0
    while d1 % 4 == 0:
        d1 >>= 2
        s1 += 1
    while d2 % 4 == 0:
        d2 >>= 2
        s2 += 1

    f1 = check_odd if d1 % 2 else check_even
    f2 = check_odd if d2 % 2 else check_even
    lim = n ** 2

    buf = []
    cnt = 0
    for i, j in product(range(2 * n), repeat=2):
        if f1(i >> s1, j >> s1) and f2(i >> s2, j >> s2):
            buf.append('{} {}'.format(i, j))
            cnt += 1
            if cnt == lim:
                break
    return buf


n, d1, d2 = map(int, input().split())
print('\n'.join(solve(n, d1, d2)))

Solution (Python):
from itertools import product

N, D1, D2 = map(int, input().split())


Ds = []
Ms = []
for D in [D1, D2]:
    M = 1
    while D % 4 == 0:
        D //= 4
        M *= 2

    Ds += [D]
    Ms += [M]


num = 0
for i, j in product(range(2 * N), repeat=2):
    
    for D, M in zip(Ds, Ms):
        if D % 2:
            
            if (i // M + j // M) % 2: break
        else:
            
            if (i // M) % 2: break
    else:
        print(i, j)
        num += 1
        if num == N * N: break

Solution (Python):


import sys,collections,itertools,re,math,fractions,decimal,random,array,bisect,heapq





def isqrt(n):
    """integer sqrt: the maximum integer less than or equal to sqrt(n)"""
    x = n
    y = (x + 1) // 2
    while y < x:
        x = y
        y = (x + n // x) // 2
    return x

def solve(f):
    n, d1, d2 = f.read_int_list()

    a = [[0] * (2*n) for _ in xrange(2*n)]
    d = [d1, d2]

    for i in xrange(2):
        dist = set()
        for dx in xrange(2*n):
            dy2 = d[i] - dx**2
            if dy2 < 0: continue
            dy = isqrt(dy2)
            if dy**2 != dy2: continue

            for sx in [1, -1]:
                for sy in [1, -1]:
                    dist.add((sx*dx, sy*dy))

        vis = [[False] * (2*n) for _ in xrange(2*n)]
        for x in xrange(2*n):
            for y in xrange(2*n):
                if not vis[x][y]:
                    c = 2**i
                    que = [(x, y)]
                    while len(que) > 0:
                        que2 = []
                        for xi, yi in que:
                            if vis[xi][yi]: continue
                            vis[xi][yi] = True
                            a[xi][yi] += c
                            for dx, dy in dist:
                                if 0 <= xi+dx < 2*n and 0 <= yi+dy < 2*n:
                                    que2.append((xi+dx, yi+dy))
                        que = que2
                        c *= -1

    cnt = collections.defaultdict(list)
    for x in xrange(2*n):
        for y in xrange(2*n):
            cnt[a[x][y]].append((x, y))

    for k, v in cnt.items():
        if len(v) >= n**2:
            return "\n".join(["%d %d"%item for item in v[:n**2]])


class Reader(object):
    def __init__(self, filename=None):
        self.file = open(filename) if filename is not None else None
        self.case = 1

    def __readline(self):
        return self.file.next().strip() if self.file else raw_input()

    def next_case(self):
        self.file.next()
        self.case += 1

    def read_int(self): return int(self.__readline())
    def read_float(self): return float(self.__readline())
    def read_long(self): return long(self.__readline())
    def read_decimal(self): return decimal.Decimal(self.__readline())
    def read_str(self): return self.__readline()

    def read_int_list(self): return map(int, self.__readline().split())
    def read_float_list(self): return map(float, self.__readline().split())
    def read_long_list(self): return map(long, self.__readline().split())
    def read_decimal_list(self): return map(decimal.Decimal, self.__readline().split())
    def read_str_list(self): return self.__readline().split()

if __name__ == '__main__':
    filename = sys.argv[1] if len(sys.argv) > 1 else None
    f = Reader(filename)
    if f.file:
        while True:
            print "Case 
            try:
                f.next_case()
            except StopIteration:
                break
    else:
        print solve(f)

Solution (Python):
def judge(D):
  n = 0
  while D%4==0:
    n += 1
    D //= 4

  return (lambda x,y: ~((x>>n)^(y>>n))&1) if D%2==1 else (lambda x,y: ~(x>>n)&1)

N,D1,D2 = map(int,input().split())
j1,j2 = judge(D1),judge(D2)
for _,(x,y) in zip(range(N*N),filter(lambda p: j1(*p) and j2(*p), ((x,y) for x in range(N*2) for y in range(N*2)))):
  print(x,y)

Solution (Python):
def judge(D):
  n = 0
  while D%4==0:
    n += 1
    D //= 4

  def h1(x,y):
    return ~((x>>n)^(y>>n))&1
  def h2(x,y):
    return ~(x>>n)&1

  return h1 if D%2==1 else h2

N,D1,D2 = map(int,input().split())
j1,j2 = judge(D1),judge(D2)

cnt = 0
for x in range(N*2):
  for y in range(N*2):
    if j1(x,y) and j2(x,y):
      print(x,y)
      cnt += 1
      if cnt >= N*N:
        exit()

Solution (Python):
from math import sqrt
from itertools import product

import numpy as np

def judge(D):
  n = 1
  while D%4==0:
    n *= 2
    D //= 4

  def h1(x,y):
    return ~(x//n+y//n)%2
  def h2(x,y):
    return ~(x//n)%2

  return h1 if D%2==1 else h2

N,D1,D2 = map(int,input().split())
j1,j2 = judge(D1),judge(D2)
cnt = 0
for x,y in product(range(N*2),repeat=2):
  if j1(x,y) and j2(x,y):
    print(x,y)
    cnt += 1
    if cnt >= N*N:
      break

Solution (Python):
from math import sqrt
from itertools import product

import numpy as np

def helper(N,D):
  if D%2 == 1:
    grid = np.zeros((N,N),dtype=bool)
    grid[::2,:] = True
    grid[:,1::2] = ~grid[:,1::2]
  elif D%4 == 2:
    grid = np.zeros((N,N),dtype=bool)
    grid[::2,:] = True
  else:
    M = (N-1)//2+1
    D //= 4
    g = helper(M, D)
    grid = np.zeros((2*M,2*M),dtype=bool)
    grid[::2,::2] = g
    grid[1::2,::2] = g
    grid[::2,1::2] = g
    grid[1::2,1::2] = g
    grid = grid[:N,:N]
  return grid


N,D1,D2 = map(int,input().split())
M = N*2
g1 = helper(N*2,D1)
g2 = helper(N*2,D2)
ans = g1&g2


for x,y in np.transpose(np.nonzero(ans))[:N*N]:
  print(x,y)

Solution (Python):
from sys import exit, setrecursionlimit, stderr
from functools import reduce
from itertools import *
from collections import defaultdict
from bisect import bisect

def read():
  return int(input())

def reads():
  return [int(x) for x in input().split()]

N, D1, D2 = reads()
NN = 2 * N

def deg(x):
  return next(i for i in count() if x >> i & 1 == 1)

p1 = deg(D1)
p2 = deg(D2)

def group(x, y, p):
  hp = p >> 1
  if p & 1 == 0:
    return (x >> hp & 1) ^ (y >> hp & 1)
  else:
    return x >> hp & 1

pts = [(x, y) for x in range(NN) for y in range(NN) if group(x, y, p1) == 0 and group(x, y, p2) == 0]
for x, y in pts[:N**2]:
  print(x, y)

Solution (Python):
from sys import exit, setrecursionlimit, stderr
from functools import reduce
from itertools import *
from collections import defaultdict
from bisect import bisect

def read():
  return int(input())

def reads():
  return [int(x) for x in input().split()]

N, D1, D2 = reads()
NN = 2 * N

def deg(x):
  return next(i for i in count() if x >> i & 1 == 1)

def group(x, y, d):
  p = deg(d)
  hp = p >> 1
  if p & 1 == 0:
    return (x >> hp & 1) ^ (y >> hp & 1)
  else:
    return x >> hp & 1

pts = [(x, y) for x in range(NN) for y in range(NN) if group(x, y, D1) == 0 and group(x, y, D2) == 0]
for x, y in pts[:N**2]:
  print(x, y)