Arithmetization in STARKs (第一part作业）！

[Xor0v0]

Add-on

Stark101课程中前2节主要叙述整个课程所需要的基本组件和算术化过程，为后续FRI打基础。非常抱歉分享的时候出现点小问题，利用 Singleton bound证明RS code的最小汉明距离时脑子有点抽，应该是次数为 $k-1$ 的多项式，这里重新推理一下：

考虑RS码的最小汉明距离： $d=\Delta(u,v)$。其实是长度为 n 的两个多项式「值向量」的差异。由于RS码的消息（多项式） $p(X)$ 的次数 $\le k-1$，根据代数基本定理我们知道，其在 $\mathbb{F}_p$ 上的零点个数也要 $\le k-1$。那么在大小为 $n$ 的有限集合上取值，非零点的个数 $\ge n-k+1$。于是码字的汉明重量 $w\ge n-k+1$。又由于线性码的最小汉明距离和最小汉明重量是等价的，所以码字的汉明最小汉明距离 $d\ge n-k+1$。而根据 Singleton Bound $[n,k]$-线性码的最小汉明距离 $d\le n-k+1$，因此 RS 码的最小汉明距离等于 $n-k+1$。

（分享的时候误说成 $\le k$ 次多项式，因此少了个 + 1。非常抱歉！

这是分享的内容笔记，希望共学的同志和老师们指正，非常感谢！

Quiz

卡特兰数（Catalan number）是一种经典的组合数学问题，它在计算机科学、统计学和物理学等领域中有着广泛的应用。

它的其中一种定义为：有一个大小为 $n\times n$ 方格图左下角为 (0, 0)， 右上角为 n, n，从左下角开始每次都只能向右或者向上走一单位，不走到对角线上方（但可以触碰）的情况下到达右上角有多少可能的路径？

递推公式为：

$$ C_0=1, C_{n+1}=\frac{2(2n+1)}{n+2}C_n $$

Hint: 注意上述递推公式的 n ，仅在自然数取值时有效。当插值到多项式 $f(X)$ 的时候定义域使用的是乘法子群元素，而 $\frac{2(2n+1)}{n+2}$ 部分仍然还是自然数取值，因此此处需要特别处理。

现在，在STARK101 toturial 同样的有限域设置下，Prover 想给 Verifier 证明卡特兰数列第100项的数字为1555040254。请您对上述计算约束做算术化。

文字说明思路，以及给出最终 composition poly 在 $h^{111}$ 上的取值 $p(h^{111})$ （或者吉利点 $p(h^{888})$ 也行） 作为 proof of work。（小域 $\langle g\rangle$ size取 128，大域 $\langle h\rangle$ size取 1024）

[neuroney]

Clearly, $p_0(x)=\frac{f(x)-1}{x-g^0}$ and $p_1(x)=\frac{f(x)-1555040254}{x-g^{100}}$.

Note that $C_{n+1}=\frac{2(2n+1)}{n+2}C_{n} \Leftrightarrow C_{n+1}=\sum_{i=0}^{n}C_i \cdot C_{n-i}$,

$p_{2,1}(g^x)=\frac{f(g^{x+1})- \frac{(4x+2)f(g^x)}{x+2} }{\prod_{i=0}^{99}(x-g^i)}$ or

$p_{2,2}(g^x)=\frac{f(g^{x+1})- \sum_{i=0}^{x} f(g^i)\cdot f(g^{x-i}) }{\prod_{i=0}^{99}(x-g^i)}$

If $|G|=128$, $x^{128}-1=\prod_{i=0}^{127}(x-g^i)$, then $\prod_{i=0}^{99}(x-g^i)=\frac{x^{128}-1}{\prod_{i=100}^{127}(x-g^i)}$.

If we use p_{2,1}(g^x), ${\rm CP}(h^{888})={\rm CP}(g^{111})=962534416$.
If we use p_{2,2}(g^x), ${\rm CP}(h^{888})={\rm CP}(g^{111})=1523195286$.

[nishuzumi]

from field import FieldElement
from channel import Channel
from polynomial import interpolate_poly, X,prod

公式： $$ C_0=1, C_{n+1}=\frac{2(2n+1)}{n+2}C_n $$

a = [FieldElement(1),FieldElement(1)]  # C(1) = 1
for i in range(1,100):
    numerator = FieldElement(2) * (FieldElement(2)*i + 1)
    denominator = i + 2
    a.append((numerator * a[-1]) / denominator)

g = FieldElement.generator() ** ((3*2**30) // 128)
G = [g**i for i in range(128)]
f = interpolate_poly(G[:101], a)
assert f(g**100) == 1555040254

w = FieldElement.generator()
h = w ** ((2 ** 30 * 3) // 1024)
H = [h ** i for i in range(1024)]
eval_domain = [w * h for h in H]
f_eval = [f(d) for d in eval_domain]

p0 = (f - 1) / (X - 1)
p1 = (f - 1555040254) / (X - g**100)

M_points = [(g**i, FieldElement(i)) for i in range(100)]
M = interpolate_poly([p[0] for p in M_points], [p[1] for p in M_points])

def catalan_constraint(x):
  return f(g * x) * (M + 2) - (FieldElement(2)*(FieldElement(2)*M+1) * f(x))

p2_numerator = catalan_constraint(X)
p2_denominator = prod([X - g**i for i in range(100)])
p2 = p2_numerator / p2_denominator

def get_CP(channel):
    alpha0 = channel.receive_random_field_element()
    alpha1 = channel.receive_random_field_element()
    alpha2 = channel.receive_random_field_element()
    # ignore alpha0 because it will always be 0
    return p0 + alpha1*p1 + alpha2*p2

test_channel = Channel()
CP_test = get_CP(test_channel)
CP_test(h**888)

595485781

[tmors]

请问下：
alpha0 = channel.receive_random_field_element()
这个好像随机出来一直是0，组合多项式的时候允许alpha为0吗？

[justinzjj]

请问为什么这个M的映射为什么是成立的啊 就是为什么是这样映射的

[justinzjj]

请问为什么这个M的映射为什么是成立的啊 就是为什么是这样映射的

我懂了 其实只要是任意一个可以将元素映射到自然数0-100之间的映射关系就都是ok的 这个是最直观简单的 方式

[nishuzumi]

请问下： alpha0 = channel.receive_random_field_element() 这个好像随机出来一直是0，组合多项式的时候允许alpha为0吗？

我感觉，好像是有点不对劲

[Xor0v0]

理论上可以为0，但是这样会损耗掉第一个多项式的贡献，不可取。比较合适的做法是：在做三个多项式的 linear random combination 时第一个数一般就取 1，然后FS 生成 2 个随机数。这个repo把第一个值取到0，应该是非预期的行为。

[wz14]

Select a public ploynomial $ID(g^x) = x$ for $x$ in $[0, 127]$.
For Catalan number, there are three constraints for the trace. Suppose $f$ is a trace ploynomial.

1. $f(g^0) = 1$
2. $f(g^{100}) = 1555040254$
3. $( ID(x) + 2 ) * f(g * x) - f(x) * (4 * ID(x) + 2) = 0$

These contrait can be expressed by quotient ploynomials.

1. $p_0 = (f(x)-1) / (x-g**0)$
2. $p_1 = (f(x)- 1555040254) / (x - g^{100})$
3. $p_2 = ((ID(x) + 2 ) * f(g * x) - f(x) * (4 * ID(x) + 2)) / ((x^{128} - 1) / (x - g^{127}) )$, here for $g^{127}$ should not satisfy this constraint.

$CP$ is a combination of $p_0$, $p_1$ and $p_2$.

Here is the code.

t = [FieldElement(1)]
for i in range(127):
    n = i
    t.append(t[-1] * (4 * FieldElement(n) + 2) / (FieldElement(n) + 2))
print(t[0], t[1], t[2])
print(t[100].val)  # 1555040254
g = FieldElement.generator() ** (3 * 2**23)  # for 128 points
points = [g**i for i in range(2**7)]
f = interpolate_poly(points, t)
idx = [FieldElement(i) for i in range(2**7)]  # ploy(g^n) -> n
i = interpolate_poly(points, idx)

# constrait 1
p0 = (f - 1) / (X - g**0)

# constrait 2
p1 = (f - 1555040254) / (X - g**100)

# constriat 3
numer = f * (4 * i(X) + 2) - f(g * X) * (i(X) + 2)
print(i(g).val)
print(f(g).val)
print(f(g**2).val)
p2 = numer / ((X**128 - 1) / (X - g**127))

print("deg p0 =", p0.degree())
print("deg p1 =", p1.degree())
print("deg p2 =", p2.degree())

h_gen = FieldElement.generator() ** ((2**30 * 3) // 1024)  # extend to 1024
h = [h_gen**i for i in range(1024)]
domain = [FieldElement.generator() * x for x in h]
ev = [f.eval(d) for d in domain]
mt = MerkleTree(ev)

channel = Channel()
channel.send(mt.root)

alpha0 = channel.receive_random_field_element()
alpha1 = channel.receive_random_field_element()
alpha2 = channel.receive_random_field_element()
CP = alpha0 * p0 + alpha1 * p1 + alpha2 * p2

print(CP(h_gen**111).val)

[sinajia]

最开始，将卡特兰数的递推公式转化为适合 STARK 协议的算术约束。接着，构建多项式并计算 composition polynomial 在指定点的值。
步骤如下：

算术化约束：
将递推公式 Cn+1 = 2(2n+1)/(n+2) * Cn 转化为多项式关系。因为在有限域中工作，除法操作需要转换为乘法。
设 f(X) 为插值多项式，对应卡特兰数序列。一定有如下约束：
f(g^(i+1)) * (i+2) = 2(2i+1) * f(g^i)
其中 g 是小域 ⟨g⟩ 的生成元，i 是从0到99的索引。
边界条件：
设置 f(1) = 1 作为初始条件（C0 = 1）。
同时，必须确保 f(g^100) = 1555040254 作为最终结果。
构建 composition polynomial：
p(X) = (f(Xg) * (log_g(X) + 2) - 2(2log_g(X) + 1) * f(X)) *
(1 - (X-1)(X-g^100)) +
(f(1) - 1) * (X-g) * (X-g^2) * ... * (X-g^100) +
(f(g^100) - 1555040254) * X * (X-g) * ... * (X-g^99)
其中 log_g(X) 是在有限域中的离散对数函数。
计算 p(h^111) 或 p(h^888)：
这一步需要在大域 ⟨h⟩ 中进行具体的数值计算。一下是应该有的步骤：
a) 计算 g 和 h 在各自域中的具体值
b) 计算 f(X) 在所有需要的点上的值
c) 计算 log_g(X) 在需要的点上的值
d) 将所有这些值代入 composition polynomial 的表达式
e) 在指定点（h^111 或 h^888）评估最终的多项式

