Problema al generar primos en Delphi

[[nlsgarcia]]

Club Delphi,

Los siguientes programas son un Compendio de Cribas de Generación de Números Primos implementadas en Delphi:

Criba de Eratóstenes:

Código Delphi [-]
program GeneratorPrimeNumbers;

{
Cálculo de Números Primos por el Algoritmo: Criba de Eratóstenes (Versión 2).
}

uses
  Windows, SysUtils, Classes, Dialogs;

var
   Limit, RLimit : Integer;
   i, j: Integer;
   Numbers: TBits;
   F : TextFile;
   NumberPrime : Integer;
   TI, TF: TDateTime;
   Mnsj : String;

begin

   repeat
      try
         Limit := StrToInt(InputBox('Generador de Números Primos',
                                    'Número Primo Máximo a Calcular:', '1000'));
      except
         Limit := 0;
      end;
   until (Limit >= 2) and (Limit <= 2147483615);

   TI := Now;

   try

      Numbers := TBits.Create;
      Numbers.Size := Limit+1;
      RLimit := Trunc(Sqrt(Limit));
      NumberPrime := 0;

      for i := 2 to RLimit do
         if not Numbers[i] then
            for j := i to (Limit div i) do
               Numbers[i*j] := True;

      FileMode := fmOpenWrite;
      AssignFile(F, 'NumberPrime.txt');
      Rewrite(F);

      for i := 2 to Limit do
         if not Numbers[i] then
         begin
            // Writeln(F, Format('%.10d',[i]));
            Writeln(F, i);
            Inc(NumberPrime);
         end;

   finally

      Numbers.Free;
      CloseFile(F);

   end;

   TF := Now - TI;

   ThousandSeparator := '.';
   DecimalSeparator := ',';

   Mnsj := Format('Con el Número %s como Límite, Se Generaron %s Números Primos en %s',
                 [FormatFloat('#,###,###,###,##0',Limit),
                  FormatFloat('#,###,###,###,##0',NumberPrime),
                  FormatDateTime('hh:mm:ss:zzz', TF)
                 ]);

   MessageBox(0, PChar(Mnsj), 'Algoritmo: Criba de Eratóstenes', MB_OK + MB_ICONINFORMATION);

end.

Criba de Atkin:

Código Delphi [-]
program GeneratorPrimeNumbers;

{
Cálculo de Números Primos por el Algoritmo: Criba de Atkin Optimizada
}

uses
  Windows, SysUtils, Classes, Dialogs;

var
   Limit, RLimit : LongWord;
   i, j : Integer;
   n,k : LongWord;
   R1 : LongWord;
   R2 : LongWord;
   xStepsize : LongWord;
   y_limit : LongWord;
   s, min_y, yy : LongWord;

   Numbers : TBits;
   F : TextFile;
   NumberPrime : LongWord;
   TI, TF : TDateTime;
   Msg : String;

