汉明距离之算法和实现总结

内容简介

汉明距离,通过比较向量每一位是否相同,求出不同位的个数。用来表示两个向量之间的相似度。

汉明距离计算的步骤,即对两个向量首先进行异或操作,然后对异或的结果的每一位bit进行统计,最后合计出有多少bit的值为1。

本文主要内容为,列举出9种计算汉明距离的算法及其C++代码的实现,并使用本文的测试方法测试得出不同算法的性能。再对不同的算法进行分析比较。

计算机在进行异或操作中,CPU的指令集可以提供多种实现。比如cpu固有指令 xor 和 SSE指令集 、AVX指令集。后两种指令集都是为了提升CPU对向量的处理速度而扩展的指令集。

汉明距离计算的后一步,计算一个变量中所有bit的值为1的个数。可以使用多种算法和实现方式来实现。比如算法上,逐位移位统计、查表法、分治法等;实现方式上,可以使用前面所说的算法,也可以使用cpu的指令popcnt直接求得。

实现算法

本文算法计算的向量为二值向量;

本文中包含算法函数的Algorithm类声明如下,不同计算汉明距离的算法的类都继承这个基类,计算汉明距离由成员函数 uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size)实现,函数输入为3个参数,其中p和q分别为两个向量的指针;size为向量维度的大小,并以64为倍数; 向量最小为64个bit。

class Algorithm {
  public:
    ~Algorithm() {}
    virtual void init() {}
    virtual std::string getName() {
      return "Undefined Algorithm Name.";
    }
    virtual uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) = 0;
};

一般算法

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      for (uint64_t i = 0; i < size; ++i) {
        uint64_t r = (*(p + i)) ^ (*(q + i));
        while (r) {
          res += r & 1;
          r = r >> 1;
        }
      }
      return res;
    }

使用gcc内建函数优化一般算法

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      for (uint64_t i = 0; i < size; ++i) {
        uint64_t r = (*(p + i)) ^ (*(q + i));
        res += __builtin_popcountll(r);
      }
      return res;
    }

查表法-按8bit查询

class HammingDistanceTable8Bit : public Algorithm {
  public:
    std::string getName() {
      return "HammingDistanceTable8Bit";
    }
    void init() {
      pop_count_table_ptr = NULL;
      pop_count_table_8bit_init(&pop_count_table_ptr);
    }
    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      for (uint64_t i = 0; i < size; ++i)
      {
        uint64_t r = (*(p + i)) ^ (*(q + i));
        res += pop_count_table_8bit(r);
      }
      return res;
    }
  private:
    uint8_t *pop_count_table_ptr; 
    void pop_count_table_8bit_init(uint8_t **pop_count_table_ptr) {
      *pop_count_table_ptr = new uint8_t[256];
      for (int i = 0; i < 256; ++i) {
        (*pop_count_table_ptr)[i] = __builtin_popcount(i);
      }  
    }

    uint64_t pop_count_table_8bit(uint64_t n) {
      int res = 0;
      uint8_t *p = (uint8_t *)&n;
      for (int i = 0; i < 8; ++i) {
        res += pop_count_table_ptr[*(p + i)];
      }
      return res;
    }
};

查表法-按16bit查询

class HammingDistanceTable16Bit : public Algorithm {
  public:
    std::string getName() {
      return "HammingDistanceTable16Bit";
    }
    void init() {
      pop_count_table_ptr = NULL;
      pop_count_table_16bit_init(&pop_count_table_ptr);
    }

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      for (uint64_t i = 0; i < size; ++i)
      {
        uint64_t r = (*(p + i)) ^ (*(q + i));
        res += pop_count_table_16bit(r);
      }
      return res;
    }
  private:
    uint8_t *pop_count_table_ptr; 

    void pop_count_table_16bit_init(uint8_t **pop_count_table_ptr) {
      *pop_count_table_ptr = new uint8_t[65536];
      for (int i = 0; i < 65536; ++i) {
        (*pop_count_table_ptr)[i] = __builtin_popcount(i);
      }  
    }

    uint64_t pop_count_table_16bit(uint64_t n) {
      int res = 0;
      uint16_t *p = (uint16_t *)&n;
      for (int i = 0; i < 4; ++i) {
        res += pop_count_table_ptr[*(p + i)];
      }
      return res;
    }
};

