Anatomy of a STARK, Part 2: Basic Tools（译）

source: https://aszepieniec.github.io/stark-anatomy/basic-tools

目录

• Finite Field Arithmetic
• Univariate Polynomials(单变量多项式)
• Multivariate Polynomials（多元多项式）
• The Fiat-Shamir Transform
• merkle tree

Finite Field Arithmetic

Finite fields are ubiquitous throughout cryptography because they are natively compatible with computers. For instance, they cannot generate overflow or underflow errors, and their elements have a finite representation in terms of bits.

The easiest way to build a finite field is to select a prime number \(p\), use the elements \({\rm F}_{p} = \{0,1...,p\}\) , and define the usual addition and multiplication operations in terms of their counterparts for the integers, followed by reduction modulo \(p\) .
Subtraction is equivalent to addition of the left hand side to the negation of the right hand side, and negation represents multiplication by \(-1 \equiv p-1\) mod \(p\) .
Similarly, division is equivalent to multiplication of the left hand side by the multiplicative inverse of the right hand side. This inverse can be found using the extended Euclidean algorithm, which on input two integers \(x\) and \(y\), returns their greatest common divisor \(g\) along with matching Bezout coefficients \(a\), \(b\) such that \(ax+by=g\) . Indeed, whenever \(gcd(x,p)=1\) the inverse of \(x \in {\rm F}_p\) is \(a\) because \(ax+bp \equiv 1\) mod \(p\).
Powers of field elements can be computed with the square-and-multiply algorithm, which iterates over the bits in the expansion of the exponent, squares an accumulator variable in each iteration, and additionally multiplies it by the base element if the bit is set.

为了构建 STARK，我们需要具有特定结构的有限域: 对于一些足够大的 \(k\)，它需要包含一个 \(2^k\) 阶的子结构。我们考虑模 \(p\) 的定义为如下格式的素数域 \(p = f * 2^k + 1\), 其中 \(f\) 是使数成为素数的辅因子。In this case, the group \(F_p\) has a subgroup of order \(2^k\). 我们可以在复数单位圆上用 \(2^k\) 个均匀分布的点来识别这个子群。

An implementation starts with the extended Euclidean algorithm, for computing multiplicative inverses.

# 扩展欧几里德算法 the extended Euclidean algorithm
# a, _, _ = xgcd(x, p), a 为 x 在模为 p 的有限域上的逆
def xgcd( x, y ):
    old_r, r = (x, y)
    old_s, s = (1, 0)
    old_t, t = (0, 1)

    while r != 0:
        quotient = old_r // r
        old_r, r = (r, old_r - quotient * r)
        old_s, s = (s, old_s - quotient * s)
        old_t, t = (t, old_t - quotient * t)

    return old_s, old_t, old_r # a, b, g

It makes sense to separate the logic concerning the field from the logic concerning the field elements. To this end, the field element contains a field object as a proper field; this field object implements the arithmetic. Furthermore, python supports operator overloading, so we can repurpose natural arithmetic operators to do field arithmetic instead.

# 有限域的元素定义
class FieldElement:
    def __init__( self, value, field ):
        self.value = value
        self.field = field

    def __add__( self, right ):
        return self.field.add(self, right)

    def __mul__( self, right ):
        return self.field.multiply(self, right)

    def __sub__( self, right ):
        return self.field.subtract(self, right)

    def __truediv__( self, right ):
        return self.field.divide(self, right)

    def __neg__( self ):
        return self.field.negate(self)

    def inverse( self ):
        return self.field.inverse(self)

    # modular exponentiation -- be sure to encapsulate in parentheses!
    def __xor__( self, exponent ):
        acc = FieldElement(1, self.field)
        val = FieldElement(self.value, self.field)
        for i in reversed(range(len(bin(exponent)[2:]))):
            acc = acc * acc
            if (1 << i) & exponent != 0:
                acc = acc * val
        return acc

    def __eq__( self, other ):
        return self.value == other.value

    def __neq__( self, other ):
        return self.value != other.value

    def __str__( self ):
        return str(self.value)

    def __bytes__( self ):
        return bytes(str(self).encode())

    def is_zero( self ):
        if self.value == 0:
            return True
        else:
            return False