begin

   repeat
      try
         Limit := StrToInt(InputBox('Generador de Números Primos',
                                    'Número Primo Máximo a Calcular:', '1000'));
      except
         Limit := 0;
      end;
   until (Limit >= 2) and (Limit <= 2147483615);

   TI := Now;

   Numbers := TBits.Create;

   try

      Numbers.Size := Limit+1;
      RLimit := Trunc(Sqrt(Limit));
      NumberPrime := 0;

      // Inicio del cálculo de 3x^2 + y^2
      xStepsize := 3;
      y_limit := 0;
      n := 0;
      R1 := Trunc(Sqrt((Limit - 1) / 3));

      i := 0;
      while (i < 12 * R1) do
      begin
         xStepsize := xStepsize + i;
         y_limit := 12 * Trunc(Sqrt(Limit - xStepsize)) - 36;
         n := xStepsize + 16;
         j := -12;
         while(j < y_limit + 1) do
         begin
            n := n + j;
            Numbers[n] := not Numbers[n];
            Inc(j,72);
         end;

         n := xStepsize + 4;

         j := 12;
         while (j < y_limit + 1) do
         begin
            n := n + j;
            Numbers[n] := not Numbers[n];
            inc(j,72);
         end;
         inc(i,24);
      end;
      // Fin del cálculo de 3x^2 + y^2

      // Inicio del cálculo de 4x^2 + y^2
      xStepsize := 0;
      R1 := 8 * Trunc(Sqrt((Limit - 1) / 4)) + 4;

      i := 4;
      while (i < R1) do
      begin
         xStepsize := xStepsize + i;
         n := xStepsize + 1;

         if (xStepsize mod 3 <> 0) then
         begin
            R2 := 4 * Trunc(Sqrt(Limit - xStepsize)) - 3;
            j := 0;
            while (j < R2) do
            begin
               n := n + j;
               Numbers[n] := not Numbers[n];
               Inc(j,8);
            end;
         end
         else
         begin
            y_limit := 12 * Trunc(Sqrt(Limit - xStepsize)) - 36;
            n := xStepsize + 25;
            j := -24;
            while (j < y_limit + 1) do
            begin
               n := n + j;
               Numbers[n] := not Numbers[n];
               Inc(j,72);
            end;

            n := xStepsize + 1;

            j := 24;
            while (j < y_limit + 1) do
            begin
               n := n + j;
               Numbers[n] := not Numbers[n];
               inc(j,72);
            end;
         end;
         inc(i, 8);
      end;
      // Fin del cálculo de 4x^2 + y^2

      // Inicio del cálculo de 3x^2 - y^2
      xStepsize := 1;
      R1 := Trunc(Sqrt(Limit/2))+1;

      i := 3;
      while (i < R1) do
      begin
         xStepsize := xStepsize + (4 * i - 4);
         n := 3 * xStepsize;
         s := 4;
         if (n > Limit) then
         begin
            min_y := (Trunc(Sqrt(n - Limit)) shr 2) shl 2;
            yy := min_y * min_y;
            n := n - yy;
            s := 4 * min_y + 4;
         end
         else
            s := 4;

         j := s;
         while (j < 4 * i) do
         begin
            n := n - j;
            if (n <= Limit) and (n mod 12 = 11) then
               Numbers[n] := not Numbers[n];
            Inc(j,8);
         end;
         inc(i,2);
      end;

      xStepsize := 0;
      i := 2;
      while (i < R1) do
      begin
         xStepsize := xStepsize + (4 * i - 4);
         n := 3 * xStepsize;
         s := 0;
         if (n > Limit) then
         begin
            min_y := ((Trunc(Sqrt(n - Limit)) shr 2) shl 2) - 1;
            yy := min_y * min_y;
            n := n - yy;
            s := 4 * min_y + 4;
         end
         else
         begin
            n := n - 1;
            s := 0;
         end;

         j := s;
         while(j < 4 * i) do
         begin
            n := n - j;
            if (n <= Limit) and (n mod 12 = 11) then
               Numbers[n] := not Numbers[n];
            inc(j,8);
         end;
         inc(i,2);
      end;
      // Fin del cálculo de 3x^2 - y^2

      i := 5;
      while (i <= RLimit+1) do
      begin
         if Numbers[i] then
         begin
            k := i*i;
            while (k < Limit) do
            begin
               Numbers[k] := False;
               Inc(k,i*i);
            end;
         end;
         Inc(i,2);
      end;

      FileMode := fmOpenWrite;
      AssignFile(F, 'NumberPrime.txt');
      Rewrite(F);

      Numbers[2] := True;
      Numbers[3] := True;

      for i := 2 to Limit do
         if Numbers[i] then
         begin
            Writeln(F, i);
            Inc(NumberPrime);
         end;

   finally

      Numbers.Free;
      CloseFile(F);

   end;

   TF := Now - TI;

   ThousandSeparator := '.';
   DecimalSeparator := ',';

   Msg := Format('Con el Número %s como Límite, Se Generaron %s Números Primos en %s',
                 [FormatFloat('#,###,###,###,##0',Limit),
                  FormatFloat('#,###,###,###,##0',NumberPrime),
                  FormatDateTime('hh:mm:ss:zzz', TF)
                 ]);

   MessageBox(0, PChar(Msg), 'Algoritmo: Criba de Atkin', MB_OK + MB_ICONINFORMATION);

end.

Criba de Sundaram:

Código Delphi [-]
program GeneratorPrimeNumbers;

{
Cálculo de Números Primos por el Algoritmo: Criba de Sundaram
}

uses
  Windows, SysUtils, Classes, Dialogs;

var
   Limit, RLimit : LongWord;
   i, j, k : LongWord;
   Numbers: TBits;
   F : TextFile;
   NumberPrime : LongWord;
   TI, TF: TDateTime;
   Mnsj : String;
   maxVal : LongWord;
   denominator : LongWord;
   Prime : LongWord;