分治法

class HammingDistanceDivideConquer : public Algorithm {
  public:
    std::string getName() {
      return "HammingDistanceDivideConquer";
    }

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      for (uint64_t i = 0; i < size; ++i)
      {
        uint64_t r = (*(p + i)) ^ (*(q + i));
        res += pop_count_divide_conquer(r);
      }
      return res;
    }

    uint64_t pop_count_divide_conquer(uint64_t n) {
      n = (n & 0x5555555555555555) + ((n >> 1 ) & 0x5555555555555555); 
      n = (n & 0x3333333333333333) + ((n >> 2 ) & 0x3333333333333333);
      n = (n & 0x0F0F0F0F0F0F0F0F) + ((n >> 4 ) & 0x0F0F0F0F0F0F0F0F);
      n = (n & 0x00FF00FF00FF00FF) + ((n >> 8 ) & 0x00FF00FF00FF00FF);
      n = (n & 0x0000FFFF0000FFFF) + ((n >> 16) & 0x0000FFFF0000FFFF);
      n = (n & 0x00000000FFFFFFFF) + ((n >> 32) & 0x00000000FFFFFFFF);
      return n;
    }
  private:

};

改进分治法

class HammingDistanceDivideConquerOpt : public Algorithm {
  public:
    std::string getName() {
      return "HammingDistanceDivideConquerOpt";
    }

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      for (uint64_t i = 0; i < size; ++i)
      {
        uint64_t r = (*(p + i)) ^ (*(q + i));
        res += pop_count_divide_conquer_opt(r);
      }
      return res;
    }

    uint64_t pop_count_divide_conquer_opt(uint64_t n) {
      n = n - ((n >> 1)  & 0x5555555555555555); 
      n = (n & 0x3333333333333333) + ((n >> 2 ) & 0x3333333333333333);
      n = (n + (n >> 4 )) & 0x0F0F0F0F0F0F0F0F;
      n = n + (n >> 8 );
      n = n + (n >> 16);
      n = n + (n >> 32);
      return (uint64_t)(n & 0x7F);
    }
  private:

};

使用SSE指令集

class HammingDistanceSSE : public Algorithm {
  public:
    std::string getName() {
      return "HammingDistanceSSE";
    }

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      uint64_t temp_res[2] = {0, 0};
      for (uint64_t i = 0; i < size; i += 2) { 
        __m64 *x1 = (__m64*)(p);
        __m64 *x2 = (__m64*)(p + 1);
        __m64 *y1 = (__m64*)(q);
        __m64 *y2 = (__m64*)(q + 1);
        __m128i z1 = _mm_set_epi64 (*x1, *x2);
        __m128i z2 = _mm_set_epi64 (*y1, *y2);
        __m128i xor_res =  _mm_xor_si128(z1 , z2);
        _mm_store_si128((__m128i*)temp_res, xor_res);
        res += _mm_popcnt_u64(temp_res[0]);
        res += _mm_popcnt_u64(temp_res[1]);
      }
      return res;
    }
};

使用AVX2指令集

class HammingDistanceAVX : public Algorithm {
  public:
    std::string getName() {
      return "HammingDistanceAVX";
    }

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      uint64_t temp_res[4] = {0, 0, 0, 0};
      for (uint64_t i = 0; i < size - 1; i += 4) { 
        long long int *x1 = (long long int*)(p + i);
        long long int *x2 = (long long int*)(p + i + 1);
        long long int *x3 = (long long int*)(p + i + 2);
        long long int *x4 = (long long int*)(p + i + 3);
        long long int *y1 = (long long int*)(q + i);
        long long int *y2 = (long long int*)(q + i + 1);
        long long int *y3 = (long long int*)(q + i + 2);
        long long int *y4 = (long long int*)(q + i + 3);
        __m256i z1 = _mm256_set_epi64x (*x1, *x2, *x3, *x4);
        __m256i z2 = _mm256_set_epi64x (*y1, *y2, *y3, *y4);
        __m256i xor_res =  _mm256_xor_si256(z1 , z2);
        _mm256_store_si256((__m256i*)temp_res, xor_res);
        res += _mm_popcnt_u64(temp_res[0]);
        res += _mm_popcnt_u64(temp_res[1]);
        res += _mm_popcnt_u64(temp_res[2]);
        res += _mm_popcnt_u64(temp_res[3]);
      }
      return res;
    }
};