下面是代码

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, Polynomial, X

def generate_catalan_sequence(n):
    c = [1]
    for i in range(n):
        c.append(c[-1] * 2 * (2*i + 1) // (i + 2))
    return c

def setup_domain(size, subgroup_size):
    generator = FieldElement.generator()
    g = generator ** (3 * 2**20 // subgroup_size)
    return [g**i for i in range(size)]

def create_polynomial(points, values):
    return interpolate_poly(points, [FieldElement(v) for v in values])

def create_composition_polynomial(f, g, target_value):
    p0 = (f - 1) / (X - 1)
    p1 = (f - target_value) / (X - g**100)
    
    M_x = [g**i for i in range(100)]
    M_y = [2*(2*FieldElement(i)+1)/(FieldElement(i)+2) for i in range(100)]
    M = interpolate_poly(M_x, M_y)
    
    p2 = (f(g*X) - M(X)*f(X)) / Polynomial.prod([X - g**i for i in range(100)])
    
    return p0, p1, p2

def get_composition_polynomial(channel, p0, p1, p2):
    return sum(channel.receive_random_field_element() * p for p in [p0, p1, p2])

def main():
    # Setup
    catalan = generate_catalan_sequence(128)
    g_domain = setup_domain(128, 128)
    h_domain = setup_domain(1024, 1024)
    
    # Create and verify polynomial
    f = create_polynomial(g_domain, catalan)
    assert f(g_domain[100]) == catalan[100]
    
    # Evaluate polynomial and create Merkle tree
    evaluations = [f.eval(d) for d in [FieldElement.generator()*x for x in h_domain]]
    merkle_tree = MerkleTree(evaluations)
    
    # Create composition polynomial
    p0, p1, p2 = create_composition_polynomial(f, g_domain[1], catalan[100])
    
    # Channel operations
    channel = Channel()
    channel.send(merkle_tree.root)
    CP = get_composition_polynomial(channel, p0, p1, p2)
    
    # Test results
    print('deg CP:', CP.degree())
    print('CP(h^111):', CP(h_domain[111]))
    print('CP(h^888):', CP(h_domain[888]))

if __name__ == "__main__":
    main()

[dyzz]

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, Polynomial, X, prod

$C_0=1$

$C_{n+1}=\frac{2(2n+1)}{n+2}C_n$

c=[1]
t=[FieldElement(1)]
for i in range(127):
    c.append(c[-1]*2*(2*i+1)//(i+2))
    t.append(FieldElement(c[-1]))

assert t[100]==1555040254

g = FieldElement.generator() ** (3 * 2**20)
points = [ g**i for i in range(128) ]
f = interpolate_poly(points, t)
assert f(g**100)==t[100]

h_gen = FieldElement.generator() ** ((3 * 2**20)//1024)
h = [ h_gen**i for i in range(1024) ]
domain = [FieldElement.generator()*x for x in h]
ev = [f.eval(d) for d in domain]
mt = MerkleTree(ev)
ch = Channel()
ch.send(mt.root)

$p_0=\frac{f-1}{X-g^0}$

numer0 = f - 1
denom0 = X - 1
p0 = numer0 / denom0

$p_1=\frac{f-1555040254}{X-g^{100}}$

numer1 = f - 1555040254
denom1 = X - g**100
p1 = numer1 / denom1

$f(gx)-\frac{2*(2*i+1)}{i+2}*f(x) = 0, for\ x=g^i,0\leq i\leq100$

# find a function M to remove `i` in the constraint
M_x = points[0:100]
M_y = [2*(2*FieldElement(i)+1)/(FieldElement(i)+2) for i in range(100)]
M = interpolate_poly(M_x, M_y)

$f(gx)-M(x)*f(x) = 0, for\ x=g^i,0\leq i\leq100$

$p_2 = \frac{f(gx)-M(x)*f(x)}{\displaystyle{\prod_{0}^{100}(X-g^i)}} $

numer2 = f(g*X)-M(X)*f(X)
denom2 = prod([(X - g**i) for i in range(100)])
p2 = numer2 / denom2

print('deg p0:', p0.degree())
print('deg p1:', p1.degree())
print('deg p2:', p2.degree())

deg p0: 126
deg p1: 126
deg p2: 126

def get_CP(channel):
    alpha0 = channel.receive_random_field_element()
    alpha1 = channel.receive_random_field_element()
    alpha2 = channel.receive_random_field_element()
    return alpha0*p0 + alpha1*p1 + alpha2*p2

test_channel = Channel()
CP_test = get_CP(test_channel)

print('deg CP_test:', CP_test.degree())
print('CP_test(111):', CP_test(h_gen**111))
print('CP_test(888):', CP_test(h_gen**888))

deg CP_test: 126
CP_test(111): 15337772
CP_test(888): 296400549

[dyzz]

https://github.com/readygo67/stark101?tab=readme-ov-file#issue-in-original-stark101

原tutorial的例子test_channel没有初始化，alpha0=0，应该使用之前定义的ch

def get_CP(channel):
    alpha0 = channel.receive_random_field_element()
    alpha1 = channel.receive_random_field_element()
    alpha2 = channel.receive_random_field_element()
    print('alpha0:', alpha0)
    print('alpha1:', alpha1)
    print('alpha2:', alpha2)
    return alpha0*p0 + alpha1*p1 + alpha2*p2

CP_test = get_CP(ch)

print('deg CP_test:', CP_test.degree())
print('CP_test(111):', CP_test(h_gen**111))
print('CP_test(888):', CP_test(h_gen**888))

alpha0: -1332874275
alpha1: 965077154
alpha2: 626357484
deg CP_test: 126
CP_test(111): -848946994
CP_test(888): 324484119

[geyu210]

我不是数学专业，基础太差，不了解很多名词需要恶补，最初是看不懂题目的，好在结合几位高人已经交出的作业，以及不断回看过去的几节课，加上把 stark101 给的 jupyterLab 例题过了几遍，终于搞明白题目在问什么，以及为什么要用这些步骤，收获很大！

由于水平限制，一定有不少错误以及概念的理解不到位，请老师们见谅

首先导入所需库，这些库是 stark101 仓库中 Stark101-part1.ipynb 与 Stark101-part2.ipynb 导入过的库，也是后面编程中需要用到的。

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, Polynomial, X, prod

卡特兰数的递推公式是：

$$ C_0 = 1 $$

$$ C_{n+1} = \frac{2(2n + 1)}{n + 2} C_n $$

由此可以递推出所有的值，由于后续会用到 g 的阶是 128 ，便按照卡特兰数公式按序计算出 128 个值写入一个 list ，最后验证了第 100 项值为 1555040254

a=[FieldElement(1)]
for i in range(127):
    a.append(a[-1]*(4* FieldElement(i)+2)/(FieldElement(i)+2))
assert a[100]==1555040254

将 a 列表在域 g 上转化成一个多项式，用到了拉格朗日插值定理库（除以 128 的目的是希望使 g 在乘以自己 128 次后回到 1 ）

$$ g = \text{generator}^{3 \times 2^{23}} $$

$$ \left( \text{generator}^{3 \times 2^{23}} \right)^{128} = g^{128} $$

$$ 2^7 = 128 $$

g = FieldElement.generator() ** (3 * 2**30//128)
print(g**128)
G = [ g**i for i in range(128) ]
f = interpolate_poly(G, a)
assert f(g**100)==a[100]

按教程验证阶数是否正确

# Checks that g and G are correct.
assert g.is_order(128), 'The generator g is of wrong order.'
b = FieldElement(1)
for i in range(127):
    assert b == G[i], 'The i-th place in G is not equal to the i-th power of g.'
    b = b * g
    assert b != FieldElement(1), f'g is of order {i + 1}'
if b * g == FieldElement(1):
    print('Success!')
else:
    print('g is of order > 128')

在 stark101 part2 中,我们需要将一个数列的三个约束转化为 三个多项式

# p0
p0 = (f - 1) / (X - g**0)
# p1
p1 = (f - 1555040254) / (X - g**100)

由于第三个多项式中有其他变量，需要用拉格朗日插值定理将数列转化为多项式

$$ f(gx) - \frac{2 \times (2 \cdot i + 1)}{i + 2} \cdot f(x) = 0, \text{ for } x = g^i, 0 \leq i \leq 100 $$

$$ f(gx) - M(x) \cdot f(x) = 0, \text{ for } x = g^i, 0 \leq i \leq 100 $$

$$ p_2 = \frac{f(gx) - M(x) \cdot f(x)}{\prod_{i=0}^{100} (x - g^i)} $$

# p2
## find a function M to remove `i` in the constraint

M_x = G[0:100]
M_y = [(4*FieldElement(i)+2)/(FieldElement(i)+2) for i in range(100)]
M = interpolate_poly(M_x, M_y)

p2 = (f(g*X)- M(X)*f(X)) / prod([(X- g**i) for i in range(100)])

最后，验证composition polynomials 在 ( h^{111} ) 上的取值 ( p(h^{111}) )

def get_CP(channel):
    alpha0 = channel.receive_random_field_element()
    alpha1 = channel.receive_random_field_element()
    alpha2 = channel.receive_random_field_element()
    return alpha0*p0 + alpha1*p1 + alpha2*p2

test_channel = Channel()
CP_test = get_CP(test_channel)

h = FieldElement.generator() ** (3 * 2**30//1024)
print('deg CP_test:', CP_test.degree())
print('CP_test(111):', CP_test(h**111))
print('CP_test(888):', CP_test(h**888))

#deg CP_test: 126
#CP_test(111): 1046493864
#CP_test(888): 1053646968

[justinzjj]

解题思路

和STARK101中的整体思路类似：

• 构建数列
• 选取小域g并构建f
• 选取大域h并扩张

构建数列：

# 1. get C
c = [FieldElement(1)]
for i in range(128):
	c.append(c[i] * 2 * (2 * i + 1) / (i + 2))
print("c[100] = " + str(c[100]))
# c[100]=1555040254

选取小域g并构建f

# 2. get G, G size is 128
g = FieldElement.generator() ** (3 * 2**30 / 128)
G = []
for i in range(128):
	G.append(g**i)

# 3. get f
f = interpolate_poly(G, c[:-1])

选取大域h并扩张

# 4. get H, H size is 1024
w = FieldElement.generator()
h = w ** ((3 * (2**30)) // 1024)
H = [h**i for i in range(1024)]
eval_domain = [w * x for x in H]

在这里因为引入了额外的变量，所以需要通过一个function 对其进行一个映射，映射到g内。
然后我们再对多项式进行约束即可

# 5. ploy(g^n) -> n
T = interpolate_poly([g**i for i in range(128)], [FieldElement(p) for p in range(128)]) 

# 6. constrait p0,p1,p2
numer0 = f - 1
denom0 = X - 1
p0 = numer0 / denom0
numer1 = f - 1555040254
denom1 = X - g**100
p1 = numer1 / denom1
numer2 = f(g * X) * (T(X) + 2) - (2 * (2 * T(X) + 1) * f(X))
denom2 = (X**128 - 1) / (X - g**127)
p2 = numer2 / denom2

最后打印结果即可

def get_cp(channel, p0, p1, p2):
	alpha0 = channel.receive_random_field_element()	
	alpha1 = channel.receive_random_field_element()	
	alpha2 = channel.receive_random_field_element()	
	return alpha0 * p0 + alpha1 * p1 + alpha2 * p2

channel = Channel()
cp = get_cp(channel, p0, p1, p2)

print("cp(h*888) = " + str(cp(h**888)))
print("cp(g*111) = " + str(cp(g**111)))

[0x-stan]

$p_0(x) = \frac{f(x) - 1}{x - g^0}$

$p_1(x) = \frac{f(x) - 1555040254}{x - g^{100}}$

$p_2(g^x)=\frac{f(g^x + 1) - \frac{2*(2x + 1)f(g^x)}{x + 2}}{\prod_{i=0}^{99} (x - g^i)}$

其中 $p_2$ 的 $\frac{2*(2x + 1)}{x + 2}$ 特殊处理，将自然数其映射到有限域，（由点集合插值得到多项式 M ）

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, X, prod

# 小域 g size 取 128，大域 h size 取 1024
g = FieldElement.generator() ** (3 * 2**30 // 128)
G = [g ** i for i in range(128)]

w = FieldElement.generator()
h = w ** ((3 * 2**30) // 1024)
H = [h ** i for i in range(1024)]
eval_domain = [w * x for x in H]

f = interpolate_poly(G[:101], a)
f_eval = [f(d) for d in eval_domain]
mt = MerkleTree(f_eval)

test_channel = Channel()
test_channel.send(mt.root)

assert f(g**100) == 1555040254

# constrait 1
numer0 = f - 1
denom0 = X - 1
p0 = numer0 / denom0

# constrait 2
numer1 = f - 1555040254
denom1 = X - g**100
p1 = numer1 / denom1

# constrait 3

# 将自然数其映射到有限域，（由点集合插值得到多项式 M ）
M = interpolate_poly(
        [g**i for i in range(100)],
        [2*(2*FieldElement(i)+1)/(FieldElement(i)+2) for i in range(100)])

p2_numerator = f(g * X) - M(X) * f(X)
p2_denominator = prod([X - g**i for i in range(100)])
p2 = p2_numerator / p2_denominator

def get_CP(channel):
    alpha0 = channel.receive_random_field_element()
    alpha1 = channel.receive_random_field_element()
    alpha2 = channel.receive_random_field_element()
    return alpha0 * p0 + alpha1*p1 + alpha2*p2
    
CP_test = get_CP(test_channel)
assert CP_test(g**111) == CP_test(h**888)
print(f'g**111 {CP_test(g**111).val}')
print(f'h**888 {CP_test(h**888).val}')

g^111 723710509
h^888 723710509

[sunjianping-alt]

//拷贝 @0x-stan程序，并尝试修改
from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, X, prod

//小域 g size 取 128，大域 h size 取 1024
g = FieldElement.generator() ** (3 * 2**27)
G = [g ** i for i in range(128)]

w = FieldElement.generator()
h = w ** (3 * 2**20）
H = [h ** i for i in range(1024)]
eval_domain = [w * x for x in H]

f = interpolate_poly(G[:100], a)
f_eval = [f(d) for d in eval_domain]
mt = MerkleTree(f_eval)

test_channel = Channel()
test_channel.send(mt.root)

assert f(g**99) == 1555040254

// constrait 1
numer0 = f - 1
denom0 = X - 1
p0 = numer0 / denom0

// constrait 2
numer1 = f - 1555040254
denom1 = X - g**99
p1 = numer1 / denom1

// constrait 3

// 将自然数其映射到有限域，（由点集合插值得到多项式 M ）
M1 = interpolate_poly(
[gi for i in range(99)],
(FieldElement(i)+2) for i in range(99)])
M2 = interpolate_poly(
[gi for i in range(99)],
[2*(2*FieldElement(i)+1) for i in range(99)])

p2_numerator = M1(x)*f(g * X) - M2(X) * f(X)
p2_denominator = prod([X - g**i for i in range(99)])
p2 = p2_numerator / p2_denominator

def get_CP(channel):
alpha= channel.receive_random_field_element()
beta = channel.receive_random_field_element()

return  p0 + alpha*p1 + beta*p2

CP_test = get_CP(test_channel)
assert CP_test(g111) == CP_test(h888)
print(f'g111 {CP_test(g111).val}')
print(f'h888 {CP_test(h888).val}')

[AmuroToru]

$$\begin{align} c_0=1 \rightarrow P_0=\frac{f(x)-1}{x-g^0}\\ c_{100}=1555040254 \rightarrow P_1=\frac{f(x)-1555040254}{x-g^{100}}\\ c_n= \frac{2*(2*n+1)}{n+2}c_{n+1}\rightarrow p_2= \frac{(f(gx)-M(x))}{\prod_{n=1}^{99}x-g^n}\\ \end{align}$$

其中$M(x)$是$g^i$和$\frac{2*(2*n+1)}{n+2}$的拉格朗日插值多项式

from field import FieldElement
from polynomial import interpolate_poly,X,prod
from channel import Channel

#construct Cattelan series
a = [FieldElement(1)]
t=[]
for n in range(127):
    a.append(FieldElement(2)*(FieldElement(2)*n+FieldElement(1))*FieldElement.inverse(n+FieldElement(2))*a[-1])
assert a[100]== 1555040254


#construct big field H 
c=FieldElement.generator()
h=c**((2**30*3)//1024)
H=[h**i for i in range(1024)]
eval_domain_h = [c * x for x in H]

#construct small field G
c=FieldElement.generator()
g=c**((3*2**30)//128)
G=[g**i for i in range(128)]
eval_domain_g=[c*x for x in G]


#get interpolate_polynomial
f = interpolate_poly(G, a)
assert f(G[100]==a[100])

#extend
f_eval = [f(d) for d in eval_domain_h]

#get p0 and p1 
p0 = (f - 1) / (X - 1)
p1 = (f - 1555040254) / (X - g**100)

#get p2

M_y = [FieldElement(2)*(2*FieldElement(i)+FieldElement(1))*FieldElement.inverse((FieldElement(i)+2)) for i in range(127)]
M = interpolate_poly(G[:-1], M_y)

numerator = f(g*X)-M(X)*f(X)
denominator =prod([X-g**n for n in range(127)])
p2 = numerator / denominator
channel=Channel()
#get CP(x)
alpha0 = channel.receive_random_field_element()
alpha1 = channel.receive_random_field_element()
alpha2 = channel.receive_random_field_element()
CP=alpha0*p0+alpha0*p1+alpha2*p2
print(CP(h**888))
#-501161392

[0xxEric]

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, Polynomial,X,prod

def part0():
# compute the number list
a=[1]
t=[FieldElement(1)]
for i in range(127):
a.append(2*(2*i+1)*a[i]//(i+2))
t.append(FieldElement(a[i+1]))
assert t[100]==1555040254, "Caculate wrong"
print("Building trace success!")

#construct the domain,poly
g = FieldElement.generator() ** ((FieldElement.k_modulus - 1)//128)
trace_domain = [g ** i for i in range(128)]
assert g**128==1
trace_poly = interpolate_poly(trace_domain, t)
assert trace_poly(g**100)==t[100]
h_gen = FieldElement.generator() ** ((FieldElement.k_modulus - 1) // 1024)
h = [h_gen ** i for i in range(8192)]
lde_domain = [FieldElement.generator() * x for x in h]

lde_eva=[trace_poly(d) for d in lde_domain]

xi= [FieldElement(i) for i in range(128)] 
interpolate_x=interpolate_poly(trace_domain, xi)

#construct the cp poly
p0=(trace_poly-FieldElement(1))/(X-1)
p1=(trace_poly-t[100])/(X-g**100)
p2=(interpolate_x+2)*trace_poly(g*X)-2*(2*interpolate_x+1)*trace_poly(X)
p2 = p2 / ((X**128 - 1) / (X - g**127))

mt = MerkleTree(lde_eva)
ch = Channel()
ch.send(mt.root)
alpha0 = ch.receive_random_field_element()
alpha1 = ch.receive_random_field_element()
alpha2 = ch.receive_random_field_element()
cp=alpha0*p0+alpha1*p1+alpha2*p2

#evaluate at g*w**111
cp_eva=cp(FieldElement.generator() *h_gen**111)
print(cp_eva)

[0xxEric]

得到的结果为：1359482188

[minionsjay]

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, Polynomial,prod
t=[FieldElement(1)]

生成卡特兰数

for i in range(0,100):
numerator = FieldElement(2) * (FieldElement(2)*i + 1)
denominator = i + 2
t.append((numerator * t[-1]) / denominator)

print("t:",t)

print("len t:",len(t))

print(FieldElement.k_modulus)

生成我们需要插值的x轴上的点的生成元，是一个乘法循环群，阶为128

g=FieldElement.generator()(3*2(23))
print("阶为128的乘法循环群生成元:",g)

print(FieldElement.generator())

points = [g**i for i in range(128)]

print("乘法循环群中的元素",points)

print(len(points[:-1]))

插值得到代数次数为 100的多项式

p = interpolate_poly(points[:101], t)

h_gen=FieldElement.generator()(3*2(20))
print("阶为1024的乘法循环群生成元:",h_gen)
h_points=[g**i for i in range(1024)]

print(h**1024)

生成陪集

domain = [FieldElement.generator() * x for x in h_points]

print(domain)

将证明转化为另一种形式

生成f(x)-1

numer0 = p - Polynomial([FieldElement(1)])

生成x-1

demo0 = Polynomial.gen_linear_term(FieldElement(1))

构造多项式 f(x)-1/(x-g^0)

q0, r0 = numer0.qdiv(demo0)

print("q0 degree:", q0.degree())

构造多项式 f(x)-1555040254

numer1 = p - Polynomial([FieldElement(1555040254)])

numer1

print(t[100])

print(numer1.degree())

构造多项式x-g^(100)

demo1 = Polynomial.gen_linear_term(points[100])

demo1

构造多项式 f(x)-2338775057/(x-g^1022)

q1, r1 = numer1.qdiv(demo1)
print("q1 degree:", q1.degree())

构造单项式g*x

inner_poly1 = Polynomial([FieldElement(0), points[1]])

计算f(g*x)

composition = p.compose(inner_poly1)

构造多项式 2*(2*x+1)/(x+2)

temp_x=points[0:100]

transition_x

temp_t=[2*(2*FieldElement(i)+1)/(FieldElement(i)+2) for i in range(0,100)]

temp_f=interpolate_poly(temp_x,temp_t)
print(temp_f.degree())

transition_f

构造 f(gx)-transition_ff(x)

T_f=composition-temp_f*p

T_f

构造多项式 (x-g^0)...(x-g^99)

factors=[Polynomial.gen_linear_term(points[i]) for i in range(0,100)]
prod_factors=prod(factors)

构造多项式f(gx)-2(2*x+1)/(x+2)*f(x)

pp,pr=T_f.qdiv(prod_factors)

pp=T_f/prod_factors

ch = Channel()
cp0 = q0.scalar_mul(ch.receive_random_field_element())
cp1 = q1.scalar_mul(ch.receive_random_field_element())
cp2 = pp.scalar_mul(ch.receive_random_field_element())
cp = cp0 + cp1 + cp2

print("p(h^{111})=",cp.eval(domain[111]))

print("p(h^{888})=",cp.eval(domain[888]))

[MarkYnx]

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, Polynomial, X, prod

t = [FieldElement(1)]
# calc catalan function
small_domain_size = 128
val_100 = 1555040254
for i in range(small_domain_size - 1):
    t.append(t[-1] * 2 * (FieldElement(2)*i+1) / (2+i))

assert t[100] == val_100
large_domain_size = 1024
# prepare small domain
g = FieldElement.generator() ** (3 * 2 ** 30//small_domain_size)
G = [g ** i for i in range(small_domain_size)]
# prepare large domain
h = FieldElement.generator() ** ((3 * 2 ** 30) // large_domain_size)
H = [h ** i for i in range(large_domain_size)]
domain = [FieldElement.generator() * x for x in H]
f = interpolate_poly(G, t)
ev = [f.eval(d) for d in domain]
assert f(g**100) == val_100

# make constraints
p0 = (f - 1) / (X - g**0) 
p1 = (f - val_100) / (X - g**100)

# special: make 'n' polynomial
N = [FieldElement(i) for i in range(100)]
fn = interpolate_poly(G[:100], N)
assert fn(g**99) == 99 

# make a transfer
# f(gx) = 2 * (2n+1) / (n + 2) * f(x) => (n + 2) * f(gx) = 2 * (2n+1) * f(x) 
# let n = fn(X)
numer  = (fn(X) + 2) * f(g*X) - (4*fn(X)+2) * f

lst = [(X - g**i) for i in range(100)]
denom = prod(lst) 
p2 = numer / denom

print("deg p0 = ", p0.degree())
print("deg p1 = ", p1.degree())
print("deg p2 = ", p2.degree())

# get  Composition Polynomial
mt = MerkleTree(ev)
ch = Channel()

ch.send(mt.root)
alpha0 = ch.receive_random_field_element()
alpha1 = ch.receive_random_field_element()
alpha2 = ch.receive_random_field_element()
print('alpha0 = ', alpha0)
print('alpha1 = ', alpha1)
print('alpha2 = ', alpha2)
cp = alpha0 * p0 + alpha1 * p1 + alpha2 * p2

print('deg cp:', cp.degree())
print('cp(g^111):', cp(g**111).val)
print('cp(h^888):', cp(h**888).val)
assert cp(g**111).val == cp(h**888).val
print('cp(h^111):', cp(h**111).val)

[MarkYnx]

cp(g^111): 2394400514
cp(h^888): 2394400514
cp(h^111): 370571673

[Draculabo]

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, X, prod

计算卡特兰数列

def calculate_catalan_sequence(size, val_100):
t = [FieldElement(1)]
for i in range(size - 1):
t.append(t[-1] * 2 * (FieldElement(2) * i + 1) / (i + 2))
assert t[100] == val_100
return t
c=[1]
t=[FieldElement(1)]
for i in range(127):
c.append(c[-1]2(2*i+1)//(i+2))
t.append(FieldElement(c[-1]))

assert t[100]==1555040254

准备域

def prepare_domain(size):
generator = FieldElement.generator()
g = generator ** (3 * 2 ** 30 // size)
return [g ** i for i in range(size)]

插值并求值

def interpolate_and_evaluate(G, t, domain):
f = interpolate_poly(G, t)
return f, [f.eval(d) for d in domain]

构造约束条件

def make_constraints(f, g, g0, g100, val_100):

构造第一个约束多项式

p0 = (f - 1) / (X - g0)

构造第二个约束多项式

p1 = (f - val_100) / (X - g100)

构造辅助多项式fn

N = [FieldElement(i) for i in range(100)]
fn = interpolate_poly(G[:100], N)
assert fn(g ** 99) == 99

构造第三个约束多项式

numer = (fn(X) + 2) * f(g * X) - (4 * fn(X) + 2) * f
alpha = prod([(X - g ** i) for i in range(100)])
p2 = numer / alpha

return p0, p1, p2, fn

获取组合多项式

def get_composition_polynomial(p0, p1, p2, ch):
ch.send(mt.root)

接收随机系数

alpha0 = ch.receive_random_field_element()
alpha1 = ch.receive_random_field_element()
alpha2 = ch.receive_random_field_element()

构造组合多项式

return alpha0 * p0 + alpha1 * p1 + alpha2 * p2

if name == "main":
small_domain_size = 128
large_domain_size = 1024
expect_value_100 = 1555040254

计算卡特兰数列

t = calculate_catalan_sequence(small_domain_size, expect_value_100)

准备小域G和大域H

G = prepare_domain(small_domain_size)
H = prepare_domain(large_domain_size)
domain = [FieldElement.generator() * x for x in H]

f, ev = interpolate_and_evaluate(G, t, domain)
assert f(G[100]) == expect_value_100

构造约束条件

p0, p1, p2, fn = make_constraints(f, G[1], G[0], G[100], expect_value_100)

打印约束多项式的度

print(f"deg p0 = {p0.degree()}")
print(f"deg p1 = {p1.degree()}")
print(f"deg p2 = {p2.degree()}")

构建Merkle树

mt = MerkleTree(ev)
ch = Channel()

获取组合多项式

cp = get_composition_polynomial(p0, p1, p2, ch)

打印结果

assert cp(G[111]).val == cp(H[888]).val
print(f'cp(g111): {cp(G[111]).val}')
print(f'cp(h888): {cp(H[888]).val}')
print(f'cp(h**111): {cp(H[111]).val}')

[ahy231]

思路是构造三个约束：

1. C_0 == 1
2. C_100 == 1555040254
3. C_(n+1) = 2 * (2 * n + 1) / (n + 2) * C_n

以下给出代码和每一步的注释：

# import field lib
from field import FieldElement
# import polynomial lib
from polynomial import interpolate_poly,X,prod
# import channel lib
from channel import Channel
# import merkle tree
from merkle import MerkleTree

# transcript
channel=Channel()

## trace

# initial state: C_0 = 1
a = [FieldElement(1)]
t=[]
# length of a must equals length of G
for n in range(128 - 1):
    # C_(n+1) = 2 * (2 * n + 1) / (n + 2) * C_n
    a.append(FieldElement(2)*(FieldElement(2)*n+FieldElement(1))*FieldElement.inverse(n+FieldElement(2))*a[-1])
assert a[100]== 1555040254

## domain

# generate large domain H of size 1024
c=FieldElement.generator()
h=c**((3*2**30)//1024)
H=[h**i for i in range(1024)]
eval_domain_h = [c * x for x in H]

# generate small domain G of size 128
c=FieldElement.generator()
g=c**((3*2**30)//128)
G=[g**i for i in range(128)]
eval_domain_g=[c * x for x in G]

## polynomial and evaluations

# interpolate a on G
f = interpolate_poly(G, a)
assert f(G[100]) == a[100] # check f(100) = a[100]

# LDE by RS encoding on coset of large domain H
f_eval = [f(d) for d in eval_domain_h]
eval_oracle = MerkleTree(f_eval) # commit f_eval
channel.send(eval_oracle.root) # write commmitment to transcript

## constraints

# constrain f(0) == 1 and f(100) == 1555040254
p0 = (f - 1) / (X - 1)
p1 = (f - 1555040254) / (X - g**100)

# constrain C_(n+1) = 2 * (2 * n + 1) / (n + 2) * C_n
M_y = [FieldElement(2)*(2*FieldElement(i)+FieldElement(1))*FieldElement.inverse((FieldElement(i)+2)) for i in range(127)]
M = interpolate_poly(G[:-1], M_y)

# quotient of p2
numerator = f(g*X)-M(X)*f(X)
denominator = (X ** 128 - 1) / (X - g ** 127) # due to (X - 1)(X - g)...(X - g ** 127) == (X ** 128 - 1)
p2 = numerator / denominator

# get CP(x)
alphas = [channel.receive_random_field_element() for _ in range(3)]
CP = alphas[0] * p0 + alphas[1] * p1 + alphas[2] * p2
print(CP(h**111))
# 831670508

[iopzhu]

$p_o(x)=\frac{f(x)-1}{x-g^0}$
$p_{1}(x)=\frac{f(x)-1555040254}{x-g^{100}}$
$p_2(x)=\frac{f(gx)-\frac{2*(2i+1)f(x)}{i+2}}{\prod_{i=0}^{99}(x-g^i)}$

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, X, prod

c=[1]
t=[FieldElement(1)]
for i in range(127):
    c.append(c[-1]*2*(2*i+1)//(i+2))
    t.append(FieldElement(c[-1]))

assert t[100]==1555040254

# 小域
g = FieldElement.generator() ** (3 * 2**30 // 128)
G = [g ** i for i in range(128)]

# 大域
w = FieldElement.generator()
h = w ** ((3 * (2**30)) // 1024)
H = [h**i for i in range(1024)]
eval_domain = [w * x for x in H]
mt = MerkleTree(eval_domain)
ch = Channel()
ch.send(mt.root)

# 构建函数
f = interpolate_poly(G, t)
f_eval = [f(d) for d in eval_domain]

# constrait p0
numer0 = f - 1
denom0 = X - 1
p0 = numer0 / denom0

# constrait p1
numer1 = f - 1555040254
denom1 = X - g**100
p1 = numer1 / denom1

# constrait p2
## AssertionError: Polynomials are not divisible.
# T = interpolate_poly([g**i for i in range(100)], [FieldElement(p) for p in range(100)])
# numer2 = f(g * X) - (2 * (2 * T(X) + 1) * f(X) / (T(X) + 2))
M = interpolate_poly(G[0:100], [2*(2*FieldElement(i)+1)/(FieldElement(i)+2) for i in range(100)])
numer2 = f(g * X) - M(X) * f(X)
denom2 = prod([(X - g**i) for i in range(100)])
p2 = numer2 / denom2

def get_CP(ch):
    alpha0 = ch.receive_random_field_element()
    alpha1 = ch.receive_random_field_element()
    alpha2 = ch.receive_random_field_element()
    return alpha0 * p0 + alpha1 * p1 + alpha2 * p2
CP_test = get_CP(ch)

print('deg CP_test:', CP_test.degree())
print('CP_test(111):', CP_test(h**111))
print('CP_test(888):', CP_test(h**888))

# deg CP_test: 126
# CP_test(111): -360732386
# CP_test(888): 913804724

[keroro520]

from field import FieldElement
import polynomial
from polynomial import X
from channel import Channel

# 1. Initialize the domain
#
# G = <g>, is a smooth multiplicative subgroup of size 128
# H = <h>, is a smooth multiplicative subgroup of size 1024
g = FieldElement.generator() ** (3 * 2**23)
h = FieldElement.generator() ** (3 * 2**20)
G = [ g**i for i in range(0, 2**7) ]

# 2. Calculate the sequence of Catalan Number
#
#    a[0] = 1
#    a[i+1] = a[i] * (2*(i+1)/(i+2))
catalan = [ FieldElement(1) ]
for i in range(0, 2**7-1):
    catalan.append( catalan[i] * FieldElement(2) * (FieldElement(2) * FieldElement(i) + FieldElement(1)) / FieldElement(i+2) )


# 3. Interpolate the polynomial _f_, using the X_axis as _G_ domain and Y_axis as the sequence of Catalan numbers
f = polynomial.interpolate_poly(G, catalan)

# 4. In order to map from _X_ to integer number, we define a polynomial integerize
#    integerize(g**0) == 0
#    integerize(g**1) == 1
#    integerize(g**2) == 2
#    integerize(g**n) == n
integerize = polynomial.interpolate_poly(G, [FieldElement(i) for i in range(128)])


# 5. Compose the CP using random linear combination
#
#    CP = alpha0 * (f(x) / (x-0))
#       + alpha1 * (f(x) / (x-100))
#       + alpha2 * (f(gx)*(n+2) / (f(x)*2*(n+1))) / ((x**128-1)/(x-127))
channel = Channel()
alphas = [ channel.receive_random_field_element() for _ in range(3) ]
CP =  alphas[0] * ((f - catalan[0])   / (X - g**0))     \
    + alphas[1] * ((f - catalan[100]) / (X - g**100))   \
    + alphas[2] * ((f(g * X) * (integerize(X) + 2) - f(X) * (2 * (2 * integerize(X) + 1))) / ((X**128 - 1) / (X-g**127)))


print('f.degree():'. f.degree())   # 127
print('CP.degree():', CP.degree()) # 127
print('random alphas:', alphas)    # [0, 787618507, -1067186547]
print('CP(h**111):', CP(h**111))   # -776067114
print('CP(h**888):', CP(h**888))   # -50623335

[Rikki424]

Note:

首先打开https://github.com/starkware-industries/stark101 安装相关库或者在线使用Jupyter Notebook.

个人不理解的部分特别标注：

1. 小域扩张所用的次数
   STARK中定义的FieldElement是以 $3*20^{30}$ 为底数进行运算的，所以小域运算中幂次也需要进行单位化。

You can construct instances of from integers, and then add, multiply, divide, get inverse, and so on. The underlying field of this class is FieldElementFieldElement 𝔽3221225473 ( 3221225473=3*2^30+1), so all operations are done modulo 3221225473.

例如：
小域 $⟨g⟩$ size取 128，g = FieldElement.generator() ** (3 * 2**30//128)
大域 $⟨h⟩$ size取 1024，h = FieldElement.generator() ** (3 * 2**30//1024)

2. 大域扩张在后续的运算中起到什么作用？？
   在更大域上创建RS纠错码，因为作业一主要是在前面算数化，所以暂时不需要用到大域 H（或者说还没用到大域的核心功能）。

3. 数列递推式n的边界条件决定了 $p_2(x)$ 中的 $n$ 在循环中取到多少，斐波那契数列因为递推式有三个前后值所以跟Catalan Series不同。关于range函数在python中的左开右闭特性，可能端点处会出现interpolate_poly函数报错，调试一下即可。

4. 特殊处理，Catalan Series仅在自然数取值时有效。当插值到多项式 $f(x)$ 的时候定义域使用的是乘法子群元素，所以要把 $p_2(x)=\frac{f（g^{x+1}）-\frac{(4x+2)f(g^x)}{x+2}}{\prod_{i=0}^{99} （x-g_i）}$ 中的 $\frac{4x+2}{x+2}$ 特殊处理，将自然数其映射到有限域。（由点集合插值得到多项式 M ）

另：

$p_0(x)=\frac{f(x)-1}{x-g^0}$
$p_1(x)=\frac{f(x)-1555040254}{x-g^{100}}$

Code:

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, X, prod


# generate catalan series
c=[FieldElement(1)]
for i in range(127):
   c.append(c[-1]*(4* FieldElement(i)+2)/(FieldElement(i)+2))

# generate small size=128  domain G
# generate big size=1024  domain H
## interpolate with the small size domain

g = FieldElement.generator() ** (3 * 2**30 // 128)
G = [g ** i for i in range(128)]
f = interpolate_poly(G, c)

w = FieldElement.generator()
h = w ** ((3 * 2**30) // 1024)
H = [h ** i for i in range(1024)]
eval_domain = [w * x for x in H]
f_eval = [f(d) for d in eval_domain]
mt = MerkleTree(f_eval)

test_channel = Channel()
test_channel.send(mt.root)


# px 0+1+2
numer0 = f - 1
denom0 = X - 1
p0 = numer0 / denom0
numer1 = f - 1555040254
denom1 = X - g**100
p1 = numer1 / denom1

# 将自然数其映射到有限域，（由点集合插值得到多项式 M ）
M = interpolate_poly(
        [g**i for i in range(100)],
        [2*(2*FieldElement(i)+1)/(FieldElement(i)+2) for i in range(100)])

p2_numerator = f(g * X) - M(X) * f(X)
p2_denominator = prod([X - g**i for i in range(100)])
p2 = p2_numerator / p2_denominator

def get_CP(channel):
    alpha0 = channel.receive_random_field_element()
    alpha1 = channel.receive_random_field_element()
    alpha2 = channel.receive_random_field_element()
    return alpha0 * p0 + alpha1*p1 + alpha2*p2
    
CP_test = get_CP(test_channel)
print(f'h**111 {CP_test(h**111).val}')
print(f'h**888 {CP_test(h**888).val}')

[Rikki424]

Outcome：
h111 657459872
h888 1514244244

[NeesuNergy]

1.求trace

from field import FieldElement
f = [FieldElement(1)]
for i in range(1,128):
    numer = FieldElement(2) * FieldElement(2 * i - 1) * f[i - 1]
    denom = FieldElement(i + 1)
    f.append(numer / denom)
f[100]

得到f[100] = 1555040254

2.在128domain上插值

from polynomial import interpolate_poly, X
g = FieldElement.generator() ** ((2 ** 30 * 3) // 128)
G = [g ** i for i in range(128)]
c = interpolate_poly(G, f)

3.将domain blowup到1024

w = FieldElement.generator()
h = w ** ((2 ** 30 * 3) // 1024)
H = [h ** i for i in range(1024)]
eval_domain = [w * x for x in H]
f_eval = [c(d) for d in eval_domain]

4.求三个约束的polynomial

numer0 = c - 1
denom0 = X - g ** 0
p0 = numer0 / denom0

numer1 = c - 1555040254
denom1 = X - g ** 100
p1 = numer1 / denom1

N = [FieldElement(i) for i in range(128)]
n = interpolate_poly(G, N)
numer2 = c * (4 * n(X) + 2) - c(g * X) * (n(X) + 2)
denom2 = (X**128 - 1) / (X - g**127)
p2 = numer2 / denom2

5.组合约束求出CP

from channel import Channel
channel = Channel()
alpha0 = channel.receive_random_field_element()
alpha1 = channel.receive_random_field_element()
alpha2 = channel.receive_random_field_element()
CP = alpha0 * p0 + alpha1 * p1 + alpha2 * p2
print(CP(w * h ** 111))

得到h**111点的值为1125323761

[hanbu97]

准备工作

# 参考stark101，引入要用到的库

from channel import Channel
from field import FieldElement
from polynomial import interpolate_poly, X, prod
from merkle import MerkleTree

# 初始化 channel
channel = Channel()

# 找到 100对应的 2的幂
target_n = 100
expand = 8

nearest_pow2 = 1 << (target_n-1).bit_length() if target_n > 0 else 1

small_size = nearest_pow2
large_size = nearest_pow2 * expand
print("Nearest power of 2:", nearest_pow2)
print("small size:", small_size, "large size:", large_size)

Nearest power of 2: 128
small size: 128 large size: 1024

计算Catalan数列

• 根据公式：
  $C_{n+1} = \frac{2(2n+1)}{n+2}C_n$

catalan_series = [FieldElement(1)]
fe1 = FieldElement(1)
fe2 = FieldElement(2)

for i in range(nearest_pow2 - 1):
    fen = FieldElement(i)
    Cn = catalan_series[i]

    Cnadd1 = fe2 * (fe2 * fen + fe1) / (fen + fe2) * Cn
    catalan_series.append(Cnadd1)

print("C_100 =", catalan_series[100])
assert catalan_series[100]==1555040254

C_100 = 1555040254

计算用到的域

小域G: size = 128

大域H: size = 1024

乘一个协因子$w$，原因（待Review）：

1. 打破子群结构：直接使用子群 $H$ 可能导致安全性问题，乘以 $w$ 打破了这种结构。
2. 扩大域范围：将 $H$ "shift"到field的不同部分，扩大了评估域。
3. 增加随机性：$w$ 作为额外的随机因子，增强了整个系统的安全性。

g = FieldElement.generator() ** ((2 ** 30 * 3) // small_size)
print('small g:', g, 'g ^',small_size, '=', g**small_size)

h = FieldElement.generator() ** ((2 ** 30 * 3) // large_size)
print('large h:', h, 'h ^',large_size, '=', h**large_size)

G = [g ** i for i in range(small_size)]
H = [h ** i for i in range(large_size)]

w = FieldElement.generator()
H = [w * hh for hh in H]

small g: 1292405718 g ^ 128 = 1
large h: -1365964089 h ^ 1024 = 1

拉格朗日插值并计算扩展域上的值

f = interpolate_poly(G, catalan_series)
f_eval = [f(hh) for hh in H]

设置随机数种子

f_merkle = MerkleTree(f_eval)
channel.send(f_merkle.root)

边界约束

$p_0(x) = \frac{f(x) - C_0}{x = g^0}$

$p_{1}(x) = \frac{f(x) - C_{100}}{x - g_{100}}$

p0 = (f - catalan_series[0]) / (X - G[0])
p1 = (f - catalan_series[100]) / (X - G[100])

递推关系约束

根据公式：
$C_{n+1} = \frac{2(2n+1)}{n+2}C_n$

即:
$(n+2)C_{n+1} = 2(2n+1)C_n$

即:
$(n+2)C_{n+1} - 2(2n+1)C_n = 0$

因此，需要获得一个n的表达式，满足:

$n_{g^0} = 0$
$n_{g^1} = 1$
$n_{g^2} = 2$
$n_{g^3} = 3$ ...

$p_2(x) = \frac{(n(x)+2)f(gx)-2(2*n(x)+1)*f(x)}{\prod \limits_{i=0}^{100}(x - g^{i})}$

另一种方式，扩展到全域，分母可以简化:

$p_2(x) = \frac{(n(x)+2)f(gx)-2(2*n(x)+1)*f(x)}{（x^{128}-1)/(x-g^{127})}$

# part
n = interpolate_poly(G[:target_n], [FieldElement(i) for i in range(target_n)])

p2_numerator = (n(X) + fe2) * f(g*X) - fe2 * (fe2 * n(X) + fe1) * f(X)  
p2_denominator = prod([X - g**i for i in range(target_n)])
p2_part = p2_numerator / p2_denominator


# full
n = interpolate_poly(G, [FieldElement(i) for i in range(small_size)])

p2_numerator = (n(X) + fe2) * f(g*X) - fe2 * (fe2 * n(X) + fe1) * f(X)   
# p2_numerator = fe2 * (fe2 * n(X) + fe1) * f(X) - (n(X) + fe2) * f(g*X) 
p2_denominator = (X**128 - 1) / (X-g**127)
p2_full = p2_numerator / p2_denominator

print("p0 degree:", p0.degree())
print("p1 degree:", p1.degree())
print("p2 part degree:", p2_part.degree())
print("p2 full degree:", p2_full.degree())

p0 degree: 126
p1 degree: 126
p2 part degree: 126
p2 full degree: 127

构造CP

alpha0 = channel.receive_random_field_element()
alpha1 = channel.receive_random_field_element()
alpha2 = channel.receive_random_field_element()

print('alpha0 =', alpha0)
print('alpha1 =', alpha1)
print('alpha2 =', alpha2)

CP_part = p0 + alpha1 * p1 + alpha2 * p2_part
CP_full = p0 + alpha1 * p1 + alpha2 * p2_full

print('CP_part degree:', CP_part.degree())
print('CP_full degree:', CP_full.degree())

alpha0 = -873511777
alpha1 = -178162817
alpha2 = -442740964
CP_part degree: 126
CP_full degree: 127

计算结果

print('part h111:', CP_part(h**111).val)
print('part h888:', CP_part(h**888).val)
print('full h111:', CP_full(h**111).val)
print('full h888:', CP_full(h**888).val)

part h111: 2664842244
part h888: 2759235016
full h111: 498748567
full h888: 1923900075

[tmors]

from field import FieldElement
a = [FieldElement(1)]
for i in range(1, 100):
    a.append(a[-1] * 2 * (2 * i + 1) / (i + 2))
assert len(a) == 100, 'The trace must consist of exactly 1023 elements.'
assert a[0] == FieldElement(1), 'The first element in the trace must be the unit element.'
for i in range(1, 100):
    assert a[i] == a[i-1] * 2 * (2 * i + 1) / (i + 2), f'The FibonacciSq recursion rule does not apply for index {i}'
assert a[99] == FieldElement(1555040254), 'Wrong last element!'
print('Success!')
g = FieldElement.generator() ** (3 * 2 ** 23)
G = [g ** i for i in range(128)]
assert g.is_order(128), 'The generator g is of wrong order.'
b = FieldElement(1)
for i in range(127):
  assert b == G[i], 'The i-th place in G is not equal to the i-th power of g.'
  b = b * g
  assert b != FieldElement(1), f'g is of order {i + 1}'
if b * g == FieldElement(1):
  print('Success!')
else:
  print('g is of order > 1024')
from polynomial import interpolate_poly
f = interpolate_poly(G[:-28], a)
w = FieldElement.generator()
h = w ** ((2 ** 30 * 3) // 1024)
H = [h ** i for i in range(1024)]
eval_domain = [w * x for x in H]
from hashlib import sha256
assert len(set(eval_domain)) == len(eval_domain)
w = FieldElement.generator()

f_eval = [f(d) for d in eval_domain]

from polynomial import interpolate_poly, X, prod
numer0 = f - FieldElement(1)
denom0 = X - g ** 0
p0,r0  = numer0.qdiv(denom0) // constraint 0

numer1 = f - FieldElement(1555040254)
denom1 = X - g**99
p1, r1 = numer1.qdiv(denom1) // constraint 1

lst = [(X - g**i) for i in range(99)]
pr = prod(lst) * (X + 2)

M_points = [(g**i, FieldElement(i)) for i in range(99)]
M = interpolate_poly([p[0] for p in M_points], [p[1] for p in M_points])

def catalan_constraint(x):
  return f(g * x) * (M + 3) - f(x) * (FieldElement(2)*(FieldElement(2)*M+3))

p2_numerator = catalan_constraint(X) // constrain 2
p2_denominator = prod([X - g**i for i in range(99)])
p2, r22 = p2_numerator.qdiv(p2_denominator)

def get_CP(channel):
  alpha0 = channel.receive_random_field_element()
  alpha1 = channel.receive_random_field_element()
  alpha2 = channel.receive_random_field_element()
  return alpha0 * p0 + alpha1 * p1 + alpha2 * p2

def CP_eval(channel):
  CP = get_CP(channel)
  return [CP(d) for d in eval_domain]

from channel import Channel
test_channel = Channel()
CP_test = get_CP(test_channel)


print(CP_test(111)) //963694142
print(CP_test(888)) //-1381902967

[YoulongDing]

有手就行，注意channel一定要先send merkel root，再用它的randomness来构造CP

from field import FieldElement
from channel import Channel
from polynomial import interpolate_poly, X, prod
from merkle import MerkleTree
from channel import Channel


a = [FieldElement(1), FieldElement(3141592)]
while len(a) < 1023:
    a.append(a[-2] * a[-2] + a[-1] * a[-1])

a = [FieldElement(1)]

for i in range(100):
    a.append(FieldElement(2) * (FieldElement(2) * FieldElement(i) + 1) * FieldElement(i+2).inverse() * a[i])

assert a[-1] == 1555040254

g = FieldElement.generator() ** ((3*2**30) // 128)
G = [g**i for i in range(128)]
f = interpolate_poly(G[:len(a)], a)
assert f(g**100) == 1555040254

w = FieldElement.generator()
h = w ** ((2 ** 30 * 3) // 1024)
H = [h ** i for i in range(1024)]
eval_domain = [w * h for h in H]
f_eval = [f(d) for d in eval_domain]

f_merkle = MerkleTree(f_eval)
channel = Channel()
channel.send(f_merkle.root)

assert (f - 1) % (X - G[0]) == 0
assert (f - 1555040254) % (X - G[100]) == 0

p0 = (f - 1) / (X - G[0])
p1 = (f - 1555040254) / (X - G[100])

l = interpolate_poly([g**i for i in range(100)], [FieldElement(i) for i in range(100)])

numerator = f(g * X) * (l + 2) - (FieldElement(2) * (FieldElement(2) * l + 1) * f(X))
denominator = prod([X - g**i for i in range(100)])

assert numerator % denominator == 0

p2 = numerator / denominator

def get_CP(channel):
    alpha0 = channel.receive_random_field_element()
    alpha1 = channel.receive_random_field_element()
    alpha2 = channel.receive_random_field_element()
    return alpha0*p0 + alpha1*p1 + alpha2*p2

cp = get_CP(channel)
cp(h ** 111)

732987908

[yiyanwannian]

from channel import Channel
from field import FieldElement
from merkle import MerkleTree
from polynomial import interpolate_poly, X, prod

// 计算卡特兰数列
def calculate_catalan_sequence(size, val_100):
t = [FieldElement(1)]
for i in range(size - 1):
t.append(t[-1] * 2 * (FieldElement(2) * i + 1) / (i + 2))
assert t[100] == val_100
return t
c=[1]
t=[FieldElement(1)]
for i in range(127):
c.append(c[-1]2(2*i+1)//(i+2))
t.append(FieldElement(c[-1]))
assert t[100]==1555040254

// 准备域
def prepare_domain(size):
generator = FieldElement.generator()
g = generator ** (3 * 2 ** 30 // size)
return [g ** i for i in range(size)]

// 插值并求值
def interpolate_and_evaluate(G, t, domain):
f = interpolate_poly(G, t)
return f, [f.eval(d) for d in domain]

// 构造约束条件
def make_constraints(f, g, g0, g100, val_100):

// 构造第一个约束多项式
p0 = (f - 1) / (X - g0)

// 构造第二个约束多项式
p1 = (f - val_100) / (X - g100)

// 构造辅助多项式fn
N = [FieldElement(i) for i in range(100)]
fn = interpolate_poly(G[:100], N)
assert fn(g ** 99) == 99

// 构造第三个约束多项式
numer = (fn(X) + 2) * f(g * X) - (4 * fn(X) + 2) * f
alpha = prod([(X - g ** i) for i in range(100)])
p2 = numer / alpha

return p0, p1, p2, fn

// 获取组合多项式
def get_composition_polynomial(p0, p1, p2, ch):
ch.send(mt.root)

// 接收随机系数
alpha0 = ch.receive_random_field_element()
alpha1 = ch.receive_random_field_element()
alpha2 = ch.receive_random_field_element()

// 构造组合多项式
return alpha0 * p0 + alpha1 * p1 + alpha2 * p2

if name == "main":
small_domain_size = 128
large_domain_size = 1024
expect_value_100 = 1555040254

// 计算卡特兰数列
t = calculate_catalan_sequence(small_domain_size, expect_value_100)

// 准备小域G和大域H
G = prepare_domain(small_domain_size)
H = prepare_domain(large_domain_size)
domain = [FieldElement.generator() * x for x in H]

f, ev = interpolate_and_evaluate(G, t, domain)
assert f(G[100]) == expect_value_100

// 构造约束条件
p0, p1, p2, fn = make_constraints(f, G[1], G[0], G[100], expect_value_100)

// 打印约束多项式的度
print(f"deg p0 = {p0.degree()}")
print(f"deg p1 = {p1.degree()}")
print(f"deg p2 = {p2.degree()}")

// 构建Merkle树
mt = MerkleTree(ev)
ch = Channel()

// 获取组合多项式
cp = get_composition_polynomial(p0, p1, p2, ch)

// 打印结果
assert cp(G[111]).val == cp(H[888]).val
print(f'cp(g111): {cp(G[111]).val}')
print(f'cp(h888): {cp(H[888]).val}')
print(f'cp(h**111): {cp(H[111]).val}')

[wjlkk1]

太难了，之前没有这方面的基础，我尽量在做，也参考了其他人的作业和课件，可能有错误，没必要太参考我的答案
大饼10w刀以后我的想法又清晰了一些，未来是属于区块链、zkp这些的没错，先打好基础很重要，各位别嫌我菜不带我学🥺
已知
$C_0=1$, $C_{n+1}=\frac {2(2n + 1)}{n+2}C_n$
要给verifier证明
$C_{100}=1555040254$
只需证明约束(分式)的组合是一个多项式即可

思路

导入sagemath库

from sage.all import *

copy stark101有限域相同设置并计算n=127内catalan数

# 定义有限域的模数
k_modulus = 3 * 2**30 + 1
# 创建有限域 F_(3 * 2**30 + 1)
F = GF(k_modulus, name='a')

def Catalan(pre, n):
    return F(2)*(F(2)*n+1)/(n+2)*pre

# 定义一个生成元
generator_val = 5
generator = F(generator_val)

C_0 = F(1)
elements = [C_0]
for i in range(127):
    elements.append(Catalan(elements[-1], i))
assert elements[100] == 1555040254

使用lagrange插值得到插值多项式

# 创建多项式环
R = PolynomialRing(F, 'x')
X = R.gen()

def interpolate_poly(x_values, y_values):
    """
    Returns a polynomial of degree < len(x_values) that evaluates to y_values[i] on x_values[i] for
    all i.
    """
    assert len(x_values) == len(y_values)
    points = list(zip(x_values, y_values))
    poly = R.lagrange_polynomial(points)
    return poly

g = generator **((k_modulus - 1)//128)
points = [g**i for i in range(128)]

# 得到插值多项式
f = interpolate_poly(points, elements)

# Catalan约束所需多项式
i = interpolate_poly(points, [F(i) for i in range(128)])

定义并打印约束，约束是分式形式，需要证明他们都是多项式

# 证明满足多项式的3个约束
# constraint 1
p0 = (f - 1)/ (X - g**0)
# constraint 2
p1 = (f - 1555040254)/ (X - g**100)

prod = 1
for j in range(128):
    prod *= (X - g**j)
assert (prod == (X**128 - 1)) 
# constraint 3
nu = (f(g*X)*(i(X) + 2) - f(X)*(2*(2*i(X) + 1)))
de = (X**128 - 1)/((X - g**127)*(X - g**126))
p2 = nu/de

print(p0)
print(p1)
print(p2)

"""
output:
2241441300*x^126 + 969514308*x^125 + 286576041*x^124 + ...
2241441300*x^126 + 3032696846*x^125 + 3043428958*x^124 + ...
2151462736*x^128 + 1234103371*x^127 + 3011284354*x^126 + ...
"""

转换为证明CP是一个多项式，使用有限域中随机值初始化三个约束，组合起来得到Composition polynomial并在扩展域上验证取值

alpha0 = R.random_element()
alpha1 = R.random_element()
alpha2 = R.random_element()
CP = alpha0*p0 + alpha1*p1 + alpha2*p2

h = generator**((k_modulus - 1)//1024)

print(CP(h**111))
print(CP(h**888))
print(CP(g**111))
assert (CP(g**77) == CP(h**(77*8)))

"""
output:
1304884559
1758307133
1758307133
"""