begin

   repeat
      try
         Limit := StrToInt(InputBox('Generador de Números Primos',
                                    'Número Primo Máximo a Calcular:', '1000'));
      except
         Limit := 0;
      end;
   until (Limit >= 2) and (Limit <= 2147483615);

   TI := Now;

   try

      Numbers := TBits.Create;
      k := Limit div 2;
      Numbers.Size := k + 1;
      NumberPrime := 0;

      maxVal := 0;
      denominator := 0;

      i := 1;
      while(i < k) do
      begin
         denominator := (i shl 1) + 1;
         maxVal := (k - i) div denominator;
         j := i;
         while(j <= maxVal) do
         begin
            Numbers[i + j * denominator] := True;
            inc(j);
         end;
         inc(i);
      end;

      FileMode := fmOpenWrite;
      AssignFile(F, 'NumberPrime.txt');
      Rewrite(F);
      Writeln(F, 2);
      Inc(NumberPrime);

      Prime := 0;
      i := 1;
      while(i < k) do
      begin
         if not Numbers[i] then
         begin
            Prime := (i shl 1) + 1;
            Writeln(F, Prime);
            Inc(NumberPrime);
         end;
         inc(i);
      end;

   finally

      Numbers.Free;
      CloseFile(F);

   end;

   TF := Now - TI;

   ThousandSeparator := '.';
   DecimalSeparator := ',';

   Mnsj := Format('Con el Número %s como Límite, Se Generaron %s Números Primos en %s',
                 [FormatFloat('#,###,###,###,##0',Limit),
                  FormatFloat('#,###,###,###,##0',NumberPrime),
                  FormatDateTime('hh:mm:ss:zzz', TF)
                 ]);

   MessageBox(0, PChar(Mnsj), 'Algoritmo: Sieve of Sundaram', MB_OK + MB_ICONINFORMATION);

end.

Todas las anteriores implementaciones permite calcular con el Número 2.147.483.615 como cota límite, un máximo de 105.097.563 Números Primos en tiempos variables, sobre una máquina con un Procesador Phenom II X6 1090T, 4 GB RAM, 3 TB HDD y Windows 7 Profesional x32, como se muestra en la siguiente imagen a continuación:



El tiempo indicado en la imagen para los diferentes algoritmos es el tiempo total de proceso desde que inicia el cálculo hasta que finaliza la generación del archivo de 1.12 GB con los 105.097.563 Números Primos, lo cual implica que la verificación de los números primos: se hace en un tiempo menor al indicado.

Espero sea útil

Nelson.

[Victor Luis]

Holas Nelson...

La busqueda con el algoritmo de la criba de Eratostenes veo que es muy rapido... 2 minutos 25 segundos es un buen tiempo para un rango de poco mas de 2.000 millones...

○ Mi metodo PRI-BASE lo haria en 41 minutos, aunque ya encontre la forma de archivarlos en la mitad del tiempo que lo hacia, pero como dices el tiempo de busqueda es mas corto y lo que demora es archivar millones de numeros primos.
○ Una consulta al respecto, este tiempo es constante en cada busqueda del limite de 2.147.483.615 y por qué usan esta cantidad y no un limite entero... como 2.000 o 5.000 millones ?
Pregunto porque como indique el tiempo de mi metodo es constante, se cuando terminara no importa si busco en numeros grandes de mas de 18 o 24 digitos...

◘ Queria comentarte que intente poner el codigo en Delphi y he avanzado poco, creo que profundice mas VB, no encuentro las funciones compatibles como: Erase para inicializar un array, los Exponentes como 10^6, redondear hasta ciertos decimales y la funcion Format para dar formato a los numeros con puntos de mil o con decimales y me surge la duda si podre abrir los archivos de primos como otros con los Types que declare en Delphi, caso contrario tendria que iniciar desde 0 la busqueda y lo que queria es comparar los tiempos de las busquedas actuales...

◘ Por otro lado estoy analizando un metodo de Factorizacion para evaluar si un numero grande es primo o no, es una idea que me surgio y si resulta podria implementarlo al programa y ver que depure mayor cantidad de numeros grandes en menor tiempo.


Bueno amigo, para finalizar te dire que pronto tendre mi equipo actualizado con mas GB en el disco, mayor RAM y procesador, con lo que espero realizar busquedas con rangos mas grandes y reducir el tiempo del proceso....

Suerte y Gracias por la informacion que compartes....

[[nlsgarcia]]

Victor Luis,

Empezado por Victor Luis