使用AVX512指令集

截止本文完成时,市场上支持avx-512指令的cpu并没有普及,但是gcc已经提供了avx-512的c/c++ 语言接口。本文先把代码实现,后续购得支持avx-512指令的cpu后,再进行测试。

class HammingDistanceAVX : public Algorithm {
  public:
    std::string getName() {
      return "HammingDistanceAVX";
    }

    uint64_t cal(const uint64_t* p, const uint64_t* q, const uint64_t size) {
      uint64_t res = 0;
      uint64_t temp_res[8] = {0, 0, 0, 0, 0, 0, 0, 0};
      for (uint64_t i = 0; i < size - 1; i += 8) { 
        long long int *x1 = (long long int*)(p + i);
        long long int *x2 = (long long int*)(p + i + 1);
        long long int *x3 = (long long int*)(p + i + 2);
        long long int *x4 = (long long int*)(p + i + 3);
        long long int *x5 = (long long int*)(p + i + 4);
        long long int *x6 = (long long int*)(p + i + 5);
        long long int *x7 = (long long int*)(p + i + 6);
        long long int *x8 = (long long int*)(p + i + 7);
        long long int *y1 = (long long int*)(q + i);
        long long int *y2 = (long long int*)(q + i + 1);
        long long int *y3 = (long long int*)(q + i + 2);
        long long int *y4 = (long long int*)(q + i + 3);
        long long int *y5 = (long long int*)(q + i + 4);
        long long int *y6 = (long long int*)(q + i + 5);
        long long int *y7 = (long long int*)(q + i + 6);
        long long int *y8 = (long long int*)(q + i + 7);
        __m512i z1 = _mm512_set_epi64x (*x1, *x2, *x3, *x4, *x5, *x6, *x7, *x8);
        __m512i z2 = _mm512_set_epi64x (*y1, *y2, *y3, *y4);
        __m512i xor_res =  _mm512_xor_si512(z1 , z2);
        _mm512_store_si512((void*)temp_res, xor_res);
        res += _mm_popcnt_u64(temp_res[0]);
        res += _mm_popcnt_u64(temp_res[1]);
        res += _mm_popcnt_u64(temp_res[2]);
        res += _mm_popcnt_u64(temp_res[3]);
        res += _mm_popcnt_u64(temp_res[4]);
        res += _mm_popcnt_u64(temp_res[5]);
        res += _mm_popcnt_u64(temp_res[6]);
        res += _mm_popcnt_u64(temp_res[7]);
      }
      return res;
    }
};

性能测试

测试代码框架类:

class AlgorithmBench {
  public:
    void init() {
      startTimer();
      vector = new uint64_t[size];
      ifstream f("sample.txt", ios::in);
      for(int i = 0; i < size;i++ ) {
        uint64_t c;
        f >> vector[i];
      }
      stopTimer();
      getInterval("load sample cost:");
      f.close();
    }
    void setSize(uint64_t size) { this->size = size;};
    void push_back(Algorithm* algorithm) { _algorithm_vector.push_back(algorithm);}

    void start() {
      for (std::vector<Algorithm*>::iterator iter = _algorithm_vector.begin();
          iter != _algorithm_vector.end();
          ++iter) {
        Algorithm* ptr = *iter;
        ptr->init();
        startTimer();
        for (int i = 0; i < size - 3; ++i) {
          ptr->cal(vector + i, vector + i + 1, 4);
        }
        stopTimer();
        getInterval(ptr->getName() + " cost:");
      }
    }

