C语言 – Baocheng's home page

2018年1月6日

编译器中和64位编程有关的预定义宏

本文对分别测试VC，MinGW，GCC 三种编译器，32位和64位模式，共6种情况下，和64位编程有关的预定义宏的值。对跨平台编程具有参考意义。

Agner Fog 在他的《Calling conventions for different C++ compilers and operating systems》提到一些预定义宏。这里摘录如下。

注：下面的内容来自《Calling conventions for different C++ compilers and operating systems》,Last updated 2012-02-29，作者：By Agner Fog. Copenhagen University College .如与原文内容不符，以原文为准。

Most C++ compilers have a predefined macro containing the version number of the compiler. Programmers can use preprocessing directives to check for the existence of these macros in order to detect which compiler the program is compiled on and thereby fix problems with incompatible compilers.

表格23. 编译器版本预定义宏

Compiler	Predefined macro
Borland	__BORLANDC__
Codeplay VectorC	__VECTORC__
Digital Mars	__DMC__
Gnu	__GNUC__
Intel	__INTEL_COMPILER
Microsoft	_MSC_VER
Pathscale	__PATHSCALE__
Symantec	__SYMANTECC__
Watcom	__WATCOMC__

不幸的是，并不是所有编译器都有良好的文档来告诉你哪些宏是用来说明正在编译的硬件平台和操作系统的。以下宏可以被定义也可以不被定义

表 24. 硬件平台预定义宏

platform	Predefined macro
x86	_M_IX86	__INTEL__	__i386__
x86-64	_M_X64	__x86_64__	__amd64
IA64	__IA64__
DEC Alpha	__ALPHA__
Motorola Power PC	__POWERPC__
Any little endian	__LITTLE_ENDIAN__
Any big endian	__BIG_ENDIAN__

Table 25. 操作系统预定义宏

Operating system	Predefined macro
DOS 16 bit	__MSDOS__	_MSDOS
Windows 16 bit	_WIN16
Windows 32 bit	_WIN32	__WINDOWS__
Windows 64 bit	_WIN64	_WIN32
Linux 32 bit	__unix__	__linux__
Linux 64 bit	__unix__	__linux__	__LP64__	__amd64
BSD	__unix__	__BSD__	__FREEBSD__
Mac OS	__APPLE__	__DARWIN__	__MACH__
OS/2	__OS2__

下面的代码主要测试和64编程有关的宏

void test()
{
	int len=sizeof(int)*8;
	printf("sizeof(int)=%d\n",len);

	len=sizeof(int *)*8;
	printf("sizeof(int*)=%d\n",len);

#ifdef _MSC_VER
	printf("_MSC_VER is defined\n");
#endif

#ifdef __GNUC__
	printf("__GNUC__ is defined\n");
#endif

#ifdef __INTEL__ 
	printf("__INTEL__  is defined\n");
#endif

#ifdef __i386__ 
	printf("__i386__  is defined\n");
#endif

#ifdef __x86_64__ 
	printf("__x86_64__  is defined\n");
#endif

#ifdef _WIN32 
	printf("_WIN32 is defined\n");
#endif

#ifdef _WIN64 
	printf("_WIN64 is defined\n");
#endif


#ifdef __linux__ 
	printf("__linux__ is defined\n");
#endif

#ifdef __LP64__ 
	printf("__LP64__ is defined\n");
#endif


#ifdef __amd64 
	printf("__amd64 is defined\n");
#endif
}

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

void test()

{

int len=sizeof(int)*8;

printf("sizeof(int)=%d\n",len);

len=sizeof(int *)*8;

printf("sizeof(int*)=%d\n",len);

#ifdef _MSC_VER

printf("_MSC_VER is defined\n");

#endif

#ifdef __GNUC__

printf("__GNUC__ is defined\n");

#endif

#ifdef __INTEL__

printf("__INTEL__ is defined\n");

#endif

#ifdef __i386__

printf("__i386__ is defined\n");

#endif

#ifdef __x86_64__

printf("__x86_64__ is defined\n");

#endif

#ifdef _WIN32

printf("_WIN32 is defined\n");

#endif

#ifdef _WIN64

printf("_WIN64 is defined\n");

#endif

#ifdef __linux__

printf("__linux__ is defined\n");

#endif

#ifdef __LP64__

printf("__LP64__ is defined\n");

#endif

#ifdef __amd64

printf("__amd64 is defined\n");

#endif

}

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

{

test();

return 0;

}

最后给出结果.

1. 在VC2010， 32位模式下编译，输出下面的信息
sizeof(int)=32
sizeof(int*)=32
_MSC_VER is defined
_WIN32 is defined