# 有限域的定义, 有限域上的计算需要 mod p
class Field:
    def __init__( self, p ):
        self.p = p

    def zero( self ):
        return FieldElement(0, self)

    def one( self ):
        return FieldElement(1, self)

    def multiply( self, left, right ):
        return FieldElement((left.value * right.value) % self.p, self)

    def add( self, left, right ):
        return FieldElement((left.value + right.value) % self.p, self)

    def subtract( self, left, right ):
        return FieldElement((self.p + left.value - right.value) % self.p, self)

    def negate( self, operand ):
        return FieldElement((self.p - operand.value) % self.p, self)

    def inverse( self, operand ):
        a, b, g = xgcd(operand.value, self.p)
        return FieldElement(a, self)

    def divide( self, left, right ):
        assert(not right.is_zero()), "divide by zero"
        a, b, g = xgcd(right.value, self.p)
        return FieldElement(left.value * a % self.p, self)

Implementing fields generically is nice. However, in this tutorial we will not use any other field than the one with \(1+407⋅2^{119}\) elements. This field has a sufficiently large subgroup of power-of-two order.

def main():
    # 教程中 p 是写死的
    p = 1 + 407 * ( 1 << 119 ) # 1 + 11 * 37 * 2^119
    return Field(p)

Besides ensuring that the subgroup of power-of-two order exists, the code also needs to supply the user with a generator for the entire multiplicative group, as well as the power-of-two subgroups. A generator for such a subgroup of order \(n\) will be called a primitive \(n\)th root.

def generator( self ):
    assert(self.p == 1 + 407 * ( 1 << 119 )), "Do not know generator for other fields beyond 1+407*2^119"
    return FieldElement(85408008396924667383611388730472331217, self)

# 这里的root满足 root^n == 1    
def primitive_nth_root( self, n ):
    if self.p == 1 + 407 * ( 1 << 119 ):
        assert(n <= 1 << 119 and (n & (n-1)) == 0), "Field does not have nth root of unity where n > 2^119 or not power of two."
        root = FieldElement(85408008396924667383611388730472331217, self)
        order = 1 << 119
        while order != n:
            root = root^2
            order = order/2
        return root
    else:
        assert(False), "Unknown field, can't return root of unity."

Lastly, the protocol requires the ability to sample field elements randomly and pseudorandomly. To do this, the user supplies random bytes and the field logic turns them into a field element. The user should take care to provide enough random bytes.

# 随机数生成：输入随机的byte数组，输出有限域内的对应的数
def sample( self, byte_array ):
    acc = 0
    for b in byte_array:
        acc = (acc << 8) ^ int(b)
    return FieldElement(acc % self.p, self)

Univariate Polynomials(单变量多项式)

形如 \(f(X) = c0 + c1⋅X + ... + c_{d}X^{d}\) 或者 \(f(X) = \sum_{i=0}^d c_{i}X^{i}\) , 其中 \(c_i\) 为 系数, \(d\) 为多项式的 阶。

Univariate polynomials are immensely useful in proof systems because relations that apply to their coefficient vectors extend to their values on a potentially much larger domain. 如果多项式相等，则它们处处相等； 而如果它们不相等，那么它们几乎在所有地方都是不相等的。通过这个特性，单变量多项式将关于大向量的声明减少为在一小部分足够随机的点中声明其对应多项式的值。

from algebra import *

class Polynomial:
    def __init__( self, coefficients):
        # 由高阶系数到低阶系数, coefficients[0]是最高阶，不可能为0
        self.coefficients = [c for c in coefficients]

    def degree( self ):
        if self.coefficients == []:
            return -1
        zero = self.coefficients[0].field.zero()
        if self.coefficients == [zero] * len(self.coefficients):
            return -1
        maxindex = 0
        for i in range(len(self.coefficients)):
            if self.coefficients[i] != zero:
                maxindex = i
        return maxindex

    def __neg__( self ):
        return Polynomial([-c for c in self.coefficients])

    def __add__( self, other ):
        if self.degree() == -1:
            return other
        elif other.degree() == -1:
            return self
        field = self.coefficients[0].field
        coeffs = [field.zero()] * max(len(self.coefficients), len(other.coefficients))
        for i in range(len(self.coefficients)):
            coeffs[i] = coeffs[i] + self.coefficients[i]
        for i in range(len(other.coefficients)):
            coeffs[i] = coeffs[i] + other.coefficients[i]
        return Polynomial(coeffs)

    def __sub__( self, other ):
        return self.__add__(-other)

    # 多项式相乘
    def __mul__(self, other ):
        if self.coefficients == [] or other.coefficients == []:
            return Polynomial([])
        zero = self.coefficients[0].field.zero()
        buf = [zero] * (len(self.coefficients) + len(other.coefficients) - 1)
        for i in range(len(self.coefficients)):
            if self.coefficients[i].is_zero():
                continue # optimization for sparse polynomials
            for j in range(len(other.coefficients)):
                buf[i+j] = buf[i+j] + self.coefficients[i] * other.coefficients[j]
        return Polynomial(buf)

    def __eq__( self, other ):
        if self.degree() != other.degree():
            return False
        if self.degree() == -1:
            return True
        return all(self.coefficients[i] == other.coefficients[i] for i in range(len(self.coefficients)))

    def __neq__( self, other ):
        return not self.__eq__(other)

    def is_zero( self ):
        if self.degree() == -1:
            return True
        return False

    def leading_coefficient( self ):
        return self.coefficients[self.degree()]

