Buggy python submission (verdict: Runtime Error):
import sys
sys.setrecursionlimit(10**7) #再帰関数の上限,10**5以上の場合python
import math
from copy import copy, deepcopy
from copy import deepcopy as dcp
from operator import itemgetter
from bisect import bisect_left, bisect, bisect_right#2分探索
#bisect_left(l,x), bisect(l,x)#aはソート済みである必要あり。aの中からx未満の要素数を返す。rightだと以下
from collections import deque
#deque(l), pop(), append(x), popleft(), appendleft(x)
#q.rotate(n)で → にn回ローテート
from collections import Counter#文字列を個数カウント辞書に、
#S=Counter(l),S.most_common(x),S.keys(),S.values(),S.items()
from itertools import accumulate,combinations,permutations#累積和
#list(accumulate(l))
from heapq import heapify,heappop,heappush
#heapify(q),heappush(q,a),heappop(q) #q=heapify(q)としないこと、返り値はNone
#import fractions#古いatcoderコンテストの場合GCDなどはここからimportする
from functools import lru_cache#pypyでもうごく
#@lru_cache(maxsize = None)#maxsizeは保存するデータ数の最大値、2**nが最も高効率
from decimal import Decimal

def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printl(li): _=print(*li, sep="\n") if li else None
def argsort(s, return_sorted=False): 
    inds=sorted(range(len(s)), key=lambda k: s[k])
    if return_sorted: return inds, [s[i] for i in inds]
    return inds
def alp2num(c,cap=False): return ord(c)-97 if not cap else ord(c)-65
def num2alp(i,cap=False): return chr(i+97) if not cap else chr(i+65)
def matmat(A,B):
    K,N,M=len(B),len(A),len(B[0])
    return [[sum([(A[i][k]*B[k][j]) for k in range(K)]) for j in range(M)] for i in range(N)]
def matvec(M,v):
    N,size=len(v),len(M)
    return [sum([M[i][j]*v[j] for j in range(N)]) for i in range(size)]
def T(M):
    n,m=len(M),len(M[0])
    return [[M[j][i] for j in range(n)] for i in range(m)]
def main():
    mod = 1000000007
    #w.sort(key=itemgetter(1),reversed=True)  #二個目の要素で降順並び替え

    #N = int(input())
    N, b1,b2,b3 = map(int, input().split())
    #A = tuple(map(int, input().split())) #1行ベクトル
    #L = tuple(int(input()) for i in range(N)) #改行ベクトル
    L = tuple(tuple(map(int, input().split())) for i in range(N)) #改行行列
    R = tuple(tuple(map(int, input().split())) for i in range(N)) #改行行列

    A=[[L[i][j] for i in range(N)] for j in range(N)]
    aveL=0
    aveR=0
    for li, ri in zip(L,R):
        aveL+=sum(li)
        aveR+=sum(ri)
    aveL/=N**2
    aveR/=N**2
    tave=round((aveL+aveR)*0.5)

    for i in range(0,N):
        for j in range(0,N):
            kouho=set()
            l=L[i][j]
            r=R[i][j]
            
            for i2 in range(i-1):
                yoko=0
                for k in range(i2,i):
                    yoko+=A[k][j]
                    for b in (b1,b2,b3):
                        if l<=b-yoko<=r:
                            kouho.add((b-yoko,b))
            for j2 in range(j-1):
                tate=0
                for k in range(j2,j):
                    tate+=A[i][k]
                    for b in (b1,b2,b3):
                        if l<=b-tate<=r:
                            kouho.add((b-tate,b))
            if not kouho:
                if l<=tave<=r:
                    A[i][j]=tave
                elif tave<l:
                    A[i][j]=l
                else:A[i][j]=r
                continue
            opt=-10000000
            dif=1000000
            for (x,b) in kouho:
                cdif=-abs(x-tave)-b*10
                if cdif<dif:
                    opt=x
                    dif=cdif
            A[i][j]=opt

    for a in A:
        print(*a)
    f="./out.txt"

    with open(f,"w") as f1:
        for a in A:
            for ai in a:
                f1.write(str(ai)+" ")  
            f1.write("\n")





if __name__ == "__main__":
    main()