...este tiempo es constante en cada búsqueda del limite de 2.147.483.615 y por qué usan esta cantidad...

Empezado por Victor Luis

...no encuentro las funciones...Erase para inicializar un array...los Exponentes...redondear...archivos...

Empezado por Victor Luis

...estoy analizando un método de Factorización para evaluar si un numero grande es primo o no...

Te comento:

1- El tiempo de generación de números primos nunca es constante en ningún algoritmo, varia en función de la carga del computador a nivel de recursos, pero en general se mantiene con pocas variaciones. El limite de 2.147.483.615 esta fijado en función de la clase TBits, la cual tiene como tamaño máximo dicho límite, algo similar ocurre con los arrays en Delphi y C#.

2- En lo referente a las funciones equivalentes entre VB6 y Delphi 7 te sugiero revisar este link:

Delphi Basics : http://www.delphibasics.co.uk/

3- Los métodos de factorización y primalidad no son adecuados para ser incluidos en un generador de números primos, dado que los tiempos se incrementarían notablemente, como ejemplo de un método de primalidad incluyo el de Miller-Rabin realizado en C# por tener el tipo BigInteger el cual facilita el uso de números enteros de gran tamaño, adecuado para su uso en algoritmos de números primos:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;
using System.Security.Cryptography;
using System.Windows.Forms;

namespace MillerRabin
{
    class Program
    {

        static bool IsProbablePrime(BigInteger source)
        {
            int certainty = 10;
            if (source == 2 || source == 3)
                return true;
            if (source < 2 || source % 2 == 0)
                return false;

            BigInteger d = source - 1;
            int s = 0;

            while (d % 2 == 0)
            {
                d /= 2;
                s += 1;
            }

            RandomNumberGenerator rng = RandomNumberGenerator.Create();
            byte[] bytes = new byte[source.ToByteArray().LongLength];
            BigInteger a;

            for (int i = 0; i < certainty; i++)
            {
                do
                {
                    rng.GetBytes(bytes);
                    a = new BigInteger(bytes);
                }
                while (a < 2 || a >= source - 2);

                BigInteger x = BigInteger.ModPow(a, d, source);
                if (x == 1 || x == source - 1)
                    continue;

                for (int r = 1; r < s; r++)
                {
                    x = BigInteger.ModPow(x, 2, source);
                    if (x == 1)
                        return false;
                    if (x == source - 1)
                        break;
                }

                if (x != source - 1)
                    return false;
            }

            return true;
        }

        static void Main(string[] args)
        {
            
            BigInteger Number;
            String AuxNumber;

            do
            {
               AuxNumber = Microsoft.VisualBasic.Interaction.InputBox("Número a Verificar Primalidad", 
                                                                       "Test de Primalidad MillerRabin", "");
               BigInteger.TryParse(AuxNumber, out Number);
               if (IsProbablePrime(Number))
                   MessageBox.Show("Es un Número Primo", "Test de Primalidad MillerRabin", 
                                   MessageBoxButtons.OK, MessageBoxIcon.Information);
               else if (AuxNumber != "")
                   MessageBox.Show("Es un Número Compuesto", "Test de Primalidad MillerRabin", 
                                   MessageBoxButtons.OK, MessageBoxIcon.Information);
            } while (AuxNumber != "");
            
        }
    }
}

El código anterior en C#, permite verificar la primalidad de un número por medio del algoritmo Miller-Rabin aplicado 10 veces por cada Primo Probable, a fin de garantizar su confiabilidad. Estoy revisando librerías que me permitan manejar grandes números enteros para trasladar el código a Delphi.

Pregunto: ¿Podrías explicar como manejas números enteros de 18 y 24 dígitos en VB6?

Espero sea útil

Nelson.

[Victor Luis]

Hola Amigo NELSON...

Realmente tienes una preparacion amplia y me quedo atras frente a vos, yo soy mas empirico y autodidacta... tengo un hermano mayor que estudio Ingenieria de Sistemas; sabe mucho pero como me conoce y dice siempre estoy queriendo descubrir la polvora, me ha enseñado algunas cosas pero la mayor parte de lo que se es por cuenta propia, pues como digo Todo se puede hacer en programacion y los programas no se equivocan, sino el programador....