This always get a little tricky when implementing division of polynomials. The intuition behind the schoolbook algorithm is that in every iteration you multiply the dividend by the correct term so as to generate a cancellation of leading terms. Once no such term exists, you have your remainder.

# divide(分子, 分母)
def divide( numerator, denominator ):
    if denominator.degree() == -1:
        return None
    if numerator.degree() < denominator.degree():
        return (Polynomial([]), numerator)
    field = denominator.coefficients[0].field
    remainder = Polynomial([n for n in numerator.coefficients])
    quotient_coefficients = [field.zero() for i in range(numerator.degree()-denominator.degree()+1)]
    for i in range(numerator.degree()-denominator.degree()+1):
        if remainder.degree() < denominator.degree():
            break
        coefficient = remainder.leading_coefficient() / denominator.leading_coefficient()
        shift = remainder.degree() - denominator.degree()
        subtractee = Polynomial([field.zero()] * shift + [coefficient]) * denominator
        quotient_coefficients[shift] = coefficient
        remainder = remainder - subtractee
    quotient = Polynomial(quotient_coefficients)
    return quotient, remainder

def __truediv__( self, other ):
    quo, rem = Polynomial.divide(self, other)
    assert(rem.is_zero()), "cannot perform polynomial division because remainder is not zero"
    return quo

def __mod__( self, other ):
    quo, rem = Polynomial.divide(self, other)
    return rem

In terms of basic arithmetic operations, it is worth including a powering map, although mostly for notational easy rather than performance.

def __xor__( self, exponent ):
    if self.is_zero():
        return Polynomial([])
    if exponent == 0:
        return Polynomial([self.coefficients[0].field.one()])
    acc = Polynomial([self.coefficients[0].field.one()])
    for i in reversed(range(len(bin(exponent)[2:]))):
        acc = acc * acc
        if (1 << i) & exponent != 0:
            acc = acc * self
    return acc

如果多项式不允许在给定的任意点计算其值，则它是毫无意义的。 对于 STARK，我们需要一些更通用的东西——对一个值域进行多项式评估，而不是在一个点上。 此时我们不考虑性能的问题，因此以下实现遵循简单的迭代方法。此外，STARK 也需要多项式插值，其中 x 坐标是另一个已知的值范围。重申一次，我们暂时不考虑性能的问题，所以目前标准的拉格朗日插值就足够了。

# 计算单个点处的取值
def evaluate( self, point ):
    xi = point.field.one()
    value = point.field.zero()
    for c in self.coefficients:
        value = value + c * xi
        xi = xi * point
    return value

# 计算一个domain上所有点的取值
def evaluate_domain( self, domain ):
    return [self.evaluate(d) for d in domain]

# 根据domain上所有点的取值求多项式的系数
def interpolate_domain( domain, values ):
    assert(len(domain) == len(values)), "number of elements in domain does not match number of values -- cannot interpolate"
    assert(len(domain) > 0), "cannot interpolate between zero points"
    field = domain[0].field
    x = Polynomial([field.zero(), field.one()])
    acc = Polynomial([])
    for i in range(len(domain)):
        prod = Polynomial([values[i]])
        for j in range(len(domain)):
            if j == i:
                continue
            prod = prod * (x - Polynomial([domain[j]])) * Polynomial([(domain[i] - domain[j]).inverse()])
        acc = acc + prod
    return acc

说到域：一次又一次重复的事是计算它们上的消失多项式。任何这样的多项式都是 \(Z_D = \prod_{d \in D}^{}(X-d)\) 的倍数，\(Z_D\) 是唯一的 a) 最高阶系数为1; b)在 \(D\) 的所有点处取值为0的; c) 阶数最低的多项式. 该多项式通常称为消失多项式(vanishing polynomial)，有时也称为归零多项式(zerofier)。本教程更喜欢第二个术语。

