A Mersenne number is a number in the form of 2P-1.

If P is prime, the Mersenne number may be a Mersenne prime
(if P is not prime, the Mersenne number is also not prime).

In the search for Mersenne prime numbers it is advantageous to eliminate exponents by finding a small factor before starting a, potentially lengthy, Lucas-Lehmer test.

There are very efficient algorithms for determining if a number divides 2P-1 (or equivalently, if 2P mod (the number) = 1).
Some languages already have built-in implementations of this exponent-and-mod operation (called modPow or similar).

The following is how to implement this modPow yourself:

For example, let's compute 223 mod 47.
Convert the exponent 23 to binary, you get 10111. Starting with square = 1, repeatedly square it.
Remove the top bit of the exponent, and if it's 1 multiply square by the base of the exponentiation (2), then compute square modulo 47.
Use the result of the modulo from the last step as the initial value of square in the next step:

                  remove       optional   
      square      top bit   multiply by 2   mod 47
   ────────────   ───────   ─────────────   ────── 
   1*1 = 1        1  0111   1*2 = 2            2
   2*2 = 4        0   111      no              4
   4*4 = 16       1    11   16*2 = 32         32
   32*32 = 1024   1     1   1024*2 = 2048     27
   27*27 = 729    1         729*2 = 1458       1

Since 223 mod 47 = 1, 47 is a factor of 2P-1.
(To see this, subtract 1 from both sides: 223-1 = 0 mod 47.)
Since we've shown that 47 is a factor, 223-1 is not prime.
Further properties of Mersenne numbers allow us to refine the process even more.
Any factor q of 2P-1 must be of the form 2kP+1, k being a positive integer or zero. Furthermore, q must be 1 or 7 mod 8.
Finally any potential factor q must be prime.
As in other trial division algorithms, the algorithm stops when 2kP+1 > sqrt(N).

These primality tests only work on Mersenne numbers where P is prime. For example, M4=15 yields no factors using these techniques, but factors into 3 and 5, neither of which fit 2kP+1.



Task

Using the above method find a factor of 2929-1 (aka M929)



Related tasks

   count in factors
   prime decomposition
   factors of an integer
   Sieve of Eratosthenes
   primality by trial division
   trial factoring of a Mersenne number
   partition an integer X into N primes
   sequence of primes by Trial Division


See also

   Computers in 1948: 2127 - 1 
       (Note:   This video is no longer available because the YouTube account associated with this video has been terminated.)


P:	equ	929		; P for 2^P-1
	cpu	8086
	bits	16
	org	100h
section	.text
	mov	ax,P		; Is P prime?
	call	prime
	mov	dx,notprm
	jc	msg		; If not, say so and stop.
	xor	bp,bp		; Let BP hold k
test_k:	inc	bp		; k += 1
	mov	ax,P		; Calculate 2kP + 1
	mul	bp		; AX = kP
	shl	ax,1		; AX = 2kP
	inc	ax		; AX = 2kP + 1
	mov	dx,ovfl		; If AX overflows (16 bits), say so and stop
	jc 	msg
	mov	bx,ax		; What is 2^P mod (2kP+1)?
	mov	cx,P
	call	modpow
	dec	ax		; If it is 1, we're done
	jnz	test_k		; If not, increment K and try again
	mov	dx,factor	; If so, we found a factor
	call	msg
prbx:	mov	ax,10		; The factor is still in BX
	xchg	bx,ax		; Put factor in AX and divisor (10) in BX
	mov	di,number	; Generate ASCII representation of number
digit:	xor	dx,dx
	div	bx		; Divide current number by 10,
	add	dl,'0'		; add '0' to remainder,
	dec	di		; move pointer back,
	mov	[di],dl		; store digit,
	test	ax,ax		; and if there are more digits,
	jnz	digit		; find the next digit.
	mov	dx,di		; Finally, print the number.
	jmp	msg
	;;;	Calculate 2^CX mod BX
	;;;	Output: AX
	;;;	Destroyed: CX, DX
modpow:	shl	cx,1		; Shift CX left until top bit in high bit
	jnc	modpow		; Keep shifting while carry zero
	rcr	cx,1		; Put the top bit back into CX
	mov	ax,1		; Start with square = 1
.step:	mul	ax		; Square (result is 32-bit, goes in DX:AX)
	shl	cx,1		; Grab a bit from CX
	jnc	.nodbl		; If zero, don't multiply by two
	shl	ax,1		; Shift DX:AX left (mul by two)
	rcl	dx,1
.nodbl:	div	bx		; Divide DX:AX by BX (putting modulus in DX)
	mov	ax,dx		; Continue with modulus
	jcxz	.done		; When CX reaches 0, we're done
	jmp	.step		; Otherwise, do the next step
.done:	ret
	;;;	Is AX prime?
	;;;	Output: carry clear if prime, set if not prime.
	;;;	Destroyed: AX, BX, CX, DX, SI, DI, BP
prime:	mov	cx,[prcnt]	; See if AX is already in the list of primes
	mov	di,primes
	repne	scasw		; If so, return
	je	modpow.done	; Reuse the RET just above here (carry clear)
	mov	bp,ax		; Move AX out of the way
	mov	bx,[di-2]	; Start generating new primes
.sieve:	inc	bx		; BX = last prime + 2
	inc	bx
	cmp	bp,bx		; If BX higher than number to test,
	jb	modpow.done	; stop, number is not prime. (carry set)
	mov	cx,[prcnt]	; CX = amount of current primes
	mov	si,primes	; SI = start of primes
.try:	mov	ax,bx		; BX divisible by current prime? 
	xor	dx,dx
	div	word [si]
	test	dx,dx		; If so, BX is not prime.
	jz	.sieve
	inc	si
	inc	si
	loop	.try		; Otherwise, try next prime.
	mov	ax,bx		; If we get here, BX _is_ prime
	stosw			; So add it to the list
	inc	word [prcnt]	; We have one more prime
	cmp	ax,bp		; Is it the prime we are looking for?
	jne	.sieve		; If not, try next prime
	ret
	;;;	Print message in DX
msg:	mov	ah,9
	int	21h
	ret
section	.data
	db	"*****"		; Placeholder for number
number:	db	"$"
notprm:	db	"P is not prime.$"
ovfl:	db	"Range of factor exceeded (max 16 bits)."
factor:	db	"Found factor: $"
prcnt:	dw	2		; Amount of primes currently in list
primes:	dw	2, 3		; List of primes to be extended