【Miracl用户手册学习】——【第一篇：因子分解算法】

一、分解因子算法——Pollard rho算法

Pollard-Rho算法是John Pollard发明的一种能快速找到大整数的一个非1、非自身的因子的算法。
参考1
参考2
参考3

import math
from random import randint

def is_prime(number):
   if number > 1:
       if number == 2:
           return True
       if number % 2 == 0:
           return False
       for current in range(3, int(math.sqrt(number) + 1), 2):
           if number % current == 0:
               return False
       return True
   return False

def f(x,c,N):
   return (x*x +c)%N;

def Pollard_Rho(N):
   if N==4:
       return 2;
   if (is_prime(N)):
       return N;
   while 1:
       c = randint(1, N - 1);
       t = f(0,c,N)
       r = f(f(0,c,N),c,N)
       while (t != r):
           d = math.gcd(abs(t - r), N)
           if d>1:
               return d;
           t = f(t,c,N)
           r = f(f(r,c,N),c,N)
print(Pollard_Rho(7171))

import numpy as np
import matplotlib.pyplot as plt
x=0
c=24
N=9400
num = []
num.append(x)
for i in range(130):
   x=(x*x+c)%N
   num.append(x)

plt.scatter(range(131), num)
plt.plot([0,9.5],[9.5,0],linewidth = '0.5',color='#000000')
plt.legend()
plt.show()

二、brute force to factor numbers

遵循原理：任意整数可以唯一标识成素数的积

#include <iostream>
#include<big.h>
using namespace std;
extern "C"
{
#include "miracl.h"
#include "mirdef.h"

}
#define LIMIT 15000

int main()
{ /* find factors by brute force division */
	int n, p;

#ifdef MR_NOFULLWIDTH
	miracl *mip = mirsys(50, MAXBASE);
#else  
	miracl *mip = mirsys(50, 0);
#endif

	char stack_mem[mr_big_reserve(2, 50)];  // 2 bigs, 4 bytes per big...
	big x, y;

	memset(stack_mem, 0, mr_big_reserve(2, 50));
	x = mirvar_mem(stack_mem, 0);
	y = mirvar_mem(stack_mem, 1);

	/* This is an example of allocating bigs x and y from the stack rather than
	   from the heap. The "50" should be the same as that used in the prior call
	   to mirsys(.), and MUST be a constant. There is never a need to explicitly
	   delete bigs allocated from the stack, via a call to mirkill(). Therefore
	   do not use for sensitive data as these bigs are not erased, unless you
	   do it yourself via a call to zero(.)

	   Allocation from the stack is not particularly recommended, but could be
	   faster if bigs are not too big, and are required as temporary variables
	   in an often called subroutine. Allocating from the stack also avoids the
	   possibility of memory leaks.

	   NOTE: A big allocated from the stack has a random initial value

	*/
	gprime(LIMIT);
	printf("input number to be factored\n");
	cinnum(x, stdin);
	if (isprime(x))
	{
		printf("this number is prime!\n");
		return 0;
	}
	n = 0;
	p = mip->PRIMES[0];
	printf("factors are \n");
	forever //#define forever for(;;)   
	{ /* try division by each prime in turn */
		if (subdiv(x,p,y) == 0)
		{ /* factor found */
			copy(y,x);
			printf("prime factor     ");
			printf("%d\n",p);
			if (size(x) == 1) return 0;
			continue;
		}
		if (size(y) <= p)
		{ /* must be prime */
			printf("prime factor     ");
			cotnum(x,stdout);
			return 0;
		}
		p = mip->PRIMES[++n];
		if (p == 0) break;
	}
	if (isprime(x)) printf("prime factor     ");
	else            printf("composite factor ");
	cotnum(x, stdout);
	return 0;
}

gprime

isprime

size*

subdiv

三、brute force to factor numbers

#include <iostream>
#include<big.h>
using namespace std;
extern "C"
{
#include "miracl.h"
#include "mirdef.h"

}
/*
 *   Program to factor big numbers using Brent-Pollard method.
 *   See "An Improved Monte Carlo Factorization Algorithm"
 *   by Richard Brent in BIT Vol. 20 1980 pp 176-184
 *
 */

#include <stdio.h>
#ifdef MR_STATIC
#include <string.h>
#endif

#define mr_min(a,b) ((a) < (b)? (a) : (b))

int main()
{  /*  factoring program using Brents method */
	long k, r, i, m, iter;
	big x, y, z, n, q, ys, c3;
#ifndef MR_STATIC
	miracl *mip = mirsys(16, 0);         /* if allocating from the heap, specify size of bigs here */
	char *mem = (char *)memalloc(7);            /* allocate and clear memory from the heap for 7 bigs     */
#else
	miracl *mip = mirsys(MR_STATIC, 0);  /* If allocating from the stack, size of bigs is pre-defined */
	char mem[MR_BIG_RESERVE(7)];      /* reserve space on the stack for 7 bigs ...  */
	memset(mem, 0, MR_BIG_RESERVE(7));  /* ... and clear this memory */
#endif

	x = mirvar_mem(mem, 0);              /* 7 bigs have index from 0-6  */
	y = mirvar_mem(mem, 1);
	ys = mirvar_mem(mem, 2);
	z = mirvar_mem(mem, 3);
	n = mirvar_mem(mem, 4);
	q = mirvar_mem(mem, 5);
	c3 = mirvar_mem(mem, 6);
	convert(3, c3);

	printf("input number to be factored\n");
	cinnum(n, stdin);
	if (isprime(n))
	{
		printf("this number is prime!\n");
		return 0;
	}
	m = 10L;
	r = 1L;
	iter = 0L;
	do
	{
		printf("iterations=%5ld", iter);
		convert(1, q);
		do
		{
			copy(y, x);
			for (i = 1L; i <= r; i++)
				mad(y, y, c3, n, n, y);
			k = 0;
			do
			{
				iter++;
				if (iter % 10 == 0) printf("\b\b\b\b\b%5ld", iter);
				fflush(stdout);
				copy(y, ys);
				for (i = 1L; i <= mr_min(m, r - k); i++)
				{
					mad(y, y, c3, n, n, y);
					subtract(y, x, z);
					mad(z, q, q, n, n, q);
				}
				egcd(q, n, z);
				k += m;
			} while (k < r && size(z) == 1);
			r *= 2;
		} while (size(z) == 1);
		if (mr_compare(z, n) == 0) do
		{ /* back-track */
			mad(ys, ys, c3, n, n, ys);
			subtract(ys, x, z);
		} while (egcd(z, n, z) == 1);
		if (!isprime(z))
			printf("\ncomposite factor ");
		else printf("\nprime factor     ");
		cotnum(z, stdout);
		if (mr_compare(z, n) == 0) return 0;
		divide(n, z, n);
		divide(y, n, n);
	} while (!isprime(n));
	printf("prime factor     ");
	cotnum(n, stdout);
#ifndef MR_STATIC
	memkill(mem, 7);                  /* delete all 7 bigs */
#else
	memset(mem, 0, MR_BIG_RESERVE(7)); /* clear memory used for bigs */
#endif
	return 0;
}

memset(b, 0, sizeof(b))

memset是个函数，它在string.h头文件中有声明。它有三个参数，一是所要set的首地址，二是set的值，三是set的字节数。

b大部分情况下是个数组，因为数组在传参时能自动转成指向数组的首元素的指针；如果不是数组，应该写成memset(&b, 0, sizeof(b));