算法 – Baocheng's home page

2018年1月4日

一个卓有成效的汇编优化范例–使用SSE2指令优化进制转化

我的一个感兴趣的编程方向是大数计算，因此用汇编语言写了很多大数计算方面的小程序，上周突然想出一个使用SSE2指令将整数转为16进制字符串的好主意，遂付诸实现。原以为至多可提速500%，那知测试后发现，相对于最初的C语言版本，速度竟提高20倍以上，兴奋之余，遂有了这篇博客文章。

这个程序主要示范将64bit一个整数转化为16进制字符串的功能，功能和算法都比较简单。我相信许多人都写过类似的程序，但不知有没有人尝试去你优化它。这个示范程序包括3个C语言版和1个使用SSE2指令的汇编语言版。下面我们给出代码和说明。

先看这个函数最初的版本，UINT64_2_hexString_c1，为了性能起见，我们使用 __fastcall 函数约定，__fastcall 接口的函数使用寄存器来传递参数，免除了调用时压栈的开销，而且被调函数可以省去保存/恢复寄存器等指令。

//这是C语言普通版,直接使用64位整数逻辑指令和算术指令
void __fastcall UINT64_2_hexString_c1(UINT64 *p, char *buff)
{
  UINT64 x=*p;
  int i;
  for (i=15;i>=0;i--)
  {
    char c= (x & 0xf)+'0';
    if ( c>'9')
      c+=7;
    buff[i]=c;
    x>>=4;
  }
  buff[16]=0;
}

//这是C语言普通版,直接使用64位整数逻辑指令和算术指令

void __fastcall UINT64_2_hexString_c1(UINT64 *p, char *buff)

{

UINT64 x=*p;

int i;

for (i=15;i>=0;i--)

{

char c= (x & 0xf)+'0';

if ( c>'9')

c+=7;

buff[i]=c;

x>>=4;

}

buff[16]=0;

}

上面这个函数虽然简单，然而速度却仍不理想。我们知道，在32位运行环境，对64位整数计算的C语言语句要翻译成多条指令，故速度较慢，下面这个版本使用完全的32位整数处理，故速度快于上面的版本。

void __fastcall UINT64_2_hexString_c2(UINT64 *p, char *buff)
{
	DWORD *pDW=(DWORD *)p;
	DWORD x;
	int i;
	for (x=pDW[1],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		if ( c>'9')
			c+=7;
		buff[i]=c;
		x>>=4;
	}
	for (x=pDW[0],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		if ( c>'9')
			c+=7;
		buff[8+i]=c;
		x>>=4;
	}
	buff[16]=0;
}

void __fastcall UINT64_2_hexString_c2(UINT64 *p, char *buff)

{

DWORD *pDW=(DWORD *)p;

DWORD x;

int i;

for (x=pDW[1],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

if ( c>'9')

c+=7;

buff[i]=c;

x>>=4;

}

for (x=pDW[0],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

if ( c>'9')

c+=7;

buff[8+i]=c;

x>>=4;

}

buff[16]=0;

}

上面这个函数首先将64位整数地址转化32位整数的地址，然后使用32位整数运算。代码是复杂了，但速度更快了。尽管如此，程序仍有优化的空间。

我们注意到，上面的函数包括2个循环，在每个循环中又有一个if语句，站在汇编语言视角，循环和if语句都是分支语句，可编译成比较跳转指令对。分支对CPU是个麻烦事儿，由于现代CPU普遍采用管线技术，在执行当前指令时，后面的指令已经被取到甚至译码完成。CPU遇到分支时，就需要预测那个分支最可能被执行到，从而将最可能被执行到的那个分支的代码加载到管线。当分支预测成功时，所有的指令可流畅地无停顿的执行。一旦分支预测失败，则不得不将管线中已加载测的指令全部丢弃，重新从正确的分支点取指和译码。分支预测的技术很复杂，完整的讲述需要一本书的内容。我们这里仅作简单介绍。分支预测的的实现通常是这样的，在首次遇到分支时，执行非跳转分支。在每次执行分支指令时，将实际执行情况（执行那个分支）存入历史记录。以后再遇到这个分支时，则可以根据历史记录来判断那个分支最可能被执行到。最简单的一种判断算法是，总是预测执行概率比较高的那个分支，这种分支预测方案对循环引起的分支最有效。比如一个循环次数为n的for循环，前n次总是从循环体底部跳转到头，只有最后一次循环不跳转，换言之，跳转分支执行的概率远高于非跳转分支，故CPU总是预测跳转分支。就上面的例子而言，两个循环都是固定次数的循环且循环次数很少。在这种情况下，编译器可使用循环展开的方法来消除分支，但是对于语句”if ( c>’9′ ) c+=7″ 这样的分支，分支预测技术很难奏效，0-15之间的随机数，大于9的概率37.5%，即使CPU总是预测<=’9’的那个分支，也有37.5%的预测失败的概率，分支预测失败，CPU需要付出相当的代价，需要几个甚至10个额外的周期。

我的下一个版本用到的技术就是分支消除技术，通过消除分支来提高函数执行速度，为了消除分支，不得不使用额外的语句，虽然代码变多了，但函数执行速度大大加快。

//这是使用消除分支技术的C语言版，在i7-4700HQ，速度是上一个版本的2.6倍
void __fastcall UINT64_2_hexString_c3(UINT64 *p, char *buff)
{
	DWORD *pDW=(DWORD *)p;
	DWORD x;
	int i;
	for (x=pDW[1],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		char mask= 0 -( c>'9');	//the flag is 0xff when c>'9'
		buff[i]= c + ( mask & 7);
		x>>=4;
	}
	for (x=pDW[0],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		char mask= 0 -( c>'9');	//the flag is 0xff when c>'9'
		buff[i+8]= c + ( mask & 7);
		x>>=4;
	}
	buff[16]=0;
}

//这是使用消除分支技术的C语言版，在i7-4700HQ，速度是上一个版本的2.6倍

void __fastcall UINT64_2_hexString_c3(UINT64 *p, char *buff)

{

DWORD *pDW=(DWORD *)p;

DWORD x;

int i;

for (x=pDW[1],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

char mask= 0 -( c>'9'); //the flag is 0xff when c>'9'

buff[i]= c + ( mask & 7);

x>>=4;

}

for (x=pDW[0],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

char mask= 0 -( c>'9'); //the flag is 0xff when c>'9'

buff[i+8]= c + ( mask & 7);

x>>=4;

}

buff[16]=0;

}

下面是终极版本，使用SSE2指令的汇编版，直接使用ALU指令编程的优化作用有限，甚至不如编译器。我们这里直接使用SSE2指令，SSE2指令主要使用XMM寄存器来工作，1个XMM寄存器可看成是4个DWORD，8个WORD，16个BYTE，1个UINT64位数转化为字符串后有16个字符，可全部装在一个XMM寄存器中，所以这个工作正好适合用SSE2指令来做。下面的汇编版的代码，用到几个16字节数组，要求16直接对齐，尽管在汇编中也可以定义16字节对齐，但我们这里把数组的定义放在C文件中，放在C文件中的好处是易于扩展，比如可在C文件中定义32字节对齐，我曾尝试在汇编文件中定义32字节对齐时，但汇编器总是报错。这里给出C语言中的常数数组的定义。

#include <stdio.h>
#include <stdlib.h>

typedef unsigned char  BYTE;

#if defined(__GNUC__)    // GCC
    #define INTRIN_ALIGN(n)    __attribute__((aligned(n)))
#else    
    #define INTRIN_ALIGN(n)    __declspec(align(n))
#endif    // #if defined(__GNUC__)    // GCC

INTRIN_ALIGN(16) BYTE i256_num_0s[16]=
{
	'0','0','0','0','0','0','0','0',
	'0','0','0','0','0','0','0','0',
};

INTRIN_ALIGN(16) BYTE i256_num_9s[16]=
{
	'9','9','9','9','9','9','9','9',
	'9','9','9','9','9','9','9','9',
};

INTRIN_ALIGN(16) BYTE i256_num_full7[16]=
{
	7,7,7,7,7,7,7,7,
	7,7,7,7,7,7,7,7,
};

#include <stdio.h>

#include <stdlib.h>

typedef unsigned char BYTE;

#if defined(__GNUC__) // GCC

#define INTRIN_ALIGN(n) __attribute__((aligned(n)))

#else

#define INTRIN_ALIGN(n) __declspec(align(n))

#endif // #if defined(__GNUC__) // GCC

INTRIN_ALIGN(16) BYTE i256_num_0s[16]=

{

'0','0','0','0','0','0','0','0',

};

INTRIN_ALIGN(16) BYTE i256_num_9s[16]=

{

'9','9','9','9','9','9','9','9',

};

INTRIN_ALIGN(16) BYTE i256_num_full7[16]=

{

7,7,7,7,7,7,7,7,

};

下面是微软汇编器ml.exe格式的汇编语言源代码。

 	.686P
	.XMM
	.model	flat

_DATA	SEGMENT
	EXTRN _i256_num_0s:BYTE
	EXTRN _i256_num_9s:BYTE
	EXTRN _i256_num_full7:BYTE
_DATA	ENDS

PUBLIC @UINT64_2_hexString_sse2@8

_TEXT	SEGMENT

; The XMM register definition for function _UINT64_2_hexString_sse2
XMM_R_DB_STR_0  TEXTEQU <XMM2>
XMM_R_DB_STR_9  TEXTEQU <XMM3>
XMM_R_DB_7	TEXTEQU <XMM4>
XMM_R_TMP	TEXTEQU <XMM1>