    void startTimer() {
      gettimeofday(&tv,NULL);
      start_timer = 1000000 * tv.tv_sec + tv.tv_usec;
    }
    void stopTimer() {
      gettimeofday(&tv,NULL);
      end_timer = 1000000 * tv.tv_sec + tv.tv_usec;
    }
    void getInterval(std::string prefix) {
      std::cout<<std::left<<setw(40) << prefix 
        << std::right << end_timer - start_timer<<endl;
    }

  private:
    uint64_t size;
    uint64_t *vector;
    timeval tv;
    uint64_t start_timer;
    uint64_t end_timer;
    std::vector<Algorithm*> _algorithm_vector;
};

编译指令:

g++ -msse4.2 -mavx2 -O2 -o test_hamming hamming_distance.cpp

windows-gcc 测试结果

测试环境:

Windows 7
gcc version 7.4.0

HammingDistanceBase cost:               330066
HammingDistanceBuildin cost:            326566
HammingDistanceTable8Bit cost:          2381976
HammingDistanceTable16Bit cost:         1435287
HammingDistanceDivideConquer cost:      1215243
HammingDistanceDivideConquerOpt cost:   1226745
HammingDistanceSSE cost:                972695
HammingDistanceAVX cost:                680636

Linux-gcc 测试结果

测试环境:

Ubuntu Server 19.04
gcc version 8.3.0

测试结果:

load sample cost:                       78070
HammingDistanceBase cost:               145393
HammingDistanceBuildin cost:            75905
HammingDistanceTable8Bit cost:          598789
HammingDistanceTable16Bit cost:         142502
HammingDistanceDivideConquer cost:      343414
HammingDistanceDivideConquerOpt cost:   316748
HammingDistanceSSE cost:                59322
HammingDistanceAVX cost:                115784

总结

不同平台,不同的编译器版本,测试的结果有所差异。但整体表现上使用SSE指令集的性能最好,其次是使用内建函数计算popcnt性能最优。AVX指令性能略好于SSE。性能最差的是查表法8bit。分治法性能居中。

附录

所有代码都存放在github上面: https://github.com/zuocheng-liu/code-samples/blob/master/algorithm/hamming_distance.cpp

STL容器 vector list 性能对比(附带测试方法和测试数据)

最近在重构公司的一个C++模块,逻辑里有排序、过滤等操作。开发过程中,遇到下面的一些问题:

问题和猜想

  • 单链表的归并排序线性表的快速排序的时间复杂度都是 O(nlog(n)),但在实际情况中哪个更快呢?需要测试一下。
  • 猜想的发散,stl vector 和 list(双向链表)的迭代器遍历、顺序插入、clear操作, 时间复杂度好像都是相同的,但谁更快呢?
  • 看stl代码去研究vector和list的实现,并不能得到我们最终想要的结果,测试数据会更直观。

测试总结

先写测试总结,测试方法、测试结果都在下面。

  • 迭代器遍历,list 比vector 稍块,使用for_each 比使用for更快。
  • 顺序插入, vector比list 约快3倍
  • clear操作,vector 几乎不耗时,list 要耗费好多时间,vector比list至少快1千倍以上
  • 排序, vector 大约比list 快2倍。

结合项目,对比vector和list的性能

  • 在我重构前,模块使用的是vector存储数据,排序使用的是快排(大于16小于堆),过滤则是用临时数据vector转存。
  • 在我重构后,使用list存储,排序使用list的归并排序,过滤则是直接使用list的erase操作。
  • 重构前,耗时4000微秒,重构后1500微秒,时间减少60%。
  • list除了删除操作和遍历操作,排序和插入都比vector慢,但是我使用list后,却提升了性能。原因可能有2个: 1.单个元素数据结构大,vector移动这些元素代价较大,而list移动这些元素反而代价很小;2.去掉了中间临时存储,数据的转存代价也比较大。
  • 使用list后,模块的可扩展性也变得更好了。