# 归零多项式 Z
def zerofier_domain( domain ):
    field = domain[0].field
    x = Polynomial([field.zero(), field.one()])
    acc = Polynomial([field.one()])
    for d in domain:
        acc = acc * (x - Polynomial([d]))
    return acc

Another useful tool is the ability to scale polynomials. 具体来说，这意味着从 \(f(X)\) 的系数中获得 \(f(c·X)\) 的系数向量. This function is particularly useful when \(f(X)\) is defined to take a sequence of values on the powers of \(c:v_i = f(c^i)\). Then \(f(c·X)\) represents the same sequence of values but shifted by one position.

# [c3, c2, c1, c0] 变为 [factor^3 * c3, factor^2 * c2, factor^1 * c1, c0]
def scale( self, factor ):
    return Polynomial([(factor^i) * self.coefficients[i] for i in range(len(self.coefficients))])

The last function that belongs to the univariate polynomial module anticipates a key operations in the FRI protocol, namely testing whether a triple of points fall on the same line – a fancy word for which is colinearity.

# 检查点是否在同一条直线上
def test_colinearity( points ):
    domain = [p[0] for p in points]
    values = [p[1] for p in points]
    polynomial = Polynomial.interpolate_domain(domain, values)
    return polynomial.degree() == 1 # 阶为1说明点在同一条直线上

在继续下一节之前，值得注意的是所有成分都适用于有限扩展域(finite extension fields)，或者简单的称之为扩展域(extension fields)。 有限域是一个简单的集合，它配备了符合高中代数规则的加法和乘法运算符，例如 每个非零元素都有一个逆元素，或者没有两个非零元素相乘得到零。有两种获取方式：

1. 从整数集开始, 以给定素数 \(p\) 为模减少任何加法或乘法的结果
2. 从有限域上的多项式集合开始，并以给定的不可约多项式(irreducible polynomial) \(p(X)\) 为模，减少任何加法或乘法的结果。当多项式不能分解为两个较小多项式的乘积时，多项式是不可约的，类似于素数。

关键是可以在比密码编译步骤更小的域中进行算术化，只要后一步使用前一步的扩展域。例如，EthSTARK 在由 62 位素数定义的有限域上运行，但为了达到更高的安全目标，FRI的步骤在其二次扩展域(quadratic extension field)上运行。

本教程不会使用扩展域，因此对该主题的详细讨论超出了本文的范围。

Multivariate Polynomials（多元多项式）

多元多项式将单变量多项式推广到许多未知数——不仅是 \(X\)，还有 \(X\), \(Y\), \(Z\) ...。单变量多项式可用于将关于大向量的大声明减少到关于随机点中标量值的小声明，而多元多项式可用于表示整个计算满足的算术约束。

在实现单变量多项式的结构是系数列表的情况下，多元多项式的结构是将指数向量映射到系数的字典([exponents] => coefficient)。系数为0的项会被忽略掉。像往常一样，第一步是重载标准算术运算符、基本构造函数和标准功能。