@UINT64_2_hexString_sse2@8 PROC
	
	movq    XMM0, mmword ptr [ecx]	
	PSHUFD  XMM0, XMM0, 11001101b		;now XMM0 contain 2 QWORD and low 32bits in every QWORD is valid

	; We denote the value of xmm register with a string, "N" means a byte, "0" means a byte whose value is 0,
	; We use low-bytes first order,
  	; Now, the value of XMM0 is (NNNN0000,NNNN0000),contain 2 QWORD and contain 32 bits in every qword

	; the first round transform, 2 QWORD => 4 DWORD
	MOVDQA  XMM_R_TMP, XMM0
	
	PSLLQ  XMM0, 48				;toward high shift 48 bits
	; now  value of  XMM0 is (000000NN, 000000NN), now bit0-bit15 of every original QWORD is in XMM0 
	
	PSRLQ  XMM_R_TMP, 16		;toward low shift 16 bits
	; now  value of  XMM1 is (NN000000,NN000000), bit16-bit31 of every original QWORD is in XMM1
	
	PSRLQ  XMM0, 16				;toward low shift 16 bits
	; now  value of  XMM0 is  (0000NN00,0000NN00), now bit0-bit15 of every original QWORD is in XMM0 
	
	POR		XMM0, XMM_R_TMP		;merge bit0-bit15 and bit16-bit31
	; now  value of  XMM0 is (NN00,NN00,NN00,NN00), contain 4 DWORD, and 16 bits in every DWORD is valid
	 
	
	; the second round transform, 4 DWORD  => 8 WORD
	MOVDQA  XMM_R_TMP, XMM0
	PSLLD  XMM0, 24				;toworad high shift 24 bits
	; now  value of  XMM0 is (000N,000N,000N,000N), the bit0-bit7 of every original DWORD is in XMM0 
	
	PSRLD  XMM_R_TMP, 8			;toworad low shift 8 bits
	; now  value of  XMM0 is (N000,N000,N000,N000), the bit8-bit15 of every original DWORD is in XMM1

	PSRLD  XMM0, 8				;toworad low shift 8 bits
	; now  value of  XMM0 is (00N0,00N0,00N0,00N0), the bit0-bit7 of every original DWORD is in XMM0 

	POR		XMM0, XMM_R_TMP		;merge bit0-bit7 and bit8-bit16
	; now  value of  XMM0 is (N0,N0,N0,N0,N0,N0,N0,N0), contain 8 WORD, and 8 bits in every WORD is valid
	

	; the third round transform, 8 WORD  => 16 BYTE
	MOVDQA  XMM_R_TMP, XMM0
	PSLLW  XMM0, 12				;toworad high shift 12 bits
	; the bit0-bit3 of every original WORD is in XMM0 
	
	PSRLW  XMM_R_TMP, 4			;toworad low shift 4 bits
	; now  the bit4-bit7 or every original WORD is in XMM1

	PSRLW  XMM0, 4				; toworad low shift 4 bits
	; the bit0-bit3 of every original WORD is in XMM0 

	POR		XMM0, XMM_R_TMP		;merge bit0-bit3 and bit4-7
	; now  value of  XMM0 is (N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,) contain 16 BYTE, the value range of every byte is 0-15
	
	; pre-load array to XMM  register
	MOVDQA  XMM_R_DB_STR_0, xmmword ptr _i256_num_0s
	MOVDQA  XMM_R_DB_STR_9, xmmword ptr _i256_num_9s
	
	POR     XMM0, XMM_R_DB_STR_0		;now the BYTE[i] of XMM0 is '0' to '0'+15
	MOVDQA  XMM_R_TMP, XMM0			 ;now the BYTE[i] of XMM_R_TMP is '0' to '0'+15
	
	MOVDQA  XMM_R_DB_7,	xmmword ptr _i256_num_full7		;now XMM_R_DB_7 is full 7

	PCMPGTB  XMM_R_TMP, XMM_R_DB_STR_9	; BYTE[i] of XMM_R_TMP (0<=i<15) is -1 if BYTE[i]>'9'
	PAND     XMM_R_TMP, XMM_R_DB_7		; BYTE[i] of XMM_R_TMP (0<=i<15) is 7, if BYTE[i]>'9'

	PADDB    XMM0, XMM_R_TMP		;final result,now the BYTE[i] of XMM0 is '0'-'9' or 'A' to 'F'
	 
	MOVDQU   xmmword ptr [edx], XMM0
	mov    byte ptr [edx+16], 0
	ret	0
@UINT64_2_hexString_sse2@8 ENDP

_TEXT	ENDS

END

.686P

.XMM

.model flat

_DATA SEGMENT

EXTRN _i256_num_0s:BYTE

EXTRN _i256_num_9s:BYTE

EXTRN _i256_num_full7:BYTE

_DATA ENDS

PUBLIC @UINT64_2_hexString_sse2@8

_TEXT SEGMENT

; The XMM register definition for function _UINT64_2_hexString_sse2

XMM_R_DB_STR_0 TEXTEQU <XMM2>

XMM_R_DB_STR_9 TEXTEQU <XMM3>

XMM_R_DB_7 TEXTEQU <XMM4>

XMM_R_TMP TEXTEQU <XMM1>

@UINT64_2_hexString_sse2@8 PROC

movq XMM0, mmword ptr [ecx]

PSHUFD XMM0, XMM0, 11001101b ;now XMM0 contain 2 QWORD and low 32bits in every QWORD is valid

; We denote the value of xmm register with a string, "N" means a byte, "0" means a byte whose value is 0,

; We use low-bytes first order,

; Now, the value of XMM0 is (NNNN0000,NNNN0000),contain 2 QWORD and contain 32 bits in every qword

; the first round transform, 2 QWORD => 4 DWORD

MOVDQA XMM_R_TMP, XMM0

PSLLQ XMM0, 48 ;toward high shift 48 bits

; now value of XMM0 is (000000NN, 000000NN), now bit0-bit15 of every original QWORD is in XMM0

PSRLQ XMM_R_TMP, 16 ;toward low shift 16 bits

; now value of XMM1 is (NN000000,NN000000), bit16-bit31 of every original QWORD is in XMM1

PSRLQ XMM0, 16 ;toward low shift 16 bits

; now value of XMM0 is (0000NN00,0000NN00), now bit0-bit15 of every original QWORD is in XMM0

POR XMM0, XMM_R_TMP ;merge bit0-bit15 and bit16-bit31

; now value of XMM0 is (NN00,NN00,NN00,NN00), contain 4 DWORD, and 16 bits in every DWORD is valid

; the second round transform, 4 DWORD => 8 WORD

MOVDQA XMM_R_TMP, XMM0

PSLLD XMM0, 24 ;toworad high shift 24 bits

; now value of XMM0 is (000N,000N,000N,000N), the bit0-bit7 of every original DWORD is in XMM0

PSRLD XMM_R_TMP, 8 ;toworad low shift 8 bits

; now value of XMM0 is (N000,N000,N000,N000), the bit8-bit15 of every original DWORD is in XMM1

PSRLD XMM0, 8 ;toworad low shift 8 bits

; now value of XMM0 is (00N0,00N0,00N0,00N0), the bit0-bit7 of every original DWORD is in XMM0

POR XMM0, XMM_R_TMP ;merge bit0-bit7 and bit8-bit16

; now value of XMM0 is (N0,N0,N0,N0,N0,N0,N0,N0), contain 8 WORD, and 8 bits in every WORD is valid

; the third round transform, 8 WORD => 16 BYTE

MOVDQA XMM_R_TMP, XMM0

PSLLW XMM0, 12 ;toworad high shift 12 bits

; the bit0-bit3 of every original WORD is in XMM0

PSRLW XMM_R_TMP, 4 ;toworad low shift 4 bits

; now the bit4-bit7 or every original WORD is in XMM1

PSRLW XMM0, 4 ; toworad low shift 4 bits

; the bit0-bit3 of every original WORD is in XMM0

POR XMM0, XMM_R_TMP ;merge bit0-bit3 and bit4-7

; now value of XMM0 is (N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,) contain 16 BYTE, the value range of every byte is 0-15

; pre-load array to XMM register

MOVDQA XMM_R_DB_STR_0, xmmword ptr _i256_num_0s

MOVDQA XMM_R_DB_STR_9, xmmword ptr _i256_num_9s

POR XMM0, XMM_R_DB_STR_0 ;now the BYTE[i] of XMM0 is '0' to '0'+15

MOVDQA XMM_R_TMP, XMM0 ;now the BYTE[i] of XMM_R_TMP is '0' to '0'+15

MOVDQA XMM_R_DB_7, xmmword ptr _i256_num_full7 ;now XMM_R_DB_7 is full 7

PCMPGTB XMM_R_TMP, XMM_R_DB_STR_9 ; BYTE[i] of XMM_R_TMP (0<=i<15) is -1 if BYTE[i]>'9'

PAND XMM_R_TMP, XMM_R_DB_7 ; BYTE[i] of XMM_R_TMP (0<=i<15) is 7, if BYTE[i]>'9'

PADDB XMM0, XMM_R_TMP ;final result,now the BYTE[i] of XMM0 is '0'-'9' or 'A' to 'F'

MOVDQU xmmword ptr [edx], XMM0

