wNAF标量乘算法介绍

一、算法解决的问题

        椭圆曲线标量乘法是椭圆曲线密码系统中最关键（也是消耗资源和能量最多）的步骤，是椭圆曲线上一个点P与一个随机的整数k的乘积，即𝑄=𝑘𝑃=𝑃+𝑃+⋯+𝑃.因此通常做法是将k进行转换以减少其中点加次数，如二进制、NAF、wNAF等标量乘算法这些算法比可以较大降低计算的次数。今天主要介绍wNAF，wNAF是一种NAF的w进制表示形式，这里的w指的是占二进制位数，比如w=1表示1位：0b1，w=2表示2位：0b11，w=3表示3位：0b111，w=4表示4位：0B1111等等依次类推。wNAF标量乘算法通过将二进制形式的标量k转化为wNAF形式，构建一条wNAF链，即，其中∈{±1,±3,⋯,±()}，𝑚表示wNAF链中的非零数字个数， 是wNAF链中每个非零数所在的位置.实现标量乘时，一般通过调用已预计算并存储的点{±𝑃,±3𝑃,⋯,±()𝑃}，其中w是窗口宽度，可极大地减少运算时的计算量.wNAF的具体实现k的求解算法如所示。

       值得注意的是预算点的值存储的都是奇数点，原因在于在标量k除以2时的是k偶数，就余数为0，只有当除到标量k为奇数才运用wNaf法则，此时余数不可能为偶数，故预存的窗口计算值不可能为偶数，因此就没必要计算偶数的计算值。

fn compute_modify_naf_weight(a:u32,w:u32){
    //a表示k的正整数值，w表示使用的位数
    let mut i :u32 =0  ;
    let mut k = a as u32;
    let mut naf:[i32; 33] = [0_i32;33];

    while(k>=1){
        if k %2 == 1 {
            let mut e = (k as i32 % 2_i32.pow(w)) as i32;
            if e >= 2_i32.pow(w-1) {
                e = e as i32 -2_i32.pow(w);
            }
            naf[i as usize] = e as i32;            
            k = (k as i32 -naf[i as usize] as i32) as u32;
        }else{
            naf[i as usize] = 0;
        }
        k = k/2;
        i=i+1;
    }
    //println!("{:?}",naf);
    let (fac_naf,__) = naf.split_at_mut(i as usize);
    fac_naf.reverse(); 
    println!("mod {i} {:?}",fac_naf);

}

二、w-NAF的原理解析

     w-NAF的原理与NAF是一样的方式，都是进行低位向高位补位，满足高位可以整除2或者。以下对w-NAF用299这个整数进行低补高进位分析：

设k=299，w=4 ，w=4表示当值k/2为奇数时，在操作k/2之后，下面跟随是连续000个，这就可以降低k*P时非零值的密度。操作过程如下：

299= '0b100101011'

==16，

第一次 k=299

299/2 =149 余1 299是奇数，可以运用w-NAF法则，

取299/16 =18 余11，说明299+(11-16)可以整整16

故第1次299/2的计算变成,余（11-16）余-5，k=(299+(11-16))/2=152

余5

第二次 k=152

152/2 = 76 余0  k=76

余0

第三次 k=76

76/2 = 38 余0 k=38

 余0

第4次 k = 38

76/2=19  余0 k = 19

 余0

 第5次 k=19

19/2=9 余1，为奇数，运行w-NAF法则，19%16 = 3，此处是余3，就不补进位，直接算为余数，k = (19 -3)/2=8。需要注意的是否补高位取决于余数大小，如果余数大于等于16/2=8，大于就补进位，小于就补需要补进位，具体公式就是 余数>=,向高位进，余数就变成了(k% - )，否则就保持k%。

余3

第6次 k = 8

8/2=4 余0

余0

第7次 k = 4

4/2=2 余0

余0

第8次 k=2

2/2 = 1 余0

余0

第9次 k =1

1/2 = 0 余1

转化w-NAF的分解链为

w-NAF = 10003000-5

转化为299的计算公式为 -5* + 3*+1* = -5+ 3*16 + 256 = 299

三、kP的计算

import ecdsa
from ecdsa.ellipticcurve import Point,INFINITY

################################计算b的w-NAF链####################################
def comput_coeff_wnaf(b,w):
    k = b
    coeff_wnaf_val = []
    i = 0
    zeroes = 0
    used = 0
    while k>0:
        while k%2==0 :
            zeroes+=1
            k = int(k/2)
        word = k & ((1<<w) -1)
        k = int(k / pow(2,w))
        if (word & (1<<(w-1)))>0:
            k +=1
            for i in range(0,zeroes):
                coeff_wnaf_val.append(0)
                used+=1
            coeff_wnaf_val.append(word -(1<<w))
            used+=1
        else:
            for i in range(0,zeroes):
                coeff_wnaf_val.append(0)
                used+=1
            coeff_wnaf_val.append(word)
            used+=1
        zeroes = w-1

    print('coeff_wnaf_val=', coeff_wnaf_val, ' len %d' % (used))
    return coeff_wnaf_val

###############预计算 [1P,3P,5P,7P...(2n+1)P]###################
def precomP(w):
    pre = []
    count =0
    i =1
    max = 1<<(w-2)#因为只需要存储2^(w-1)的一半值，偶数值不需要存储，故预计算数组存储值需要2^(w-2)
    # 初始化P点
    x = 0x79f823642ddb09f87acae2126775613542af3fe156d1b936a813cc368ff626bd
    y = 0x316b223447295841d0d025b5c76b87b42c6ddf879ec812bcf9eaeb538565a6b2
    P: Point = Point(ecdsa.curves.SECP256k1.curve, x, y)
    pre.append(P)
    count +=1
    D: Point = P.double()
    for i in range(1,max):
        X: point = D + pre[i-1]
        pre.append(X)
        count +=1
    """
    for i in range(0,max):
        t:Point = pre[i]
        q_x = t.x()
        q_y = t.y()
        print('i=%d'%i, 'q_x', q_x.to_bytes(32, byteorder='big').hex())
        print('i=%d'%i, 'q_y', q_y.to_bytes(32, byteorder='big').hex())   
    """
    return pre
################使用w-NAF计算bP###############
def compute_wNAF_kP(b,w):

    Q: Point = INFINITY
    wnaf_arr = comput_coeff_wnaf(b,w)
    pre_P_arr = precomP(w)
    
   # print('binarray=%s' % binarray, ' len %d' % bin_len)

    for i in reversed(wnaf_arr):
        Q = 2 * Q
        if i!=0:
            if i >0 :
                Q = Q + pre_P_arr[int((i-1)/2)]
            if i < 0:
                Q = Q +pre_P_arr[int((abs(i)-1)/2)].__neg__() # 即 Q= Q -P

    q_x = Q.x()
    q_y = Q.y()

    print('q_x', q_x.to_bytes(32, byteorder='big').hex())
    print('q_y', q_y.to_bytes(32, byteorder='big').hex())

    None
#################################计算bP##########################################
w = 4 #窗口宽度4
b = 6 #标量值
print('b=%d'%b,",w=%d"%w)
compute_wNAF_kP(b,w)

#comput_coeff_wnaf(b,w)
#precomP(w)

执行结果如下：