首先分析一下这个鬼畜的函数，我们考虑
\(f(x)+2C\)
\(=f(2f(x)-x+1)+C\)
\(=f(2f(2f(x)-x+1)-(2f(x)-x+1)+1)\)
\(=f(2(f(x)+C)-2f(x)+x-1+1)\)
\(=f(x+2C)\)
也就是 \(f(x)=f(x\bmod 2C)+2C\lfloor \dfrac{x}{2C}\rfloor\)
也就是，只要决定了 \(f(x)\)，\(f(x+2mC)\) 也就被确定了。
但是剩下的就没有关系了吗？
设 \(a\) 是偶数，\(a=2a'\)
那么设 \(b'=(f(a)-a')\bmod C\)
也就是 \(f(a)=a'+b'+mC\)
那么对于 \(b=2b'+1\)
\(f(2f(a)-a+1)=f(a)+C\)
\(f(2a'+2b'-a+2mC+1)=a'+b'+mC+C\)
\(f(b)+2mC=a'+b'+mC+C\)
\(f(b)=a'+b'-mC+C\)
所以，一旦我们决定了 \(a\) 其实也就相对应的决定了 \(b\)。同时，如果我们决定 \(b\)，那么也可以从 \(b\) 反推出 \(a\)，也就是我们其实是将 \(a\) 和 \(b\) 配对了。
那么，如果我们将 \(a\) 和 \(b\) 配对，\(f(a)\) 和 \(f(b)\) 其实都变成了 \(m\) 的函数。而对于其他的值，它们之间的取值是互不相干的。
然后，我们对于所有 \(x_i\bmod 2C=a\)，
\(|f(x_i)-y_i|\)
\(=|f(a)+2C\lfloor \dfrac{x_i}{2C}\rfloor-y_i|\)
\(=|a'+b'+mC+2C\lfloor \dfrac{x_i}{2C}\rfloor-y_i|\)
\(m\) 对这个类的答案的贡献其实就是 \(mC\) 到 \(y_i-a'-b'-2C\lfloor \dfrac{x_i}{2C}\rfloor\) 的距离。
对于所有 \(x_i\bmod 2C=b\)，
\(|f(x_i)-y_i|\)
\(=|f(b)+2C\lfloor \dfrac{x_i}{2C}\rfloor-y_i|\)
\(=|a'+b'-mC+C+2C\lfloor \dfrac{x_i}{2C}\rfloor-y_i|\)
\(=|mC-a'-b'-C-2C\lfloor \dfrac{x_i}{2C}\rfloor+y_i|\)
\(m\) 对这个类的贡献就是 \(mC\) 到 \(a'+b'+C+2C\lfloor \dfrac{x_i}{2C}\rfloor-y_i\) 的距离。
那么，我们就可以把这些数依次排开，然后我们就希望 \(mC\) 能取到这些数的中位数。可惜有时取不到，那就把最靠近中位数的两个都拿出来试试，得到 \(a\) 和 \(b\) 配对情况下 \(a\) \(b\) 两类的最小贡献。
然后就是二分图最小权最大匹配，因为 \(C\) 是 \(16\)，所以用状压 \(dp\)、费用流、\(\text{KM}\) 等各种算法碾过去就可以了。
struct FunctionalEquation{
int x[10005],y[10005];ll e[16][16];
ll dp[1<<17];
ll minAbsSum(int c,int n,int xzero,int xprod,int xadd,int xmod,int yzero,int yprod,int yadd,int ymod){
x[0]=xzero,y[0]=yzero;
rep(i,1,n-1)x[i]=(1ll*x[i-1]*xprod+xadd)%xmod;
rep(i,1,n-1)y[i]=(1ll*y[i-1]*yprod+yadd)%ymod;
vt<int>v[32];
rep(i,0,n-1){
y[i]-=2*c*(x[i]/(2*c));
x[i]=x[i]%(2*c);
v[x[i]].pb(y[i]);
}
rd(i,2*c)sort(v[i].begin(),v[i].end());
for(int a=0;a<c;a++)for(int b=0;b<c;b++){
vt<int>cur;
for(auto i:v[2*a]){
cur.pb(i-a-b);
}
for(auto i:v[2*b+1]){
cur.pb(a+b+c-i);
}
sort(cur.begin(),cur.end());
int m=cur.size();
if(m==0){
e[a][b]=0;
continue;
}
m=(m-1)/2;
ll sum1=0,sum2=0,n1=floor(1.0*cur[m]/c),n2=n1+1;
for(auto i:cur)sum1+=abs(n1*c-i);
for(auto i:cur)sum2+=abs(n2*c-i);
e[a][b]=min(sum1,sum2);
}
rd(i,(1<<c))dp[i]=1e18;
dp[0]=0;
rep(msk,0,(1<<c)-1){
int x=__builtin_popcount(msk);
rd(y,c)if(!(msk>>y&1)){
dp[msk|(1<<y)]=min(dp[msk|(1<<y)],dp[msk]+e[x][y]);
}
}
return dp[(1<<c)-1];
}
};