mov byte ptr [edx+16], 0

ret 0

@UINT64_2_hexString_sse2@8 ENDP

_TEXT ENDS

END

上面我没有将英文注释翻译成中文。这是因为，对于汇编代码，高手不用讲，初学者不会看，故这里就不再给出更多的说明了。

下面给出主程序中的全部代码。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/*
这个程序示范使用优化技术来实现将64位整数转化为16进制字符串。
包括3个C语言的版本和一个使用SSE2指令的汇编版本。
经测试，在我的电脑上（I7-4700HQ），汇编版本的速度是C语言普通版20倍，是分支消除技术版本的7.5倍。

作者：Baocheng Liang
完成日期：2015-7-30, 版权所有。
*/

#define ARRAY_LEN   4096
#define LOOP_COUNT	2048

typedef unsigned long long UINT64;
typedef unsigned long DWORD;
typedef unsigned short WORD;
typedef unsigned char BYTE;

UINT64 g_nums[ARRAY_LEN];
char g_buff[ARRAY_LEN*16+16];

extern double currTime(); //使用高精度计时器

extern void __fastcall UINT64_2_hexString_sse2(UINT64 *p, char *buff);		//使用SSE2指令的汇编版本

extern void __fastcall UINT64_2_hexString_c1(UINT64 *p, char *buff);		//最普通的C语言版本

extern void __fastcall UINT64_2_hexString_c2(UINT64 *p, char *buff);		//最普通32位整数的的C语言版本

extern void __fastcall UINT64_2_hexString_c3(UINT64 *p, char *buff);		//使用分支消除技术的C语言版本

typedef void ( __fastcall *lpfn_UInt64_2_hexString)(UINT64 *p, char *buff);


//这是C语言普通版,直接使用64位整数逻辑指令和算术指令
void __fastcall UINT64_2_hexString_c1(UINT64 *p, char *buff)
{
	UINT64 x=*p;
	int i;
	for (i=15;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		if ( c>'9')
			c+=7;
		buff[i]=c;
		x>>=4;
	}
	buff[16]=0;
}

//这是C语言改进版，使用32位整数逻辑指令和算术指令
void __fastcall UINT64_2_hexString_c2(UINT64 *p, char *buff)
{
	DWORD *pDW=(DWORD *)p;
	DWORD x;
	int i;
	for (x=pDW[1],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		if ( c>'9')
			c+=7;
		buff[i]=c;
		x>>=4;
	}
	for (x=pDW[0],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		if ( c>'9')
			c+=7;
		buff[8+i]=c;
		x>>=4;
	}
	buff[16]=0;
}

//这是使用消除分支技术的C语言版，在i7-4700HQ，速度是上一个版本的2.6倍
void __fastcall UINT64_2_hexString_c3(UINT64 *p, char *buff)
{
	DWORD *pDW=(DWORD *)p;
	DWORD x;
	int i;
	for (x=pDW[1],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		char mask= 0 -( c>'9');	//the flag is 0xff when c>'9'
		buff[i]= c + ( mask & 7);
		x>>=4;
	}
	for (x=pDW[0],i=7;i>=0;i--)
	{
		char c= (x & 0xf)+'0';
		char mask= 0 -( c>'9');	//the flag is 0xff when c>'9'
		buff[i+8]= c + ( mask & 7);
		x>>=4;
	}
	buff[16]=0;
}

void test_UINT64_2_hexString(lpfn_UInt64_2_hexString fp, const char *funcName)
{
	UINT64 arr[]={
		0x0123456789abcdef,
		0x02468ACE13579BDF,
		0xaaaaaaaaaaaaaaaa,
		0xffffffffffffffff,
	};
	int i;
	char buff[17*4];
	memset(buff,0,sizeof(buff));
	printf("Test function %s\n",funcName);
	for (i=0;i<4;i++)
	{
		fp(arr+i,buff+i*17);
		printf("%s\n", buff+i*17);
	}
}

void perf_UINT64_2_hexString(lpfn_UInt64_2_hexString fp,const char *funcName)
{
	int i,j;
	UINT64 x;
	char *p;
	double t;

	//初始化全局数组g_nums 
	for (i=0;i<ARRAY_LEN;i++)
	{
		p=(char *)(&x);
		for (j=0;j<8;j++)
			p[j]= (rand() & 0xff);//get a 64bit random number
		g_nums[i]=x;
	}
	
	t=currTime();	//计时开始
	g_buff[0]=0;
	for ( i=0;i<LOOP_COUNT;i++)
	{
		p=g_buff;
		for (j=0;j<ARRAY_LEN;j++)
		{
			fp( g_nums+j,p);
			p+=16;
		}
	}
	t=(currTime()-t)*1000000000;	//转化时间到纳秒
	printf("It take %.2f ns to run function %s\n",t/(LOOP_COUNT*ARRAY_LEN),funcName);
	printf("strlen(buff) is %d\n",strlen(g_buff)); 
	
}

//功能测试
void test_function()
{
	test_UINT64_2_hexString(UINT64_2_hexString_c1,  "UINT64_2_hexString_c1");
	test_UINT64_2_hexString(UINT64_2_hexString_c2,  "UINT64_2_hexString_c2");
	test_UINT64_2_hexString(UINT64_2_hexString_c3,  "UINT64_2_hexString_c3");
	test_UINT64_2_hexString(UINT64_2_hexString_sse2,"UINT64_2_hexString_sse2");
}

//性能测试
void perf_function()
{
	perf_UINT64_2_hexString(UINT64_2_hexString_c1,  "UINT64_2_hexString_c1");
	perf_UINT64_2_hexString(UINT64_2_hexString_c2,  "UINT64_2_hexString_c2");
	perf_UINT64_2_hexString(UINT64_2_hexString_c3,  "UINT64_2_hexString_c3");
	perf_UINT64_2_hexString(UINT64_2_hexString_sse2,"UINT64_2_hexString_sse2");
}

int main(int argc, char* argv[])
{
	test_function();
	perf_function();
	return 0;
}

补充：

   这里给出跨平台的计时函数currTime的代码。

#if defined(_WIN32)

#include <windows.h>

static LARGE_INTEGER freq;

static BOOL initFreq()
{
	BOOL ret;
	if ( !QueryPerformanceFrequency( &freq) )
	{	ret=FALSE;	}
	else
	{	ret=TRUE;	}
	return ret;
}

double currTime() //使用高精度计时器
{	
	LARGE_INTEGER performanceCount;
	BOOL result;
	double time=0.0;
	BOOL bRet=TRUE;

	if (freq.QuadPart==0)
	{
		bRet=initFreq();
	}
	
	if (bRet)
	{
		result=QueryPerformanceCounter(  &performanceCount );
		time= performanceCount.HighPart * 4294967296.0 + performanceCount.LowPart;
		time=time / (   freq.HighPart * 4294967296.0 + freq.LowPart);
	}
	return time;
}

#elif defined(__linux__) 
#include <sys/time.h>
#include <stdio.h>
#include <stdlib.h>

double currTime()
{	
	struct timeval tv;
	gettimeofday(&tv, NULL);
	return (double)(tv.tv_sec) + (double)(tv.tv_usec)/1000000.00;
}
#else
 #error do not support this complier
#endif

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

这个程序示范使用优化技术来实现将64位整数转化为16进制字符串。

包括3个C语言的版本和一个使用SSE2指令的汇编版本。

经测试，在我的电脑上（I7-4700HQ），汇编版本的速度是C语言普通版20倍，是分支消除技术版本的7.5倍。

作者：Baocheng Liang

#define ARRAY_LEN 4096

#define LOOP_COUNT 2048

typedef unsigned long long UINT64;

typedef unsigned long DWORD;

typedef unsigned short WORD;

typedef unsigned char BYTE;

UINT64 g_nums[ARRAY_LEN];

char g_buff[ARRAY_LEN*16+16];

extern double currTime(); //使用高精度计时器

extern void __fastcall UINT64_2_hexString_sse2(UINT64 *p, char *buff); //使用SSE2指令的汇编版本

extern void __fastcall UINT64_2_hexString_c1(UINT64 *p, char *buff); //最普通的C语言版本

extern void __fastcall UINT64_2_hexString_c2(UINT64 *p, char *buff); //最普通32位整数的的C语言版本

extern void __fastcall UINT64_2_hexString_c3(UINT64 *p, char *buff); //使用分支消除技术的C语言版本

typedef void ( __fastcall *lpfn_UInt64_2_hexString)(UINT64 *p, char *buff);

//这是C语言普通版,直接使用64位整数逻辑指令和算术指令

void __fastcall UINT64_2_hexString_c1(UINT64 *p, char *buff)

{

UINT64 x=*p;

int i;

for (i=15;i>=0;i--)

{

char c= (x & 0xf)+'0';

if ( c>'9')

c+=7;

buff[i]=c;

x>>=4;

}

buff[16]=0;

}

//这是C语言改进版，使用32位整数逻辑指令和算术指令

void __fastcall UINT64_2_hexString_c2(UINT64 *p, char *buff)

{

DWORD *pDW=(DWORD *)p;

DWORD x;

int i;

for (x=pDW[1],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

if ( c>'9')

c+=7;

buff[i]=c;

x>>=4;

}

for (x=pDW[0],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

if ( c>'9')

c+=7;

buff[8+i]=c;

x>>=4;

}

buff[16]=0;

}

//这是使用消除分支技术的C语言版，在i7-4700HQ，速度是上一个版本的2.6倍

