[fpc-pascal]isprime()

[astlab at inwind.it]

Hi Anton :)

> Does the function return false for 0 and 1?

Yes

> How easy isit to upgrade for > 64 bit lengths?

Indeed, isprime() is part of a routine working with the 'nat' big
numbers of the unit bigint10. Its purpose is to perform a primality
check by a strong test to the base 2 and a Lucas test. But up to 32
bits it is enough a Miller-Rabin test to bases 2, 7 and 61, which
covers up to 4,759,123,141. To do this, isprime() must evaluate
some powers and squares modulo a cardinal involving 64 bit operands.
I write here the function, so it may be used and tested. If there
are problems let me know.

Franco Milani
-------------------------------------------------------------------------

CONST

arlc_ispv: array [1..18] of cardinal =
(61,59,53,47,43,41,37,31,29,23,19,17,13,11,7,5,3,2);

VAR

vloc_ispa,vloc_ispb,vloc_ispc,vloc_ispd,vloc_ispe,vloc_ispf: cardinal;

FUNCTION isprime (n: cardinal): boolean; assembler;
asm
	movl	8(%ebp),%ebx
	leal	arlc_ispv,%edi
	movl	$17,%ecx
	cmpl	$61,%ebx
	ja	.L5
.L6:    cmpl	(%edi,%ecx,4),%ebx
	je	.L7
	decl	%ecx
	jns	.L6
	jmp	.L8
.L5:    movl	%ebx,%eax
	xorl	%edx,%edx
	movl	(%edi,%ecx,4),%esi
	divl	%esi
	cmpl	$0,%edx
	je	.L8
	decl	%ecx
	jns	.L5
	decl	%ebx
	movl	%ebx,vloc_ispa
	xorl	%ecx,%ecx
	testl	$1,%ebx
	jnz	.L9
.L10:   shrl	$1,%ebx
	incl	%ecx
	testl	$1,%ebx
	jz	.L10
.L9:    movl	%ebx,vloc_ispb
	cmpl	$0,%ecx
	je	.L11
	decl	%ecx
.L11:   movl	%ecx,vloc_ispc
	movl	$2,vloc_ispd
	movl	$2,%eax
	jmp	.L12
.L13:   movl	$7,%eax
	jmp	.L12
.L14:   movl	$61,%eax
.L12:   movl	vloc_ispb,%ebx
	movl	8(%ebp),%ecx
	cmpl	$1,%ebx
	jne	.L15
	xorl	%edx,%edx
	divl	%ecx
	movl	%edx,%eax
	jmp	.L16
.L15:   movl	%eax,%edi
	bsrl	%ebx,%edx
	movl	%edx,vloc_ispe
	btl	$0,%ebx
	jnc	.L17
	movl	%eax,%esi
	jmp	.L18
.L17:   movl	$1,%esi
.L18:   movl	$1,vloc_ispf
.L19:   movl	%edi,%eax
	mull	%edi
	cmpl	%ecx,%edx
	jb	.L20
	movl	%eax,%ebx
	movl	%edx,%eax
	xorl	%edx,%edx
	divl	%ecx
	movl	%ebx,%eax
.L20:   divl	%ecx
	movl	%edx,%edi
	movl	vloc_ispb,%eax
	movl	vloc_ispf,%ebx
	btl	%ebx,%eax
	jnc	.L21
	movl	%esi,%eax
	mull	%edi
	cmpl	%ecx,%edx
	jb	.L22
	movl	%eax,%ebx
	movl	%edx,%eax
	xorl	%edx,%edx
	divl	%ecx
	movl	%ebx,%eax
.L22:   divl	%ecx
	movl	%edx,%esi
.L21:   incl	vloc_ispf
	decl	vloc_ispe
	jnz	.L19
	movl	%esi,%eax
.L16:   cmpl	$1,%eax
	je	.L23
	cmpl	vloc_ispa,%eax
	je	.L23
	cmpl	$0,vloc_ispc
	je	.L8
	movl	vloc_ispc,%esi
.L24:   mull	%eax
	cmpl	%ecx,%edx
	jb	.L25
	movl	%eax,%ebx
	movl	%edx,%eax
	xorl	%edx,%edx
	divl	%ecx
	movl	%ebx,%eax
.L25:   divl	%ecx
	movl	%edx,%eax
	cmpl	vloc_ispa,%eax
	je	.L23
	cmpl	$1,%eax
	je	.L8
.L26:	decl	%esi
	jnz	.L24
	jmp	.L8
.L23:	decl	vloc_ispd
	jz	.L14
	jns	.L13
.L7:	movl	$1,%eax
	jmp	.L27
.L8:	xorl	%eax,%eax
.L27:	end;