Factorización de Fermat en R y Maple

Los números primos son aquellos que no tienen divisores, son los “primeros” de su género. En contraste los números que no son primos son compuestos (los componen números primos). 
El problema de saber si un numero en concreto es primo o no, no es fácil de solventar. El mismo Fermat ideó un algoritmo para determinar si un número es primo o no. 
Me ha hecho gracia escribir un programa en lenguaje R que aplique este algoritmo. Para los interesados aquí este el código en basto. 

---

Código sin pulir del algoritmo de factorizaron de Fermat para R.

Raw code of Fermat’s factoritation method for R.

#============
#x^2-y^2=n
#x^2-n=y^2
#============
#INTRODUCIR EL NÚMERO A FACTORIZAR

n<-23711

fras<-c("Los factores de ")
fras1<-c(" son ")
fras2<-c(" y ")
p<-sqrt(n)
p3<-trunc(p)
s<-(p3)
fin<-(n+1)/2

for (i in p3:fin){
s=s+1
r<-((s^2)-n)
a<-trunc(sqrt(r))
b<-sqrt(r)
if (a==b){      
            break
        }
    }
#==============================
#y= raíz cuadrada de r. x es s.
#tenemos que X=(a+b)/2. y=(a-b)/2
#ya que ab=n. Y a=x+y, b=x-y
#==============================
af<-s+a
bf<-s-a

paste (fras,n,fras1,af,fras2,bf)

---

El código funciona bien, pero resulta que R es un software pensado para la estadística y según el número que queramos comprobar nos dice que el cálculo es demasiado grande. Pasa si utilizamos el número: 502.560.280.658.509.
Así pues escribí otro código para Maple (es mi primer código en Maple), que es un software especial para matemáticas. También lo adjunto por si alguien llega aquí por alguna consulta googlena.

---

Código del algoritmo de factorizaron de Fermat para Maple.

Code of Fermat’s factoritation method for Maple. 

fermatfactor:= proc()
local a, s, fin, n, p, p3, b, r, x, i, af, bf; #variables
printf(`Programa que usa el Algoritmo de Fermat para factorizar`);
printf(`   Introduzca el número a factorizar:`);
n := scanf(%f)[1];
p:=sqrt(n);
p3:=trunc(p);

s:=p3;
fin:=(n+1)/2;
a:=1;
b:=2;
#Bucle
for i from p3 to fin   
while not (a=b) do 
r:=i^2-n;
a:=trunc(sqrt(r));
b:=sqrt(r);od;

x:=i-1;
af:=x+a;
bf:=x-a; 

print(`los factores de`,n,`son`,af,bf);

end:
# NOTA: Hay que llamar a la función con fermatfactor()

---

Este programa factorizó el anterior número que esta por los 500 billones (con b y europeo –en estados unidos un billón solo son mil millones-) en tres cuartos de hora. Por un teorema (que no expongo), solo nos hace falta buscar números desde el primer número entero cuyo cuadrado sea mayor que el número que queremos factorizar (en este caso 22.417.857) hasta el número a factorizar más uno y entre dos, es decir 251.280.140.329.255.

En resumen, el programa lo que hace es probar una ecuación, número a numero desde 22.417.857 a 251.280.140.329.255 en el peor de los casos. Para este número se evaluaron los números hasta el  23.969.353, es decir, tan solo un poco más de millón y medio de comprobaciones.