# dict的key为指数列表，value为系数
# {[2,3,4]:2, [4,5,6]:5} 表示多项式 f(X,Y,Z) = 2*X^2*Y^3*Z^4 + 5*X^4*Y^5*Z^6
class MPolynomial:
    def __init__( self, dictionary ):
        self.dictionary = dictionary

    def zero():
        return MPolynomial(dict())

    def __add__( self, other ):
        dictionary = dict()
        num_variables = max([len(k) for k in self.dictionary.keys()] + [len(k) for k in other.dictionary.keys()])
        for k, v in self.dictionary.items():
            pad = list(k) + [0] * (num_variables - len(k))
            pad = tuple(pad)
            dictionary[pad] = v
        for k, v in other.dictionary.items():
            pad = list(k) + [0] * (num_variables - len(k))
            pad = tuple(pad)
            if pad in dictionary.keys():
                dictionary[pad] = dictionary[pad] + v
            else:
                dictionary[pad] = v
        return MPolynomial(dictionary)

    def __mul__( self, other ):
        dictionary = dict()
        num_variables = max([len(k) for k in self.dictionary.keys()] + [len(k) for k in other.dictionary.keys()])
        for k0, v0 in self.dictionary.items():
            for k1, v1 in other.dictionary.items():
                exponent = [0] * num_variables
                for k in range(len(k0)):
                    exponent[k] += k0[k]
                for k in range(len(k1)):
                    exponent[k] += k1[k]
                exponent = tuple(exponent)
                if exponent in dictionary.keys():
                    dictionary[exponent] = dictionary[exponent] + v0 * v1
                else:
                    dictionary[exponent] = v0 * v1
        return MPolynomial(dictionary)

    def __sub__( self, other ):
        return self + (-other)

    def __neg__( self ):
        dictionary = dict()
        for k, v in self.dictionary.items():
            dictionary[k] = -v
        return MPolynomial(dictionary)

    def __xor__( self, exponent ):
        if self.is_zero():
            return MPolynomial(dict())
        field = list(self.dictionary.values())[0].field
        num_variables = len(list(self.dictionary.keys())[0])
        exp = [0] * num_variables
        acc = MPolynomial({tuple(exp): field.one()})
        for b in bin(exponent)[2:]:
            acc = acc * acc
            if b == '1':
                acc = acc * self
        return acc

    def constant( element ):
        return MPolynomial({tuple([0]): element})

    def is_zero( self ):
        if not self.dictionary:
            return True
        else:
            for v in self.dictionary.values():
                if v.is_zero() == False:
                    return False
            return True

    def variables( num_variables, field ):
        variables = []
        for i in range(num_variables):
            exponent = [0] * i + [1] + [0] * (num_variables - i - 1)
            variables = variables + [MPolynomial({tuple(exponent): field.one()})]
        return variables

Since multivariate polynomials are a generalization of univariate polynomials, there needs to be a way to inherit the logic that was already defined for the former class. The function lift does this by lifting a univariate polynomial to the multivariate polynomials. The second argument is the index of the indeterminate that corresponds to the univariate indeterminate.

# 将单变量多项式转化为多元多项式
def lift( polynomial, variable_index ):
    if polynomial.is_zero():
        return MPolynomial({})
    field = polynomial.coefficients[0].field
    variables = MPolynomial.variables(variable_index+1, field)
    x = variables[-1]
    acc = MPolynomial({})
    for i in range(len(polynomial.coefficients)):
        acc = acc + MPolynomial.constant(polynomial.coefficients[i]) * (x^i)
    return acc

接下来是求值(evaluation)。 此方法的参数需要是一个标量元组，因为它需要为每个变量分配一个值。 然而，值得注意 STARK 中使用的一个特征，即评估是 象征性(symbolic) 的：它不是在一个标量元组中评估多元多项式，而是在一个单变量多项式元组中评估。 结果不是一个标量，而是一个新的单变量多项式。

# point是一个元组，因为需要同时提供X,Y,Z...等处的值
# 结果为f(X,Y,Z...)在这里的取值，是一个标量
def evaluate( self, point ):
    acc = point[0].field.zero()
    for k, v in self.dictionary.items():
        prod = v
        for i in range(len(k)):
            prod = prod * (point[i]^k[i])
        acc = acc + prod
    return acc 

# dict中的每一项变成了多项式的系数
# 最后输出的是一个单变量多项式
def evaluate_symbolic( self, point ):
    acc = Polynomial([])
    for k, v in self.dictionary.items():
        prod = Polynomial([v])
        for i in range(len(k)):
            prod = prod * (point[i]^k[i])
        acc = acc + prod
    return acc

The Fiat-Shamir Transform

In an interactive public coin protocol, the verifier’s messages are pure randomness sampled from a distribution that anyone can sample from. 目标是在不牺牲安全性的情况下获得证明效果相同的非交互式协议。Fiat-Shamir 变换实现了这一点。

It turns out that for generating security against malicious provers, generating the verifier’s messages randomly as the interactive protocol stipulates, is overkill. It is sufficient that the verifier’s messages be difficult to predict by the prover. Hash functions are deterministic but still satisfy this property of outputs being difficult to predict. 所以直观地说，如果验证者的真实随机性被散列函数的伪随机输出替换，协议仍然是安全的. It is necessary to restrict the prover’s control over what input goes into the hash function, because otherwise he can grind until he finds a suitable output. It suffices to set the input to the transcript of the protocol up until the point where the verifier’s message is needed.

This is exactly the intuition behind the Fiat-Shamir transform: replace the verifier’s random messages by the hash of the transcript of the protocol up until those points. The Fiat-Shamir heuristic states that this transform retains security. In an idealized model of the hash function called the random oracle model, this security is provable.