这次项目最大的收获就是,如果有过滤操作,优选使用list。vector的缺点如下,1.vector快速排序需要移动元素,如果元素占据空间大,移动代价也非常大。2.vector过滤需要借助中间临时存储,直接erase的代价更大。

测试方法

  • 随机生成10万数据量,由vector和list结构的变量分别存储
  • 对vector 和 list 的数据分别做如下, 1)迭代器遍历 2)顺序插入 3) clear操作 4)排序
  • 记录每种操作消耗的时间
  • 多次测试,记录典型的数据结果

测试结果

直接列出测试结果,单位是微秒

vector for_each : 3263
list for_each : 2556
vector traverse : 4783
vector num traverse : 666
list traverse : 3078
vector push_back : 8136
list push_back : 24581
list delete : 7428
vector push_back : 7329
list push_back : 22968
vector clear : 0
list clear : 13079
vector sort : 89512
list sort : 146529

(附)测试程序

Github 托管地址,保持更新:https://github.com/zuocheng-liu/code-samples/blob/master/linux-c/stl/vector_list_performance_test.cpp

#include<stdio.h>
#include<stdlib.h>
#include<sys/time.h>
#include <list>
#include <vector>
#include <iostream>

#define MAX_NUM 100000

using namespace std;

timeval tv;
uint64_t timer;


bool NumCmp(uint32_t a, uint32_t b) { return (a > b); }

void startTimer() {
    gettimeofday(&tv,NULL);
    timer = 1000000 * tv.tv_sec + tv.tv_usec;
}

uint64_t stopTimer() {
    gettimeofday(&tv,NULL);
    timer = 1000000 * tv.tv_sec + tv.tv_usec - timer;
    return timer;
}

int main() {
    vector<uint32_t> v;
    vector<uint32_t> v2;
    list<uint32_t> l;
    list<uint32_t> l2;
    for (int i = 0; i< MAX_NUM; ++i) {
        srand((int)time(0)  * i * i + 1);
        int num = rand();
        v.push_back(num);
        l.push_back(num);
    }

    // compare oper traverse
    startTimer();   
    for (vector<uint32_t>::iterator iter = v.begin(); iter != v.end(); ++ iter) {
        *iter;
    }
    cout<<"vector\t traverse\t :\t"<< stopTimer() << endl;  

    startTimer();   
    for (int i = 0 ; i < MAX_NUM; ++ i) {
        //v[i];
    }
    cout<<"vector\t num traverse\t :\t"<< stopTimer() << endl;  


    startTimer();   
    for (list<uint32_t>::iterator iter = l.begin(); iter != l.end(); ++ iter) {
    }
    cout<<"list\t traverse\t :\t"<< stopTimer() << endl;  


    // compare oper push_back   
    startTimer();   
    for (vector<uint32_t>::iterator iter = v.begin(); iter != v.end(); ++ iter) {
        v2.push_back(*iter);
    }
    cout<<"vector\t push_back\t :\t"<< stopTimer() << endl;  

    startTimer();   
    for (list<uint32_t>::iterator iter = l.begin(); iter != l.end(); ++ iter) {
        l2.push_back(*iter);
    }
    cout<<"list\t push_back\t :\t"<< stopTimer() << endl;  

    // compare oper delete
    startTimer();
    v2.clear();
    for (vector<uint32_t>::iterator iter = v2.begin(); iter != v2.end(); ++ iter) {
    //    v2.erase(iter);
    }
    //cout<<"vector\t delete\t :\t"<< stopTimer() << endl;  

    startTimer();   
    for (list<uint32_t>::iterator iter = l2.begin(); iter != l2.end(); ++ iter) {
        iter = l2.erase(iter);
    }
    cout<<"list\t delete\t :\t"<< stopTimer() << endl;  

    // compare oper push_back   
    startTimer();   
    for (vector<uint32_t>::iterator iter = v.begin(); iter != v.end(); ++ iter) {
        v2.push_back(*iter);
    }
    cout<<"vector\t push_back\t :\t"<< stopTimer() << endl;  

    startTimer();   
    for (list<uint32_t>::iterator iter = l.begin(); iter != l.end(); ++ iter) {
        l2.push_back(*iter);
    }
    cout<<"list\t push_back\t :\t"<< stopTimer() << endl;  


    // compare oper clear
    startTimer();   
    v2.clear();
    cout<<"vector\t clear      \t:\t"<< stopTimer() << endl;  

    startTimer();   
    l2.clear();
    cout<<"list\t clear     \t :\t"<< stopTimer() << endl;  

    // compare oper sort
    startTimer();   
    std::sort(v.begin(), v.end(), NumCmp);
    cout<<"vector\t sort    \t :\t"<< stopTimer() << endl;  

    startTimer();   
    l.sort(NumCmp);
    cout<<"list\t sort    \t :\t"<< stopTimer() << endl;  

    return 0;
}