void __fastcall UINT64_2_hexString_c3(UINT64 *p, char *buff)

{

DWORD *pDW=(DWORD *)p;

DWORD x;

int i;

for (x=pDW[1],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

char mask= 0 -( c>'9'); //the flag is 0xff when c>'9'

buff[i]= c + ( mask & 7);

x>>=4;

}

for (x=pDW[0],i=7;i>=0;i--)

{

char c= (x & 0xf)+'0';

char mask= 0 -( c>'9'); //the flag is 0xff when c>'9'

buff[i+8]= c + ( mask & 7);

x>>=4;

}

buff[16]=0;

}

void test_UINT64_2_hexString(lpfn_UInt64_2_hexString fp, const char *funcName)

{

UINT64 arr[]={

0x0123456789abcdef,

0x02468ACE13579BDF,

0xaaaaaaaaaaaaaaaa,

0xffffffffffffffff,

};

int i;

char buff[17*4];

memset(buff,0,sizeof(buff));

printf("Test function %s\n",funcName);

for (i=0;i<4;i++)

{

fp(arr+i,buff+i*17);

printf("%s\n", buff+i*17);

}

void perf_UINT64_2_hexString(lpfn_UInt64_2_hexString fp,const char *funcName)

{

int i,j;

UINT64 x;

char *p;

double t;

//初始化全局数组g_nums

for (i=0;i<ARRAY_LEN;i++)

{

p=(char *)(&x);

for (j=0;j<8;j++)

p[j]= (rand() & 0xff);//get a 64bit random number

g_nums[i]=x;

}

t=currTime(); //计时开始

g_buff[0]=0;

for ( i=0;i<LOOP_COUNT;i++)

{

p=g_buff;

for (j=0;j<ARRAY_LEN;j++)

{

fp( g_nums+j,p);

p+=16;

}

t=(currTime()-t)*1000000000; //转化时间到纳秒

printf("It take %.2f ns to run function %s\n",t/(LOOP_COUNT*ARRAY_LEN),funcName);

printf("strlen(buff) is %d\n",strlen(g_buff));

}

//功能测试

void test_function()

{

test_UINT64_2_hexString(UINT64_2_hexString_c1, "UINT64_2_hexString_c1");

test_UINT64_2_hexString(UINT64_2_hexString_c2, "UINT64_2_hexString_c2");

test_UINT64_2_hexString(UINT64_2_hexString_c3, "UINT64_2_hexString_c3");

test_UINT64_2_hexString(UINT64_2_hexString_sse2,"UINT64_2_hexString_sse2");

}

//性能测试

void perf_function()

{

perf_UINT64_2_hexString(UINT64_2_hexString_c1, "UINT64_2_hexString_c1");

perf_UINT64_2_hexString(UINT64_2_hexString_c2, "UINT64_2_hexString_c2");

perf_UINT64_2_hexString(UINT64_2_hexString_c3, "UINT64_2_hexString_c3");

perf_UINT64_2_hexString(UINT64_2_hexString_sse2,"UINT64_2_hexString_sse2");

}

int main(int argc, char* argv[])

{

test_function();

perf_function();

return 0;

}

补充：

这里给出跨平台的计时函数currTime的代码。

#if defined(_WIN32)

#include <windows.h>

static LARGE_INTEGER freq;

static BOOL initFreq()

{

BOOL ret;

if ( !QueryPerformanceFrequency( &freq) )

{ ret=FALSE; }

else

{ ret=TRUE; }

return ret;

}

double currTime() //使用高精度计时器

{

LARGE_INTEGER performanceCount;

BOOL result;

double time=0.0;

BOOL bRet=TRUE;

if (freq.QuadPart==0)

{

bRet=initFreq();

}

if (bRet)

{

result=QueryPerformanceCounter( &performanceCount );

time= performanceCount.HighPart * 4294967296.0 + performanceCount.LowPart;

time=time / ( freq.HighPart * 4294967296.0 + freq.LowPart);

}

return time;

}

#elif defined(__linux__)

#include <sys/time.h>

#include <stdio.h>

#include <stdlib.h>

double currTime()

{

struct timeval tv;

gettimeofday(&tv, NULL);

return (double)(tv.tv_sec) + (double)(tv.tv_usec)/1000000.00;

}

#else

#error do not support this complier

#endif

关于编译：本程序使用C语言和汇编语言混合编程。C语言源文件直接加到VC工程中即可。汇编语言使用VC中自带的汇编器ml.exe来编译。

方法

1. 将汇编文件加入到VC工程中。

2.选中文件，右键属性菜单，Item type选Custom Build Tool. 在Command Line一栏: 输入“ml /coff /c %(FullPath)”, 在Outputs 一栏输入”%(Filename).obj;%(Outputs)”

测试结果

函数	C1	C2	C3	SSE4
时间（纳秒）	67.95	52.71	20.16	3.18
相对速度	100%	129%	335%	2138%

2017年12月25日2018年1月4日

重磅消息–一种新的任意精度对数算法研究成功

经过许多个月的艰辛努力，我自己独立开发的高精度算术库（BCMP）终于可以告一段落了。另外，由我自己独创的对数log(x)算法也已经完成。测试报告显示，这个算法的性能是非常令人振奋的，在100-10000位精度下，可比现在最快软件Mathematica[注9]快2.8到8倍。

BCMP和GMP[注1]类似，利用BCMP可以计算任意精度整数和浮点数的算术运算，包括加，减，乘，除，开平方等。和GMP类似，我开发这套大数库的主要目标为追求极致的性能。现在，这一目标已基本达到。在中等精度（1万位十进制）以下，BCMP在大多数情况下，显著领先于GMP。当精度更高时，bcmp则慢于GMP，这是因为，大数运算库的关键是乘法，而不同大小的操作数需要使用不同的算法，没有一种算法可以适应于任意情况。BCMP中实现了basecase乘法[注2]，karatsuba乘法[注3]，Toom-Cook[注4]系列乘法，而尚未实现适应于操作数很大的SSA[注5]算法，或者FFT[注6]算法和NTT[注7]算法。和GMP不同，BCMP提供了2套接口，第一套接口支持2进制运算，内部表示采用\(2^{30}\)进制，另一套接口支持10进制运算，内部表示采用\(10^9\)进制。BCMP计划同时支持32位和64位方式编译，也计划同时支持在windows下编译和Linux下编译。现在已实现32位模式下的用Visual Studio编译，今后将逐步完善其他功能。

BCMP主要使用C语言和汇编语言开发，为了实现极高的性能，使用了大量的汇编代码，对关键部分做了精心的优化，除了普通的汇编指令外，还大量使用SSE2，SSE4和少数AVX2指令。

新的任意精度对数函数算法是我独创的，不同于任何一个公开发表的算法。新算法基于质数分解和比特突发算法，大大加速了计算速度，从算法复杂度来说，它和目前最快的算法AGM算法[注8]有相同的复杂度，都是O(log(P)M(P)),p为精度，其中M(P)代表两个P位整数相乘所需要的基本运算的次数。但是我的算法有更小的常数因子，速度更快。下表列出了几个最著名数学软件和大数库与BCMP在计算log(x)的性能对比结果。从这个表中可以看出，当精度为100到10000位时，我的新算法是MPFR[注11]的3.7到8.6倍

测试平台：
处理器: Intel(R) Core(TM) i7-2600K CPU @3.40GH
操作系统: Windows7 32-bit
Maple[注10] 版本：17.0
Mathematica版本： 10.3.0 32-bit
Pari[注12]版本： 2.9.3
mpfr版本：3.1.6

位数	时间，单位为毫秒
位数	Maple	Mathematica	Pari	MPFR	BCMP_log1	BCMP_log2
20	0.029941	0.00247424	0.003114	0.004833	0.008654	0.000762
30	0.040354	0.00305579	0.003744	0.005667	0.008218	0.000975
40	0.043768	0.003999	0.004687	0.006456	0.010237	0.001043
50	0.045405	0.00464357	0.005267	0.006782	0.010515	0.001217
60	0.061336	0.00551004	0.007068	0.008939	0.012052	0.002432
70	0.063722	0.00660306	0.008129	0.009373	0.012467	0.002622
80	0.06523	0.00748533	0.008918	0.010451	0.014011	0.002847
90	0.067488	0.00873332	0.009401	0.010859	0.013989	0.002967
100	0.069835	0.00976343	0.011404	0.012145	0.017113	0.003194
125	0.072167	0.0121835	0.013434	0.014376	0.020085	0.003843
158	0.079578	0.0186991	0.020295	0.018568	0.026015	0.004825
199	0.084953	0.0257783	0.027464	0.023096	0.032686	0.005604
251	0.101058	0.0337484	0.039437	0.030076	0.041313	0.006338
316	0.111391	0.0496319	0.057234	0.03874	0.057398	0.008075
398	0.125592	0.0655544	0.090215	0.054713	0.076384	0.009971
501	0.140205	0.0873334	0.126219	0.076801	0.100851	0.012947
630	0.169863	0.143491	0.173471	0.1227	0.126775	0.025564
794	0.201257	0.213486	0.262401	0.192741	0.175562	0.022624
1000	0.26965	0.331134	0.378275	0.29729	0.260477	0.041692
1258	0.344599	0.4512	0.545905	0.425357	0.321757	0.056777
1584	0.437988	0.615793	0.836217	0.671066	0.504561	0.077545
1995	0.570668	0.871282	1.219222	1.042976	0.671209	0.133183
2511	0.807129	1.20226	1.723848	1.448282	0.983084	0.190615
3162	1.07266	1.68033	2.561828	2.257336	1.489376	0.273305
3981	1.47754	2.42274	3.66907	2.959093	2.155975	0.510987
5011	2.15234	3.27768	5.208427	3.788049	3.193957	0.727249
6309	2.96615	4.7273	7.723928	5.709551	4.507839	1.019573
7943	4.0625	6.69787	11.071721	7.894751	6.753858	1.748821
10000	6.33854	9.14068	15.950279	10.198301	9.649934	2.538205
17782	15.1146	21.4501	38.916218	23.862088	23.569465	6.929607
31622	42.25	56.5504	94.352581	68.980759	62.309453	18.861327
56234	681	117.001	207.879794	131.337302	146.628323	45.943156
100000	2012.67	260.002	587.256462	296.982964	364.771994	117.998387

性能对比折线图，以mpfr为基准，令mpfr的运行时间为1，值越小越好。

性能对比折线图2，去掉Maple的测试结果，仍以mpfr为基准，值越小越好。

说明：BCMP用两种方法实现了两个版本的log函数，第一个使用AGM算法，在上表中标记为BCMP_log1, 第二个使用新算法，标记为BCMP_log2

测试代码下载：
mathematica,maple,pari,mpfr,bcmp,time

测试结果下载：
test_result
更多的任意精度浮点软件性能比较，请参阅Comparison of multiple-precision floating-point software

1.GMP是一个用于任意精度算术的免费库，用于对有符号整数，有理数和浮点数进行操作。GMP的主要目标应用是密码应用和研究，互联网安全应用和计算机代数系统。GMP是一个非常优秀的大数库，已有26年的历史。GMP的作者曾宣称，GMP是地球上最快的大数库,原文“the fastest bignum library on the planet!”。因为GMP非常快，所以很多应用软件使用GMP，如计算机代数系统Mathematica和Maple中的整数运算就是调用GMP来实现的。
关于GMP更多的信息，请参阅https://gmplib.org 和 https://en.wikipedia.org/wiki/GNU_Multiple_Precision_Arithmetic_Library

2.basecase乘法，也叫长乘法(Long multiplication),请参阅 https://en.wikipedia.org/wiki/Multiplication_algorithm

3.Karatsuba乘法，也称Karatsuba算法，这个算法由Anatoly Karatsuba在1960年发现的，是一种比长乘法更快的算法，请参阅https://en.wikipedia.org/wiki/Karatsuba_algorithm。

4.Toom–Cook乘法，狭义的Toom-Cook算法指Toom-3算法，可以看做是Karatsuba乘法的推广。广义的Toom-Cook乘法包括Toom4，Toom5等一大类算法。这个算法最早由俄罗斯数学家Andrei Toom提出，而Stephen Cook(是美国和加拿大计算机科学和数学家，图灵奖的获得者）给出更清晰的描述，请参阅维基百科https://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication

5.SSA算法指Schönhage–Strassen algorithm ，他是德国数学家和计算机科学家Arnold Schönhage和德国数学家Volker Strassen在1971开发出来的一个算法。关于这个算法的更多信息，请参阅维基百科https://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm 和百度百科https://baike.baidu.com/item/SSA/18691142，原始论文请参阅《A GMP-based Implementation of Schönhage-Strassen’s Large Integer Multiplication Algorithm》

6.FFT指快速傅里叶变换（fast Fourier transform），是基于复数的一种线性变换，这个算法在电子技术领域得到极其广泛的应用，是20世纪10大算法之一。另外，这个算法可用来计算大数乘法。关于FFT的更多信息，请参阅 https://en.wikipedia.org/wiki/Fast_Fourier_transform。关于20世纪10大算法，请参阅http://science.dataguru.cn/article-9643-1.html

7.NTT指数论变换(Number-theoretic transform)，它是FFT推广到有限域上的一种变换算法，和FFT不同，他不使用复数域上的加，减和乘法，而是使用模算术。apfloat中的大数乘法用的就是NTT算法。感兴趣的读者可以参阅蒋增荣的《数论变换》或者孙琦的《快速数论变换》等书。

8.AGM指算法几何平均数，关于这个算法的描述请参阅澳大利亚数学家和计算机科学家Richard P. Brent《Fast Algorithms for High-Precision Computation of Elementary Functions》和的Borwein兄弟的《The arithmetic-geometricmean and fast computation of elementary functions》和。Borwein兄弟都是数学家，他们出生在苏格兰。

9.Mathematica 是一款功能强大的科学计算软件，和MATLAB、Maple 并称为三大数学软件。请参阅百度百科
https://baike.baidu.com/item/Mathematica

10.Maple和Mathematica类似，也是一款功能强大的科学计算软件。请参阅百度百科https://baike.baidu.com/item/maple/2306572

11.MPFR是基于GMP的一套2进制浮点库，请参阅维基百科https://en.wikipedia.org/wiki/GNU_MPFR

12.PARI/GP是一个计算机代数系统，请参阅维基百科https://en.wikipedia.org/wiki/PARI/GP

2017年12月25日2017年12月25日

AGM函数近似值的估计

AGM是 Arithmetic-Geometric Mean的缩写，意为算术几何平均数。其定义为,给定两个正实数a0和b0，我们定义一个迭代过程

\( a_{i+1}= \frac {a_i+b_i} {2}, \qquad b_{i+1}= \sqrt {a_i \cdot b_i } \)

前者求两个数的算术平均数，后者求2个数的几何平均数。这两个序列有同样的极限，其值表示为A(a0,b0)，在AGM迭代过程中，序列收敛地非常快，一旦\(a_i\) 和\(b_i\)大小相近时，每迭代一次，有效数字加倍。下表显示了这两个序列的收敛速度[1]

n	\(a_n\)	\(b_n\)
0	24	6
1	15	12
2	13.5	13.41640786499873817845
3	13.45820393249936908922	13.45813903099098487720
4	13.45817148174517698321	13.45817148170605385831
5	13.45817148172561542076	13.458171481725615420761

高斯发现，AGM函数可以用来求椭圆积分的值

\( A(1,x)= \frac{\pi} {2F(x)} \)

F(x)为椭圆积分函数

\( \displaystyle \int_0^{\pi \over 2} \frac{d\theta} { \sqrt {1-(1-x^2 ) \sin^2 \theta)}} ) \)

AGM函数不但可以快速的计算椭圆积分，也可用来计算对数函数ln(x).
[2]给出一个椭圆积分F(x)和对数函数ln(s)的关系式

\(
F(\frac{4}{s})=\ln(s)+ \frac{4\ln(s)-4}{s^2}+\frac{36\ln(s)-42}{s^4} + \frac{1200 \ln(s)-1480}{3s^6} + O(\frac{1}{s^8}) \)

如果s是足够的大，F(4/s)是非常好的ln(s)的近似值。为了计算ln(s)到p位2进制数。s必须大于\(2^{p \over 2 } \)。故为了计算ln(x),我们首先需要找到一个实数s和整数m，使得

\( s=x \cdot2^m \gt 2^{p \over 2} \)

[2]给出如下求ln(x)的公式

\( \displaystyle ln(x)=\frac{\pi}{2A({4 \over s},1)} -m \ln (2) \qquad (1)\)

我们将这个公式稍作一下变形

\( \displaystyle \ln(x)=\frac{{1 \over 8}\pi s}{A({s \over 4},1)} -m \ln (2)= \frac{{1 \over 2}\pi \cdot x \cdot 2^{m-2}}{A(x \cdot 2^{m-2},1)} – m \ln(2)
\qquad (2) \)

反过来，利用这个公式，我们可以求AGM(1,y)的近似值，这里要求y>>1。在公式(2)中，我们取x为1，这样不但可以消去方程中的x，也可以消去ln(x)，这样公式(2)变为

\( \displaystyle \frac {{1 \over 2} \pi \cdot 2^{m-2}}{A(2^{m-2},1)}=m \ln(2) \quad (3)\)

令 \(y=2^{m-2}\), 则

\( \displaystyle A(y,1)= \frac { {1 \over 2} \pi \cdot y} {(log2(y)+2) \ln 2} = y \frac{{1 \over 2} \pi / \ln(2)}{log2(y)+2} \approx y \frac{2.26618007091}{log2(y)+2} \quad (4) \)

下表是y取某些值时，使用公式(4)得到的近似值和真实值的对比。我们可以发现，当y>10,这个近似公式可以精确到3位以上有效数字。

y	A(y,1)	公式(4)
10	4.25040709493	4.25819370444
\(2^4\)	6.03865834042	6.04314685577
\(2^6\)	18.1285335082	18.1294405673
\(2^8\)	58.0140204356	58.0142098154
\(2^{16}\)	8250.90983997	8250.90984041
\(2^{24}\)	1462315.09787	1462315.09787
\(2^{32}\)	286269685.042	286269685.042

参考文献
[1] agm(24, 6) at WolframAlpha. http://www.wolframalpha.com/input/?i=agm%2824%2C+6%29
[2] Muller J M. Elementary Functions: Algorithms and Implementation[M]. Springer Science+Business Media New York, 2016,130-131

2017年12月24日2017年12月25日

批量分解素因数(二)

本文中的程序使用压缩的格式存储M以内(包含M)的所有奇数的分解式。

和前一个的程序不同，这个程序侧重于使用尽可能少的内存空间来保存M以内的所有奇数的分解式.
本文中的程序可计算65759以内的所有奇数的因子表，其因子表占用的内存空间不到200KB，32881*2+65535*2=196,832字节.
和上一篇文章一样，M以内的奇数的分解式需要使用两个表格来存储，地址表addrTable和因子表factorTable.
地址表使用WORD类型，addrTable[i]==0,表示(i*2+3)为素数,否则addrTable[i]表示(i*2+3)分解式的第一项的地址.
factorTable[0]不使用，factorTable的最后一项也不使用,故其存储的因子表的总数不能超过65534.
分解式的每一项，仅用一个WORD类型变量来存储，即每一项仅需2个字节。
当奇数x=2*i+3是素数时，其分解式无需存储在factorTable，这时，置addrTable[i]==0。
当奇数x=2*i+3是合数时，且其素因子p的指数为1时，在这种情况下，置bit15=1，用bit0-bit14存储p，这里p的最大取值是32767，故M的最大取值范围不超过3*32767=98,301.
当奇数x=2*i+3是合数时，且其素因子p的指数大于1，在这种情况下，置bit15=0，用bit0-bit8存储p，用bit9-bit14存储指数e，这里p的最大取值不超过511，故M的最大取值范围不超过511*511=261,121.
更多的细节，请参阅源代码。

下面给出源代码

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

typedef unsigned char BYTE;
typedef unsigned long DWORD;
typedef unsigned short WORD;

//筛出n以内(包括n）的所有素数，并统计n以内所有自然数的非重复素因子的个数
//当函数完成后，若pBuff[i]==0, 表明i是素数，
//若pBuff[i]>0, 表明i是合数，且i的非重复素因子的个数是pBuff[i]

void sift(BYTE *pBuff,DWORD n)
{
	DWORD p,m;
	memset(pBuff,0,sizeof(BYTE)*(n+1));
	p=2;
	while (p<=n)
	{
		for (m=p*2;m<=n;m+=p)	//p素数，m是p的倍数，将p的倍数（不包括p）的非重复素因子的个数增加1
			pBuff[m]++;
		
		p++;	//调到下一个数
		while ( pBuff[p] && p<=n) //略过合数
			p++;
	}
}


// 函数 print_expression_on_zipTable输出n以内所有自然数的分解式
// addrTable是地址表，addrTable[i]==0,表示(i*2+3)是素数，否则addrTable[i]指向奇数(i*2+3)的因子表的地址,
// addrTable[i+1]-addrTable[i]表示奇数(i*2+3)的分解式的项数
// factorTable[addrTable[i]]为奇数(i*2+3)的分解式的第一项
void print_expression_on_zipTable(WORD *addrTable,WORD *factorTable,DWORD max)
{
	DWORD i,x,j,t,c;
	WORD *tab;
	
	for (x=3;x<=max;x+=2)
	{
		i=(x-3)/2;
		
		if ( addrTable[i]==0)
		{
			printf("%u=%u\n",x,x);
			continue;
		}
		
		printf("%u=",x);

		for (t=i+1;addrTable[t]==0;)
			t++; //略过素数
		c=addrTable[t]-addrTable[i];
		tab=factorTable+addrTable[i];
		for(j=0;j<c;j++)
		{
			DWORD p,e;

			if (j>0)
				printf("*");
			
			if ( tab[j] & 0x8000)
			{	p=( tab[j] & 0x7fff);  e=1;}
			else
			{	p=( tab[j] & 0x1ff); e=( (tab[j]>>9) & 0x3f);	}
			printf("(%u^%u)",p,e);
		}
		printf("\n");
	}
}

void build_zipFactorTable(DWORD max,BYTE *pBuff, const WORD *addrTable,WORD *factorTable)
{
	DWORD p,m,t;
	WORD e,item;
	int addrTableSize;
	WORD *writePointTab=NULL;
	
	//-----------------------------------------------------
	// addrTable[i]存储奇数(i*2+3)的分解式，故i的最大值为（max-3)/2
	// 另外，需要额外增加一项来存储max+2分解式的地址，以方便计算max分解式的项数
	// 故需要addrTableSize=(max-3)/2+2
	
	addrTableSize=(max-3)/2+2;	
	writePointTab=(WORD *)malloc(sizeof(WORD)*addrTableSize);
	if ( writePointTab==NULL)
	{
		printf("Alloc memory faield\n");
		return ;
	}
	memcpy(writePointTab,addrTable,sizeof(WORD)*addrTableSize);
	
	p=3;
	while (p<=max)
	{
		for (m=p*3;m<=max;m+=(p*2))	//m为p的奇数倍
		{
			e=1;	
			// m可能含有多个素因子p，如120含有3个因子2，下面的代码求出素因子p的指数
			// m是p的倍数，故p一定整除m，故素因子p的指数至少为1
			t=m/p;	
			while (t%p==0) {t/=p; e++; }
			if ( e==1)
				item= (WORD)(0x8000 | p);
			else
				item= (WORD)((e << 9) | p);
			factorTable[writePointTab[(m-3)/2]]=item;
			writePointTab[(m-3)/2]++;
		}
		
		p+=2;	//前进到下一个奇数
		while ( pBuff[p]>0 && p<=max) 
			p+=2;	//略过合数，使得p总是素数
	}
	
	free(writePointTab); writePointTab=NULL;
}

// 函数batchDecompose批量分解n以内所有奇数为素数乘积的形式，并打印其分解式
void batch_odd_Decompose(DWORD n,const char *fileName)
{
	BYTE *buff=NULL;		//用来得到n以内的所有素数，和所有合数的素因子的个数
	WORD *addrTable=NULL;	//地址表addrTable[i]表示奇数(i*2+3)的分解式在factorTable中的偏移地址
	WORD *factorTable=NULL;
	WORD nextPos;
	DWORD i,x,max;
	DWORD total,addrTableSize;

	if ( (n&1)==0)
		n++;	//前进到下一个奇数

	buff=(BYTE *)malloc(sizeof(BYTE)*(n+1));
	if (buff==NULL )
	{ printf("Alloc memory failed\n");	goto free_memory;}

	sift(buff,n);

	for (total=0,x=3;x<=n;x+=2)
	{
		if ( buff[x]>0)
		{
			//buff[x]为数x的分解式的项数，total为所有奇合数的分解式的项数和
			total += buff[x];	 
			
			// 因子表factorTable的地址使用WORD来表示，故最大容量为65536，factorTable[0]和最后一个元素不使用，
			// 故total必须<=65534
			if (total<=65534)	
				max=x; 				  
			else
			{
				total -= buff[x];
				break;
			}
		}
	}
	printf("max=%u\n",max);

	//addrTableSize[i]存储奇数i*2+3的分解式，故i的最大值为（max-3)/2
	//另外，需要额外增加一项来存储 max+2分解式的地址，故需要addrTableSize=(max-3)/2+2
	addrTableSize=(max-3)/2+2;	
	addrTable	=(WORD *)malloc(sizeof(WORD)*(addrTableSize));
	if (addrTable==NULL )
	{ printf("Alloc memory failed\n");	goto free_memory;}

	nextPos=1;
	for (x=3;x<=max;x+=2)
	{
		i=(x-3)/2;
		
		if ( buff[x]==0)		//x是素数
			addrTable[i]=0;
		else
		{
			addrTable[i]=nextPos;	
			nextPos += buff[x];
		}
	}

	i=(x-3)/2;
	addrTable[i]=nextPos;	//增加额外的一项，方便计算max分解式的项数

	factorTable=(WORD *)malloc(sizeof(WORD)*(total+2));
	if (factorTable==NULL)
	{ printf("Alloc memory failed\n");	goto free_memory; }

	build_zipFactorTable(max,buff,addrTable,factorTable);		//将max以内的奇数的分解式存储到factorTable

	print_expression_on_zipTable(addrTable,factorTable,max);	//打印将n以内所有自然数的分解式

free_memory:
	if (buff!=NULL) {free(buff);buff=NULL;}
	if (addrTable!=NULL) {free(addrTable);addrTable=NULL;}
	if (factorTable!=NULL) {free(factorTable);factorTable=NULL;}
}


int main(int argc, char* argv[])
{
	batch_odd_Decompose(3*32767,NULL);
	return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

#include <stdlib.h>

#include <stdio.h>

#include <string.h>

typedef unsigned char BYTE;

typedef unsigned long DWORD;

typedef unsigned short WORD;

//筛出n以内(包括n）的所有素数，并统计n以内所有自然数的非重复素因子的个数

//当函数完成后，若pBuff[i]==0, 表明i是素数，

//若pBuff[i]>0, 表明i是合数，且i的非重复素因子的个数是pBuff[i]

void sift(BYTE *pBuff,DWORD n)

{

DWORD p,m;

memset(pBuff,0,sizeof(BYTE)*(n+1));

p=2;

while (p<=n)

{

for (m=p*2;m<=n;m+=p) //p素数，m是p的倍数，将p的倍数（不包括p）的非重复素因子的个数增加1

pBuff[m]++;

p++; //调到下一个数

while ( pBuff[p] && p<=n) //略过合数

p++;

}

// 函数 print_expression_on_zipTable输出n以内所有自然数的分解式

// addrTable是地址表，addrTable[i]==0,表示(i*2+3)是素数，否则addrTable[i]指向奇数(i*2+3)的因子表的地址,

// addrTable[i+1]-addrTable[i]表示奇数(i*2+3)的分解式的项数

// factorTable[addrTable[i]]为奇数(i*2+3)的分解式的第一项

void print_expression_on_zipTable(WORD *addrTable,WORD *factorTable,DWORD max)

{

DWORD i,x,j,t,c;

WORD *tab;

for (x=3;x<=max;x+=2)

{

i=(x-3)/2;

if ( addrTable[i]==0)

{

printf("%u=%u\n",x,x);

continue;

}

printf("%u=",x);

for (t=i+1;addrTable[t]==0;)

t++; //略过素数

c=addrTable[t]-addrTable[i];

tab=factorTable+addrTable[i];

for(j=0;j<c;j++)

{

DWORD p,e;

if (j>0)

printf("*");

if ( tab[j] & 0x8000)

{ p=( tab[j] & 0x7fff); e=1;}

else

{ p=( tab[j] & 0x1ff); e=( (tab[j]>>9) & 0x3f); }

printf("(%u^%u)",p,e);

}

printf("\n");

}

void build_zipFactorTable(DWORD max,BYTE *pBuff, const WORD *addrTable,WORD *factorTable)

{

DWORD p,m,t;

WORD e,item;

int addrTableSize;

WORD *writePointTab=NULL;

//-----------------------------------------------------

// addrTable[i]存储奇数(i*2+3)的分解式，故i的最大值为（max-3)/2

// 另外，需要额外增加一项来存储max+2分解式的地址，以方便计算max分解式的项数

// 故需要addrTableSize=(max-3)/2+2

addrTableSize=(max-3)/2+2;

writePointTab=(WORD *)malloc(sizeof(WORD)*addrTableSize);

if ( writePointTab==NULL)

{

printf("Alloc memory faield\n");

return ;

}

memcpy(writePointTab,addrTable,sizeof(WORD)*addrTableSize);

p=3;

while (p<=max)

{

for (m=p*3;m<=max;m+=(p*2)) //m为p的奇数倍

{

e=1;

// m可能含有多个素因子p，如120含有3个因子2，下面的代码求出素因子p的指数

// m是p的倍数，故p一定整除m，故素因子p的指数至少为1

t=m/p;

while (t%p==0) {t/=p; e++; }

if ( e==1)

item= (WORD)(0x8000 | p);

else

item= (WORD)((e << 9) | p);

factorTable[writePointTab[(m-3)/2]]=item;

writePointTab[(m-3)/2]++;

}

p+=2; //前进到下一个奇数

while ( pBuff[p]>0 && p<=max)

p+=2; //略过合数，使得p总是素数

}

free(writePointTab); writePointTab=NULL;

}

// 函数batchDecompose批量分解n以内所有奇数为素数乘积的形式，并打印其分解式

void batch_odd_Decompose(DWORD n,const char *fileName)

{

BYTE *buff=NULL; //用来得到n以内的所有素数，和所有合数的素因子的个数

WORD *addrTable=NULL; //地址表addrTable[i]表示奇数(i*2+3)的分解式在factorTable中的偏移地址

WORD *factorTable=NULL;

WORD nextPos;

DWORD i,x,max;

DWORD total,addrTableSize;

if ( (n&1)==0)

n++; //前进到下一个奇数

buff=(BYTE *)malloc(sizeof(BYTE)*(n+1));

if (buff==NULL )

{ printf("Alloc memory failed\n"); goto free_memory;}

sift(buff,n);

for (total=0,x=3;x<=n;x+=2)

{

if ( buff[x]>0)

{

//buff[x]为数x的分解式的项数，total为所有奇合数的分解式的项数和

total += buff[x];

// 因子表factorTable的地址使用WORD来表示，故最大容量为65536，factorTable[0]和最后一个元素不使用，

// 故total必须<=65534

if (total<=65534)

max=x;

else

{

total -= buff[x];

break;

}

printf("max=%u\n",max);

//addrTableSize[i]存储奇数i*2+3的分解式，故i的最大值为（max-3)/2

//另外，需要额外增加一项来存储 max+2分解式的地址，故需要addrTableSize=(max-3)/2+2

addrTableSize=(max-3)/2+2;

addrTable =(WORD *)malloc(sizeof(WORD)*(addrTableSize));

if (addrTable==NULL )

{ printf("Alloc memory failed\n"); goto free_memory;}

nextPos=1;

for (x=3;x<=max;x+=2)

{

i=(x-3)/2;

if ( buff[x]==0) //x是素数

addrTable[i]=0;

else

{

addrTable[i]=nextPos;

nextPos += buff[x];

}

i=(x-3)/2;

addrTable[i]=nextPos; //增加额外的一项，方便计算max分解式的项数

factorTable=(WORD *)malloc(sizeof(WORD)*(total+2));

if (factorTable==NULL)

{ printf("Alloc memory failed\n"); goto free_memory; }

build_zipFactorTable(max,buff,addrTable,factorTable); //将max以内的奇数的分解式存储到factorTable

print_expression_on_zipTable(addrTable,factorTable,max); //打印将n以内所有自然数的分解式

free_memory:

if (buff!=NULL) {free(buff);buff=NULL;}

if (addrTable!=NULL) {free(addrTable);addrTable=NULL;}

if (factorTable!=NULL) {free(factorTable);factorTable=NULL;}

}

int main(int argc, char* argv[])

{

batch_odd_Decompose(3*32767,NULL);

return 0;

}

2017年12月24日2017年12月25日

批量分解素因数(一)

这个程序的功能是采用高效的算法，将n以内的的所有自然数分解质因数，并存储起来，最后输出其分解式。

1.自然数的分解

根据算术基本定理,一个大于1的自然数能够唯一分解成几个素数乘积的形式。
如\(120=2 \times 2 \times 2 \times 3 \times 5\)。我们这里，将自然数的分解式写成\(n=(p_1^{e_1}) \times (p_2^{e_2}) \times \cdots \times (p_m^{e_m})\)的形式，这里\(p_1,p_2,p_m\)表示素数，\(p_1^{e_1}\), \(p_2^{e_2}\)表示素数的方幂。
如\(120=(2^3)\times(3^1)\times(5^1)\)

2. 非重复素因子的个数

若一个自然数表示为几个素数乘积的形式，则去掉重复因子后，其因子的个数叫做这个自然数的非重复素因子的个数。如120含有3个非重复素因子，他们分别是2,3,5。容易看出，不同的自然数，其非重复素因子个数可能不同的，如27仅有1个非重复素因子3，而210含有2,3,5,7四个非重复素因子。\(2^{32}\)以内的所有整数中，其非重复素因子个数最多不超过9个，这个因为前10个素数的乘积\(2\times 3 \times 5 \times 7 \times 11 \times 13 \times 17 \times 19 \times 23 \times 29=6,469,693,230>2^{32}\)

3.自然数分解式的表示

若 \(n=(p_1^{e_1}) \times (p_2^{e_2}) \times \cdots \times (p_m^{e_m})\)，我们说这个自然数的分解式包含m项，每项包含2个成员,素因子p和其指数e。自然数的分解式可顺序存储，也可采用链式存储，前者其分解式的所有项的地址是相邻的，各项顺序存放。后者其分解式的每项的地址是不相邻的，使用单链表来存储。链式存需要额外的空间来存储指针，不但内存占用率高，且运行效率也低，故此文中的程序采用顺序存储方式。

4.算法

4.1 首先，函数sift采用与筛法求素数类似的方法，求出n以内所有的素数，并统计出每个合数的非重复素因子的个数。
4.2 接着，程序算出所有n以内的所有自然数分解式的项数的之和，分配内存来存储其分解式，并确定每个自然数的分解式的的地址。方法是：若c_arr[i]表示自然数i的素因子的个数，a_arr[i]表示自然数i的的分解式的存储地址，则自然数i+1的的分解式的地址a_arr[i+1]=a_arr[i]+c_arr[i]
4.3 然后，函数buildfactorTable将n以内的所有自然数分解，并存储其分解式。
4.4 最后，函数print_factor输出每个自然数的分解式。

5.复杂度分析

5.1 空间复杂度
函数batchDecompose 在运行过程中动态分配了3个数组，buff,addrtable和factorTable，前两个数组大小分别为n+1和n+2。最后一个数组的大小取决于n以内每个数的分解式的项数的和。n以内平均每个数的非重复素因子的个数约为log(log(n))，故factorTable的大小约为n*log(log(n))，故总的空间复杂度为o(n*log(log(n)))
5.2 时间复杂度
函数sift的时间复杂度度亦为o(n*log(log(n))),计算addrTable的时间复杂度为o(n),函数buildfactorTable的时间复杂度稍大于o(n*log(log(n))),故总的时间复杂度稍大于o(n*log(log(n)))

下面是C语言源代码

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

typedef unsigned char BYTE;
typedef unsigned long DWORD;

typedef struct _pe_pair
{
	DWORD p;	//p是素数
	DWORD e;	//e是指数
}PE_PAIR;

//筛出n以内(包括n）的所有素数，并统计n以内所有自然数的非重复素因子的个数
//当函数完成后，若pBuff[i]==0, 表明i是素数，
//若pBuff[i]>0, 表明i是合数，且i的非重复素因子的个数是pBuff[i]

void sift(BYTE *pBuff,DWORD n)
{
	DWORD p,m;
	memset(pBuff,0,sizeof(BYTE)*(n+1));
	p=2;
	while (p<=n)
	{
		for (m=p*2;m<=n;m+=p)	//p素数，m是p的倍数，将p的倍数（不包括p）的非重复素因子的个数增加1
			pBuff[m]++;
		
		p++;	//调到下一个数
		while ( pBuff[p] && p<=n) //略过合数
			p++;
	}
}

// 函数 print_factor输出n以内所有自然数的分解式
// addrTable是地址表，addrTable[i]表示自然数i的因子表的地址
// addrTable[i+1]-addrTable[i]表示自然数i的分解式的项数，addrTable包含n+2个元素，前3个元素不使用
// factorTable[addrTable[i]] 为自然数i的分解式的第一项
void print_expression(DWORD *addrTable,PE_PAIR *factorTable,DWORD n)
{
	DWORD i,j,c;
	PE_PAIR *tab;
	
	for (i=2;i<=n;i++)
	{
		c=addrTable[i+1]-addrTable[i];
		tab=factorTable+addrTable[i];

		printf("%u=",i);
		for(j=0;j<c;j++)
		{
			if (j>0)
				printf("*");
			printf("(%u^%u)",tab[j].p,tab[j].e);
		}
		printf("\n");
	}
}


/*函数buildTable对n以内的自然数进行分解，并将每个自然数的分解式存储到factorTable
输入参数
  n,pBuff,addrTable
输出参数
 factorTable
*/

void buildfactorTable(DWORD n,BYTE *pBuff, DWORD *addrTable,PE_PAIR *factorTable)
{
	DWORD i,p,m,t;
	PE_PAIR item;
	//-----------------------------------------------------
	p=2;
	while (p<=n)
	{
		item.p=p; item.e=1;
		factorTable[addrTable[p]]=item;
		addrTable[p]++;

		for (m=p*2;m<=n;m+=p)	//m为p的倍数
		{
			//m是p的倍数，故p一定整除m，故素因子p的指数至少为1
			item.p=p; item.e=1;	
			
			// m可能含有多个素因子p，如120含有3个因子2，下面的代码求出素因子p的指数
			t=m/p;	
			while (t%p==0) {t/=p; item.e++; }
			
			factorTable[addrTable[m]]=item;
			addrTable[m]++;
		}
		
		p++;
		while ( pBuff[p]>0 && p<=n) 
			p++;	//略过合数，使得p总是素数
	}
	
	//当完成上面的代码，addrTable[i]的值实际上为自然数i+1的因子表的首地址
	//故需要恢复addrTable为原始值
	for (i=n;i>=2;i--)
		addrTable[i+1]=addrTable[i];
	addrTable[2]=0;
}

//函数batchDecompose批量分解n以内所有自然数的，并打印其分解式
void batchDecompose(DWORD n)
{
	BYTE *buff=NULL;		//用来得到n以内的所有素数，和所有合数的素因子的个数
	DWORD *addrTable=NULL;	//地址表，addrTable[i]表示自然数i的分解式在factorTable中的偏移地址
	PE_PAIR *factorTable=NULL;
	DWORD i,c,total;

	buff		=(BYTE *)malloc(sizeof(BYTE)*(n+1));
	addrTable	=(DWORD *)malloc(sizeof(DWORD)*(n+2));
	if (buff==NULL || addrTable==NULL)
	{ printf("Alloc memory failed\n");	goto free_memory;}

	sift(buff,n);

	addrTable[0]=addrTable[1]=addrTable[2]=0;
	for (total=0,i=2;i<=n;i++)
	{
		c= (buff[i]==0 ? 1: buff[i]);	//c为自然数i的素因子的个数
		total += c;					  
		addrTable[i+1]=addrTable[i]+c;	//计算自然数i+1的分解式在factorTable中的偏移地址
	}

	factorTable=(PE_PAIR *)malloc(sizeof(PE_PAIR)*total);
	if (factorTable==NULL)
	{ printf("Alloc memory failed\n");	goto free_memory; }

	buildfactorTable(n,buff,addrTable,factorTable);	//将n以内的自然数的分解式存储到factorTable

	print_expression(addrTable,factorTable,n);	//打印将n以内所有自然数的分解式

free_memory:
	if (buff!=NULL) {free(buff);buff=NULL;}
	if (addrTable!=NULL) {free(addrTable);addrTable=NULL;}
	if (factorTable!=NULL) {free(factorTable);factorTable=NULL;}
}

int main(int argc, char* argv[])
{
	batchDecompose(100);
	return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

#include <stdlib.h>

#include <stdio.h>

#include <string.h>

typedef unsigned char BYTE;

typedef unsigned long DWORD;

typedef struct _pe_pair

{

DWORD p; //p是素数

DWORD e; //e是指数

}PE_PAIR;

//筛出n以内(包括n）的所有素数，并统计n以内所有自然数的非重复素因子的个数

//当函数完成后，若pBuff[i]==0, 表明i是素数，

//若pBuff[i]>0, 表明i是合数，且i的非重复素因子的个数是pBuff[i]

void sift(BYTE *pBuff,DWORD n)

{

DWORD p,m;

memset(pBuff,0,sizeof(BYTE)*(n+1));

p=2;

while (p<=n)

{

for (m=p*2;m<=n;m+=p) //p素数，m是p的倍数，将p的倍数（不包括p）的非重复素因子的个数增加1

pBuff[m]++;

p++; //调到下一个数

while ( pBuff[p] && p<=n) //略过合数

p++;

}

// 函数 print_factor输出n以内所有自然数的分解式

// addrTable是地址表，addrTable[i]表示自然数i的因子表的地址

// addrTable[i+1]-addrTable[i]表示自然数i的分解式的项数，addrTable包含n+2个元素，前3个元素不使用

// factorTable[addrTable[i]] 为自然数i的分解式的第一项

void print_expression(DWORD *addrTable,PE_PAIR *factorTable,DWORD n)

{

DWORD i,j,c;

PE_PAIR *tab;

for (i=2;i<=n;i++)

{

c=addrTable[i+1]-addrTable[i];

tab=factorTable+addrTable[i];

printf("%u=",i);

for(j=0;j<c;j++)

{

if (j>0)

printf("*");

printf("(%u^%u)",tab[j].p,tab[j].e);

}

printf("\n");

}

/*函数buildTable对n以内的自然数进行分解，并将每个自然数的分解式存储到factorTable

输入参数

n,pBuff,addrTable

输出参数

factorTable

void buildfactorTable(DWORD n,BYTE *pBuff, DWORD *addrTable,PE_PAIR *factorTable)

{

DWORD i,p,m,t;

PE_PAIR item;

//-----------------------------------------------------

p=2;

while (p<=n)

{

item.p=p; item.e=1;

factorTable[addrTable[p]]=item;

addrTable[p]++;

for (m=p*2;m<=n;m+=p) //m为p的倍数

{

//m是p的倍数，故p一定整除m，故素因子p的指数至少为1

item.p=p; item.e=1;

// m可能含有多个素因子p，如120含有3个因子2，下面的代码求出素因子p的指数

t=m/p;

while (t%p==0) {t/=p; item.e++; }

factorTable[addrTable[m]]=item;

addrTable[m]++;

}

p++;

while ( pBuff[p]>0 && p<=n)

p++; //略过合数，使得p总是素数

}

//当完成上面的代码，addrTable[i]的值实际上为自然数i+1的因子表的首地址

//故需要恢复addrTable为原始值

for (i=n;i>=2;i--)

addrTable[i+1]=addrTable[i];

addrTable[2]=0;

}

//函数batchDecompose批量分解n以内所有自然数的，并打印其分解式

void batchDecompose(DWORD n)

{

BYTE *buff=NULL; //用来得到n以内的所有素数，和所有合数的素因子的个数

DWORD *addrTable=NULL; //地址表，addrTable[i]表示自然数i的分解式在factorTable中的偏移地址

PE_PAIR *factorTable=NULL;

DWORD i,c,total;

buff =(BYTE *)malloc(sizeof(BYTE)*(n+1));

addrTable =(DWORD *)malloc(sizeof(DWORD)*(n+2));

if (buff==NULL || addrTable==NULL)

{ printf("Alloc memory failed\n"); goto free_memory;}

sift(buff,n);

addrTable[0]=addrTable[1]=addrTable[2]=0;

for (total=0,i=2;i<=n;i++)

{

c= (buff[i]==0 ? 1: buff[i]); //c为自然数i的素因子的个数

total += c;

addrTable[i+1]=addrTable[i]+c; //计算自然数i+1的分解式在factorTable中的偏移地址

}

factorTable=(PE_PAIR *)malloc(sizeof(PE_PAIR)*total);

if (factorTable==NULL)

{ printf("Alloc memory failed\n"); goto free_memory; }

buildfactorTable(n,buff,addrTable,factorTable); //将n以内的自然数的分解式存储到factorTable

print_expression(addrTable,factorTable,n); //打印将n以内所有自然数的分解式

free_memory:

if (buff!=NULL) {free(buff);buff=NULL;}

if (addrTable!=NULL) {free(addrTable);addrTable=NULL;}

if (factorTable!=NULL) {free(factorTable);factorTable=NULL;}

}

int main(int argc, char* argv[])

{

batchDecompose(100);

return 0;

}

2024年11月
一	二	三	四	五	六	日
« 1月
				1	2	3
4	5	6	7	8	9	10
11	12	13	14	15	16	17
18	19	20	21	22	23	24
25	26	27	28	29	30