2. 在MinGW下编译输出下面的信息
sizeof(int)=32
sizeof(int*)=32
__GNUC__ is defined
__i386__ is defined
_WIN32 is defined

3. 在64位Fedora19, 使用gcc -m32 编译，输出下面的信息
sizeof(int)=32
sizeof(int*)=32
__GNUC__ is defined
__i386__ is defined
__linux__ is defined

4. 在VC2010， 64位模式下编译，输出下面的信息
sizeof(int)=32
sizeof(int*)=64
_MSC_VER is defined
_WIN32 is defined
_WIN64 is defined

5. 在MinGW64下编译输出下面的信息
sizeof(int)=32
sizeof(int*)=64
__GNUC__ is defined
__x86_64__ is defined
_WIN32 is defined
_WIN64 is defined
__amd64 is defined

6. 在64位Fedora19, 使用gcc -m64 编译，输出下面的信息
sizeof(int)=32
sizeof(int*)=64
__GNUC__ is defined
__x86_64__ is defined
__linux__ is defined
__LP64__ is defined
__amd64 is defined

注：在VC下直接用集成环境编译，在MinGW和Linux下直接使用 gcc来编译。

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日

最短的计算大数乘法的c程序

#include <stdio.h>
char s[99],t[99];
int m,n;
void r(int i,int c)
{
  int j=0,k=i;
  while(k)c+=s[j++]*t[k---1];
  if(i)r(i-1,c/10);
  printf("%d",c%10);
}
void main()
{
  gets(s);gets(t);
  while(s[n])s[n++]-=48;
  while(t[m])t[m++]-=48;
  r(m+n-1,0);
}

#include <stdio.h>

char s[99],t[99];

int m,n;

void r(int i,int c)

{

int j=0,k=i;

while(k)c+=s[j++]*t[k---1];

if(i)r(i-1,c/10);

printf("%d",c%10);

}

void main()

{

gets(s);gets(t);

while(s[n])s[n++]-=48;

while(t[m])t[m++]-=48;

r(m+n-1,0);

}

说明：
1.这个程序接收2个从键盘输入的整数，计算他们的乘积，并输出结果。输入的两个整数的总长度不能大于99.
2.这个程序没什么大用，只是用来玩玩儿而已。
3.这个程序的主要目标是,使用尽可能短的代码来实现大数乘法。上面的代码
可在VC下编译并运行. 在GCC下编译，可省略#include语句和void关键字，
去除回车和不必要的空格，总长度仅仅194个字节。
另外，程序刻意避免使用数组来存贮中间结果和最终结果。
为此，使用了递归函数，同时，递归的使用也简化了代码。
4.在实际工作中，千万不要写这样的程序，否则会被骂死。
5.不要用这个程序考你的学生和面试者，即使他宣称精通C语言。

此类最短程序的特点
1.经常使用全局变量，全局变量的优点是
1).自动初始化数组和单变量为0，可省去某些变量初始化语句。
2).数组初始化为0也使得逻辑更简单，可省去某些边界值的判断。
3).在子程序，直接使用全局变量可省去某些参数定义和参数传递语句。

2.在表达式，大量使用“++”或者“–”之类运算符，此类语句往往起到
一箭双雕的效果，可有效的缩短代码长度.但在工作中，我强烈反对使用
这类运算符。

3.在比较语句中，很少使用if(i>0)这类语句，而是使用“if(i)”这样的
写法，这种写法比”>0″少了2个字母。

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

汇编语言编程中遇到的一个跨文件函数调用问题

我的一个项目中需要C和汇编语言混合编程，使用VS2010开发环境，汇编语言部分写在一个独立的.asm文件中，在VS工程中，汇编源文件的“Item Type“ 定义为“Custom Build Tool“。编译链接成功. 但是调试时发现，C语言文件中函数调用汇编语言文件中的函数可正确工作，但一个汇编文件中的函数调用另一个汇编源文件中的函数时却不能工作，当执行到call指令时，不能正确定位到目标函数，百度了1个多小时，参考了csdn论坛和博客的多个帖子和文章，也没有找到正确答案。最后在微软的官网https://support.microsoft.com/zh-cn/kb/86816上找到正确答案。解决方法，在调用函数所在的文件中，将被调函数声明为绝对地址。如 EXTERN _funname:abs，我之前的声明是 EXTERN _funname:DWORD。

更多的信息：

1. 我的程序生成的是32位应用程序，故在汇编源文件开头，做如下定义

.686P
.model flat
OPTION DOTNAME

2. Custom Build tool 参数定义

Command Line: ml /Zi /coff /c %(FullPath)
Outputs: %(Filename).obj %(Outputs)ml.exe

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