The Fiat-Shamir transform presents the first engineering challenge. The interactive protocol is described in terms of a channel which passes messages from prover to receiver or the other way around. The transform serializes this communication while enabling a description of the prover that makes abstraction of it. The transform does modify the description of the verifier, which becomes deterministic.

A proof stream is a useful concept to simulate this channel. The difference with respect to regular streams in programming is that there is no actual transmission to another process or computer taking place, and nor do sender and receiver need to operate simultaneously. It is not a simple queue either because the prover and the verifier have access to a function that computes pseudorandomness by hashing their view of the channel. 对于证明者，此视图是迄今为止发送的所有消息的完整列表。 对于验证者来说，这个视图是目前已读消息的子列表。验证者的消息不会添加到列表中，因为它们可以被确定性的计算出来。给定证明者的消息列表，序列化很简单。非交互式证明正是这种序列化。(可以看到Proof Stream替代了原有的交互流程中verifier的角色)

In terms of implementation, what is needed is a class ProofStream that supports 3 functionalities.

1. Pushing and pulling objects to and from a queue. The queue is simulated by a list with a read index. Whenever an item is pushed, it is appended. Whenever an item is pulled, the read index is incremented by one.
2. Serialization and deserialization. The amazing python library pickle does this.
3. Fiat-Shamir. Hashing is done below by first serializing the queue or the first part of it, and then applying SHAKE-256. SHAKE-256 admits a variable output length, which the particular application may want to set. By default the output length is set to 32 bytes.

from hashlib import shake_256
import pickle as pickle # serialization

class ProofStream:
    def __init__( self ):
        self.objects = []
        self.read_index = 0

    def push( self, obj ):
        self.objects += [obj]

    def pull( self ):
        assert(self.read_index < len(self.objects)), "ProofStream: cannot pull object; queue empty."
        obj = self.objects[self.read_index]
        self.read_index += 1
        return obj

    def serialize( self ):
        return pickle.dumps(self.objects)

    def deserialize( bb ):
        ps = ProofStream()
        ps.objects = pickle.loads(bb)
        return ps

    def prover_fiat_shamir( self, num_bytes=32 ):
        return shake_256(self.serialize()).digest(num_bytes)

    def verifier_fiat_shamir( self, num_bytes=32 ):
        return shake_256(pickle.dumps(self.objects[:self.read_index])).digest(num_bytes)

merkle tree

from hashlib import blake2b

class Merkle:
    H = blake2b

    def commit_( leafs ):
        assert(len(leafs) & (len(leafs)-1) == 0), "length must be power of two"
        if len(leafs) == 1:
            return leafs[0]
        else:
            return Merkle.H(Merkle.commit_(leafs[:len(leafs)//2]) + Merkle.commit_(leafs[len(leafs)//2:])).digest()
    
    def open_( index, leafs ):
        assert(len(leafs) & (len(leafs)-1) == 0), "length must be power of two"
        assert(0 <= index and index < len(leafs)), "cannot open invalid index"
        if len(leafs) == 2:
            return [leafs[1 - index]]
        elif index < (len(leafs)/2):
            return Merkle.open_(index, leafs[:len(leafs)//2]) + [Merkle.commit_(leafs[len(leafs)//2:])]
        else:
            return Merkle.open_(index - len(leafs)//2, leafs[len(leafs)//2:]) + [Merkle.commit_(leafs[:len(leafs)//2])]
    
    def verify_( root, index, path, leaf ):
        assert(0 <= index and index < (1 << len(path))), "cannot verify invalid index"
        if len(path) == 1:
            if index == 0:
                return root == Merkle.H(leaf + path[0]).digest()
            else:
                return root == Merkle.H(path[0] + leaf).digest()
        else:
            if index % 2 == 0:
                return Merkle.verify_(root, index >> 1, path[1:], Merkle.H(leaf + path[0]).digest())
            else:
                return Merkle.verify_(root, index >> 1, path[1:], Merkle.H(path[0] + leaf).digest())

def commit( data_array ):
    return Merkle.commit_([Merkle.H(bytes(da)).digest() for da in data_array])

def open( index, data_array ):
    return Merkle.open_(index, [Merkle.H(bytes(da)).digest() for da in data_array])

def verify( root, index, path, data_element ):
    return Merkle.verify_(root, index, path, Merkle.H(bytes(data_element)).digest())