判断平面上点和不规则多边形位置关系算法

本文列举5中计算点和不规则多边形位置关系的算法。

计算机算法实现方案,最好的是弧长法。

射线法 (铅垂线法、水平线法)

射线法是使用最广泛的算法,这是由于相比较其他算法而言,它不但可以正确使用在凹多边形上,而且不需要考虑精度误差问题。

该算法思想是从点出发向右水平做一条射线,计算该射线与多边形的边的相交点个数,当点不在多边形边上时,如果是奇数,那么点就一定在多边形内部,否则,在外部。

其中铅垂线法和水平线法是射线法的两个具体算法,更容易实现,在坐标计算中更方便。

面积法

面积法的思想是如果点在多边形内部或者边上,那么点与多边形所有边组成的三角形面积和等于多边形面积。多边形的面积公式可以用叉积计算。不过计算面积是会有一定误差的,需要设置精度的误差范围。

点线判断法

对于多边形,如果一个点它的所有边的左边,那么这个点一定在多边形内部。利用叉积正好可以判断点与给定边的关系,即点是在边的左边右边还是边上。

转角法

以被测点O为坐标原点,计算其与所有相邻两个多边形顶点点P[i]之间的角度和,所有角度相加后得到360度,即为点在多边形之内,若得到180度,则在多边形边上,0度则在多边形外。

弧长法

弧长法是改进后的转角法,解决了传统转角法的精度问题。算法思想是,以被测点O为坐标原点,将平面划分为4个象限,对每个多边形顶点P[i],计算其所在的象限,然后顺序访问多边形的各个顶点P[i],分析P[i]和P[i+1],有下列三种情况:

  • P[i+1]在P[i]的下一象限。此时弧长和加π/2;

  • P[i+1]在P[i]的上一象限。此时弧长和减π/2;

  • P[i+1]在Pi的相对象限。利用叉积f=y[i+1]x[i]-x[i+1]y[i]计算Op[i]与Op[i+1]的关系,若f=0,Op[i]与Op[i+1]共线,点在多边形边上;若f<0,Op[i+1]在Op[i]逆时针方向,弧长和减π;若f>0,Op[i+1]在Op[i]顺时针方向,弧长和加π。

由于顶点在原点和坐标轴上时需要特殊处理,为了准确性,应该在每次计算顶点时都用叉积判断P是否在当前边上,另外将π用整数代替可以避免浮点数的精度误差。

二进制运算技巧

假设都在x86或x86_64架构CPU上进行运算,二进制基本运算包括,加减乘除、与、或、异或、同或、移位

假设n 为 32 位整形数,取正整数n除以8的余数 : n & 0x07

假设n 为 32 位整形数,取正整数n除以16的余数 : n & 0x0F

假设n 为 32 位整形数,8位对齐 : (n + 7) & 0xFFFFFFF8)

假设n 为 32 位整形数,32位对齐 : (n + 31) & 0xFFFFFF80)

假设n 为 32 位整形数,求8位对齐填充位数 ((n + 7) & 0xFFFFFFF8) - n;

求int 最大值 : (1 << 31) - 1 或者干脆 0x7FFFFFFF

求int 最小值 : 1 << 31

判断奇偶数 : n & 1 == 1

交换a b 两个变量的值 : a ^= b ^= a ^= b;

n 乘以2的m次方 : n << m