◘ Respecto a lo que me comentas, de que el tiempo de generacion de numeros primos nunca es constante... no se si te refieres a la cantidad de primos encontrados lo cual es variable; pero conforme se avance esta cantidad disminuye paulatinamente. Pero si te refieres a tiempo en que se encuentran nuevos numeros primos en un rango dado, en mi caso es constante, buscando en 1.000 millones tardaba de 06:20 a 06:50 y ahora tarda entre 05:10 a 05:30 y no se si sigo mal... para mi esto es constante, no importa si busco en numeros de 12 digitos o de mas de 100 digitos, esto sigue siendo igual, claro que si aumentan o disminuyen unos segundos y eso no es constante, estoy equivocado....


◘ Referente al metodo de factorizacion, he leido un poco sobre los test de primalidad, no para copiar su metodo sino para entender hasta donde llego su logica... Es interesante por ejemplo el teorema de Fermat, la presipitacion de afirmar su descubrimiento que luego Euler comprobo las fallas o errores existentes, me falta leer mas de esa epoca; pero me pase al algoritmo AKS, el ultimo y eficiente que determina en tiempo polinominal (cosa que no me queda claro del todo) pero en una pagina indican que para la determinacion de un numero de 600 digitos emplearon unas 36 computadoras y el proceso tardo meses... su formula esta en "chino" para mi; pero ahora que avanzo y tengo mas numeros primos, se me ocurrieron dos metodos que los estoy analizando primero mediante una factorizacion directa y luego mediante una factorizacion multiple con otro enfoque, pues al evaluar numeros grandes las funciones aritmeticas no podran ser aplicadas.
☼ Se que debo estar cometiendo errores en mi expresion; pero com dije no soy matematico, es una idea que tengo, he probado unos calculos pero me falta corregirlos y ver que sea practico y lo elemental, debe ser simple, que aunque parezca ilogico es como tambien es mi metodo; pero funciona hasta ahora y ahora pude reducir el tiempo de busqueda a una (1/3) tercera parte de lo que lo hacia, por ejemplo buscar en un rango de 1.000 millones tardaba entre 21-22 min y ahora lo hace entre 7-8 min, incluido el archivado de los nuevos numeros primos, espero sea menor este tiempo luego que actualice mi equipo con mas RAM, procesador y Gigas de disco duro.

◘ Sobre tu pregunta de como manejo numeros de 18-24 digitos, es lo que pense cuando vi que mi metodo daba buenos resultados... y luego que, usando una variable Double o Decimal solo tienes para un monton de digitos; pero me cuestione si en verdad los calculos se haran de forma correcta, donde vi que desde 10 billones los ultimos digitos terminan siempre en "0" y eso me creo duda, no se, puedo estar equivocado al respecto; pero por lo simple de mi metodo, encontre un modo simple tambien, que funciona y su evaluacion final sera al pasar el cuatrillon y quintillon, donde al mismo tiempo me de datos para el analisis de la factorizacion que pienso hacer, por ahora amigo Nelson solo puedo decirte que en todo el proceso usare variables "double" y los numeros de mis calculos no pasaran de 21 a 24 digitos, osea esta variable es suficiente para mi... me preocupe por un rato e hice una verificacion inicial y esta bien.. ahora tal vez no entiendo a que te refieres con lo de "manejar numeros"... es como se hace, poner un valor a una variable y realizas calculos con ella... bueno no comprendo bien tu pregunta al respecto; pero te dije como lo estoy haciendo pero para estar seguro debo pasar con buenos resultados el cuatrillon, osea primos de mas de 24 digitos.

► Sobre el metodo de factotizacion, es una idea donde solo hice unos calculos que podrian resultar, las ecuaciones matematicas estan fuera de mi alcance, por lo que intento un metodo simple de acuerdo a mi conocimiento...

► Sobre el codigo a publicar, no avance nada, me parece ilogico pues ya mejore el rendimiento con VB; pero me dices que seria mas eficiente en Delphi... puede ser y no te contradigo, solo que como lo veo yo el tiempo se reducira y quedara estable entre 5 y maximo 6 minutos en cada busqueda de un rango de 1.000 millones, es mayor al algoritmo de Eratostenes; pero con mi equipo actualizado espero disminuya uno o dos minutos mas.


Bueno... Gracias por el enlace, esta en ingles y lo traducire... del codigo fuente que publicaste poco se de Java y peor C+, estoy actualizandome en Delphi 7, veo lo util de su compilador para evitar errores en la ejecucion, en lo que le supera a VB....
☼ Evitare no ser tan extenso; pero trato de explicarme de la forma mas clara y simple que puedo..... OK