Fixed submission:
import sys
sys.setrecursionlimit(10**7) #再帰関数の上限,10**5以上の場合python
import math
from copy import copy, deepcopy
from copy import deepcopy as dcp
from operator import itemgetter
from bisect import bisect_left, bisect, bisect_right#2分探索
#bisect_left(l,x), bisect(l,x)#aはソート済みである必要あり。aの中からx未満の要素数を返す。rightだと以下
from collections import deque
#deque(l), pop(), append(x), popleft(), appendleft(x)
#q.rotate(n)で → にn回ローテート
from collections import Counter#文字列を個数カウント辞書に、
#S=Counter(l),S.most_common(x),S.keys(),S.values(),S.items()
from itertools import accumulate,combinations,permutations#累積和
#list(accumulate(l))
from heapq import heapify,heappop,heappush
#heapify(q),heappush(q,a),heappop(q) #q=heapify(q)としないこと、返り値はNone
#import fractions#古いatcoderコンテストの場合GCDなどはここからimportする
from functools import lru_cache#pypyでもうごく
#@lru_cache(maxsize = None)#maxsizeは保存するデータ数の最大値、2**nが最も高効率
from decimal import Decimal

def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printl(li): _=print(*li, sep="\n") if li else None
def argsort(s, return_sorted=False): 
    inds=sorted(range(len(s)), key=lambda k: s[k])
    if return_sorted: return inds, [s[i] for i in inds]
    return inds
def alp2num(c,cap=False): return ord(c)-97 if not cap else ord(c)-65
def num2alp(i,cap=False): return chr(i+97) if not cap else chr(i+65)
def matmat(A,B):
    K,N,M=len(B),len(A),len(B[0])
    return [[sum([(A[i][k]*B[k][j]) for k in range(K)]) for j in range(M)] for i in range(N)]
def matvec(M,v):
    N,size=len(v),len(M)
    return [sum([M[i][j]*v[j] for j in range(N)]) for i in range(size)]
def T(M):
    n,m=len(M),len(M[0])
    return [[M[j][i] for j in range(n)] for i in range(m)]
def main():
    mod = 1000000007
    #w.sort(key=itemgetter(1),reversed=True)  #二個目の要素で降順並び替え

    N = int(input())
    #N, K = map(int, input().split())
    A = tuple(map(int, input().split())) #1行ベクトル
    #L = tuple(int(input()) for i in range(N)) #改行ベクトル
    #S = tuple(tuple(map(int, input().split())) for i in range(N)) #改行行列
    ref=0
    for i in range(2,N):
        ref^=A[i]
    
    a0=A[0]
    a1=A[1]
    ans=0
    #print(bin(a0),bin(a1),bin(ref))
    INF=10**13
    @lru_cache(None)
    def dfs(a0,a1,ref):
        inf=INF
        if (a0 & 1) ^ (a1 & 1) != (ref & 1):
            return inf
        
        if a0<0:
            return inf
        if ref == 0:
            return inf if a0<a1 else (a0-a1)//2
        a01=a0&1
        a11=a1&1
        ref1=ref&1
        ans=inf
        if a01^a11==ref1:
            ans=min(ans,2*dfs(a0>>1,a1>>1,ref>>1))
        ans=min(ans,2*dfs((a0-1)//2,(a1+1)//2,ref>>1)+1)

        return ans
    x=dfs(a0,a1,ref)
    if x>=a0:
        print(-1)
    else:
        print(x)
    # for j in range(0,41):
    #     a0j=(a0>>j)&1
    #     a1j=(a1>>j)&1
        
    #     refj=(ref>>j)&1
    #     #print(a0j,a1j,refj)
    #     if a0j^a1j==refj:
    #         nref=(ref>>(j+1))&1

    #         na0=((a0>>(j+1))&1)
    #         na1=((a1>>(j+1))&1)
            
    #         if ((a0>>(j+1))&1)^((a1>>(j+1))&1)!=nref:
    #             #print("a",na0,na1,nref)
    #             shif=1<<j
    #             if a0>0 and (((a0-shif))^(((a1+shif))==ref:
    #                 ans+=shif
    #                 a0-=shif
    #                 a1+=shif


    #         continue
    #     else:
    #         shif=1<<j

    #         ans+=shif
    #         a0-=shif
    #         a1+=shif
    #         #print(shif,a0,a1)
    #         if a0<0:
    #             print(-1)
    #             return









if __name__ == "__main__":
    main()