Experimentos con Arduino: marzo 2015

Recordando algo de matemáticas podemos decir que un número primo es un número entero que sólo es divisible entre 1 y si mismo. Los números primos son importantes para el común de los mortales porque una gran parte de la seguridad en internet se basa en utilizar un algoritmo muy dificil de descifrar ya que depende de dos números primos muy grandes.
Para mortales no tan comunes , los matemáticos, los números primos son especialmente interesantes ya que no se ha encontrado un algoritmo eficaz para determinar si un número es o no primo. El único algoritmo totalmente determinístico es la criba de eratóstenes , pero este procedimiento requiere recursos inmensos de tiempo y poder de cálculo a medida que los números son mas grandes.
Para otros mortales, como los que tenemos algo que ver con informática, desde muy temprano hemos tenido que desarrollar un programa que emule la criba de eratóstenes. Inclusive se utilizan ciertos algoritmos para determinar el poder de cálculo de un procesador (benchmark).
Es por ello que siempre encontraremos a alguien gastando el tiempo en desarrollar algún programa de búsqueda de números primos para cualquier procesador que tenga entre manos. Y claro, el Arduino no va a quedar huerfano en esta tan encomiable labor !!!
Buscando en la web me encontré con la siguiente entrada en el foro de Arduino http://forum.arduino.cc/index.php?topic=192209.0 que trata sobre la búsqueda de números primos grandes en la plataforma Arduino.
Para quien tenga unos quince minutos le recomiendo que la revisen aunque sea por encima, ya que entre otras cosas podemos ver como evoluciona el programa desarrollado por el Sr. Nick Gammon (por cierto, muy prolífico en temas de electronica y específicamente Arduino).
Como a mi esto de la teoría de números también me atrae, quise hacer un pequeño aporte.
Haciéndole unas pequeñas modificaciones al programa del Sr. Gammon logré que el programa se ejecutara en un Atmega644 a 20Mhz con una LCD donde se muestran los números primos que van siendo identificados. El programa guarda en el eeprom interno la información cada 1000 números primos, de tal forma que continúa despues de un reset o apagado. Este número puede ser cambiado en el programa.
La primera que vez que se graba el programa en la memoria flash del micro debemos eliminar los caracteres de comentarios en las líneas 43,44 y 45 , una vez que el programa ha encontrado los primeros 1000 números primos comentamos de nuevo estas líneas y grabamos de nuevo el programa en el micro.

// 64 bits primes
PROGMEM prog_uchar primes[] = {
    3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 
   37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 
   83, 89, 97,101,103,107,109,113,127,131,137,
  139,149,151,157,163,167,173,179,181,191,193,
  197,199,211,223,227,229,233,239,241,251
};
const byte nPrimes = sizeof(primes)/sizeof(primes[0]);
#include <EEPROM.h>
#include <LiquidCrystal.h>
LiquidCrystal lcd(7, 6, 5, 4, 3, 2); //Configuración de pines del LCD

#include "Pseudoprimes.h";
const uint16_t nPSP = (sizeof(psps)/sizeof(psps[0])) >> 3;
  
uint64_t candidate;
uint64_t leido;
uint64_t Barr_m;
uint8_t Barr_k;

long start;
long timerStart, timerRuntime;
long found = 0; // Number of primes found
int led = 0;
int Direccion = 0;
int Informacion;
char prime_eee[16];
int  grabar_cada=1000; //Cada 1000 primos se graba la información en el eeeprom
int  grabar_van =0;
float  tiempo_min;
float  tiempo_min0;

void setup() 
{ 
  lcd.begin(16, 2);
   
  
  leido=EEPROMReadu64(0);
  found=EEPROMReadFloat(30);
  tiempo_min0=EEPROMReadInt(40);
  
  //leido =  1999999999999999ULL;
  //found=0;
  //tiempo_min0=0;
  
  lcd.setCursor(0, 0);
  printll(leido ); 
  candidate=leido;
  lcd.setCursor(0, 1);
  lcd.print((int)tiempo_min0 ); 
  lcd.setCursor(9, 1);
  lcd.print(found); 
  delay(5000);
   
  pinMode(led, OUTPUT);   
}  // end setup()

void loop() 
{
  uint64_t nextPSP;
  int nextPSPindex;
  for (int i=0;i<nPSP;i++) {
    nextPSP = fetchPSP(i);
    if ((nextPSP == candidate) || (nextPSP > candidate)) {
      nextPSPindex = i;
      break;
    }
    nextPSP = 0;
    nextPSPindex = 0;
  }
  start = millis();
  
  if ((candidate & 1) == 0) {
    candidate++;
  }
  
  for ( ; candidate < 9999999999999999ULL; candidate += 2) {
    if (nextPSP != candidate) {
      if (isPrime(candidate)) {
        found++;
        showPrime (candidate);
        timerRuntime= 0;
      }
    }
    else {
      nextPSPindex++;
      if (nextPSPindex < nPSP) {
        nextPSP = fetchPSP(nextPSPindex);
      }
    }
  }
  found = 1;
  showPrime(1);
  found++;
  showPrime(2);
  candidate = 3;  // back to start
}  // end loop()

uint64_t mulmodBarrett(uint64_t a, uint64_t b, uint64_t c) {
  // Calculate a*b %c

  uint64_t ah; 
  uint64_t al;
  uint64_t ph; 
  uint64_t pm;
  uint64_t pl;

  multll(a, b, &ph, &pl);
  ah = ph;
  al = pl;

  if (ah == 0) {
    if (al < c) {
      return al;
    }
    else {
      if (al == c) {
        return 0;
      }
      else {
        al %= c;
        return al;
      }
    }
  }
  
  mult4lx2l(ah, al, Barr_m, &ph, &pm, &pl);
  uint8_t Bk = Barr_k << 1;
  if (Bk > 64) {
    if (Bk > 128) {
      Bk -= 128;
      pl = ph >> Bk;
    }
    else {
      Bk -= 64;
      pl = (pm >> Bk) | (ph << (64 - Bk));
    }
  }
  else {
    pl = (pl >> Bk) | (pm << (64 - Bk));
  }      
      
  multll(pl, c, &ph, &pm);
  pm = al - pm;
  if (pm > ah) {
    ph--;
  }
  ph = ah - ph;
  if (pm > c) {
    pm -= c;
  }
  return pm;
} // end mulmodBarrett()

void mult4lx2l(uint64_t ah, uint64_t al, uint64_t b, uint64_t *prodh, uint64_t *prodm, uint64_t *prodl) {
  // Multiply a and b
  multll(ah, b, prodh, prodm);
  uint64_t prodt;
  multll(al, b, &prodt, prodl);
  *prodm += prodt;
  if (*prodm < prodt) {
    (*prodh)++;
  }
}  

bool Miller(uint64_t source)
{
  if(source <= 3) {
    return true;
  }
  if((source & 1) == 0) {
    return false;
  }

  uint64_t d = source - 1;
  int s = 0;
  while((d & 1) == 0)
  {
    d >>= 1;
    s ++;
  }

  uint8_t a;
  // si se incluyen los demas primos se incrementa en 500ms 
  const byte bases[] = {2, 3, 5, 7 };// 11, 13, 17, 19, 23 };
  const byte nbases = sizeof(bases)/sizeof(bases[0]);
  for(int i = 0; i < nbases; i++)
  {
    a = bases[i];
    uint64_t x = powmodll(a, d, source);
   if((x == 1) || (x == (source - 1))) {
      continue;
    }
    bool stillOK = false;
    for(int r = 1; r < s; r++)
    {
      x = mulmodBarrett(x, x, source);
      if(x == 1)
        return false;
      if(x == (source - 1)) {
        stillOK = true;
        break;
      }
    }
    if(!stillOK)
      return false;
  }
  return true;
}  // end Miller()

bool isPrime(uint64_t n) {
  if (!check(n)) {
    return false;
  }
  Barr_k = Barrett_k(candidate);
  Barr_m = Barrett_m(candidate,Barr_k);
  if (n > 0x10000) {
    if (!Miller (n)) {
      return false;
    }
  }
  return true;
}

void showPrime (uint64_t n) 
{
  digitalWrite(led, HIGH);
  printll(n);
  grabar_van++;
  tiempo_min=tiempo_min0+(millis() - start)/60000.0;
  if (grabar_van==grabar_cada) {
    EEPROMWriteu64(0, n);  
    EEPROMWriteFloat(30, found);
    EEPROMWriteInt(40, (int)tiempo_min );
    lcd.setCursor(15, 1);
    lcd.print("*");
    grabar_van=0;
  }  else {
    lcd.setCursor(15, 1);
    lcd.print(" ");
  }
  lcd.setCursor(0, 1);
  lcd.print((int)tiempo_min );
  lcd.setCursor(4, 1);  lcd.print("m Van ");  
  lcd.setCursor(10, 1);
  lcd.print(found ); 
  digitalWrite(led, LOW);
  //start = millis();
}  // end showPrime()

bool check(uint64_t n) {
  if (n <= 3) {
    return true;
  }
  if (!(n & 1)) {
    return false;
  }
  uint8_t i = 1;
  uint8_t f = pgm_read_byte_near(primes);
  while ((i <= nPrimes) && (f < n)) {
    if ((n % f) == 0) {
      return false;
    }
    f = pgm_read_byte_near(primes+i);
    i++;
  }
  return true;
} // end check()

uint64_t fetchPSP(int i) {
  i <<= 3;
  uint64_t n = 0;
  for (int8_t k=0;k<sizeof(n);k++) {
    n = n <<= 8;
    n |= pgm_read_byte_near(psps+i+k);
  }
  return n;
} // end fetchPSP()

void printll(uint64_t n) {
  // print an unsigned long long
  if (n == 0) {  
    return;
  }
  const byte DIGITS = 21;
  char s[DIGITS+1];
  for (byte i=0;i<DIGITS;i++) {
    s[i] = 'X';
  }
  s[DIGITS] = 0;
  byte i = DIGITS-1;
  while (n > 0) {
    s[i] = (n % 10) | '0';
    n = n / 10;
    i--;
  }
  i++;
  char *z = &s[i];
  
  lcd.clear();
  lcd.setCursor(0, 0);
  lcd.print(z );
   
  // Palindrome
  
} // end printll()

uint8_t Barrett_k(uint64_t n) {
  uint8_t k = 0;
  uint8_t t;
  while (n > 0xFF) {
    n >>= 8;
    k += 8;
  }
  t = n;
  while (t > 0) {
    t >>= 1;
    k++;
  }
  return k;
}

uint64_t Barrett_m(uint64_t n, uint8_t k) {
  uint8_t t;
  uint64_t iml, imh;
  
  if (k<32) {
    t = 2 * k;
  }
  else {
    t = 2 *k - 64;
  }
  iml = 1;
  for (uint8_t i=8;i<=t;i+=8) {
    iml <<= 8;
  }
  t &= 7;
  for (uint8_t i=0;i<t;i++) {
    iml <<= 1;
  }
  imh = 0;
  if (k>=32) {
    imh = iml;
    iml = 0;
  }
  return quotll(n, imh, iml);
}
  
uint64_t powmodll(uint8_t a, uint64_t b, uint64_t c) {
  // Calculate a**b % c,
  uint64_t x = 1;
  uint64_t y = a;

  while(b > 0) {
    if(b & 1) {
        x = mulmodBarrett(x, y, c);
      }
      y = mulmodBarrett(y, y, c);
      b >>= 1;
    }
    return x;
} // end powmodll()

void multll(uint64_t a, uint64_t b, uint64_t *imh, uint64_t *iml) {
  // Multiply a and b

  *imh = (a >> 32) * (b >> 32);
  uint64_t imm = (a & 0xffffffffUL) * (b >> 32);
  *iml = (a >> 32) * (b & 0xffffffffUL);
  imm += *iml;
  if (imm < *iml) {     // Check for carry
    *imh += 1ULL << 32;
  }
  *iml = (a & 0xffffffffUL) * (b & 0xffffffffUL);
  imm += *iml >> 32;
  if (imm < (*iml >> 32)) {     // Check for carry
    *imh += 1ULL << 32;
  }
  *iml = (*iml & 0xffffffffUL) + (imm << 32);
  *imh += imm >> 32;
}

uint64_t quotll(uint64_t c, uint64_t imh, uint64_t iml) {
  uint64_t thisQuotient;
  
  if (imh == 0) {
    return (iml / c);
  }
  uint64_t quotient = 0;
  for (uint8_t i=0;i<8;i++) {
    if ((imh >> 56) != 0) {
      if (imh > c) {
        thisQuotient = (imh / c);
        quotient |= thisQuotient;
        imh = imh - (thisQuotient * c);
      }
      else {
        if (imh == c) {
          quotient |= 1;
          imh = 0;
        }
      }
    }
    imh = imh << 8;
    imh |= iml >> 56;
    iml = iml << 8;
    quotient = quotient << 8;
  }
  if (imh > c) {
    quotient |= (imh / c);
  }
  return quotient;
}

void printllHex(uint64_t n) {
  char t[19] = "0x0000000000000000";
  char s[] = "0123456789ABCDEF";

  for (int8_t i=17;i>1;i--) {
    t[i] = s[n & 0xf];
    n = n >> 4;
  }
} 

void startTimer() {
  timerStart = micros();
}

void accTimer() {
  timerRuntime += micros() - timerStart;
}

void EEPROMWriteu64(int address, uint64_t value)
       {
       byte eight = (value & 0xFF);
       byte seven = ((value >> 8)  & 0xFF);
       byte six   = ((value >> 16) & 0xFF);
       byte five  = ((value >> 24) & 0xFF);
       byte four  = ((value >> 32) & 0xFF);
       byte three = ((value >> 40) & 0xFF);
       byte two   = ((value >> 48) & 0xFF);
       byte one   = ((value >> 56) & 0xFF);

       EEPROM.write(address, eight);
       EEPROM.write(address + 1, seven);
       EEPROM.write(address + 2, six);
       EEPROM.write(address + 3, five);
       EEPROM.write(address + 4, four);
       EEPROM.write(address + 5, three);
       EEPROM.write(address + 6, two);
       EEPROM.write(address + 7, one);
       }

uint64_t EEPROMReadu64(long address)
       {
       uint64_t eight = EEPROM.read(address);
       uint64_t seven = EEPROM.read(address + 1);
       uint64_t six   = EEPROM.read(address + 2);
       uint64_t five  = EEPROM.read(address + 3);
       uint64_t four  = EEPROM.read(address + 4);
       uint64_t three = EEPROM.read(address + 5);
       uint64_t two   = EEPROM.read(address + 6);
       uint64_t one   = EEPROM.read(address + 7);

       return ((eight << 0)  & 0xFF) + 
              ((seven << 8)  & 0xFFFF) + 
              ((six << 16)   & 0xFFFFFF) + 
              ((five << 24)  & 0xFFFFFFFF) +
              ((four << 32)  & 0xFFFFFFFFFFULL) + 
              ((three << 40) & 0xFFFFFFFFFFFFULL) +
              ((two << 48)   & 0xFFFFFFFFFFFFFFULL) + 
              ((one << 56)   & 0xFFFFFFFFFFFFFFFFULL)
       ;
       }
       
void EEPROMWriteFloat(int p_address, long p_value)
{
  byte Byte1 = ((p_value >> 0) & 0xFF);
  byte Byte2 = ((p_value >> 8) & 0xFF);
  byte Byte3 = ((p_value >> 16) & 0xFF);
  byte Byte4 = ((p_value >> 24) & 0xFF);

  EEPROM.write(p_address, Byte1);
  EEPROM.write(p_address + 1, Byte2);
  EEPROM.write(p_address + 2, Byte3);
  EEPROM.write(p_address + 3, Byte4);
}       

float EEPROMReadFloat(int p_address)
{
  byte Byte1 = EEPROM.read(p_address);
  byte Byte2 = EEPROM.read(p_address + 1);
  byte Byte3 = EEPROM.read(p_address + 2);
  byte Byte4 = EEPROM.read(p_address + 3);

  long firstTwoBytes = ((Byte1 << 0) & 0xFF) + ((Byte2 << 8) & 0xFF00);
  long secondTwoBytes = (((Byte3 << 0) & 0xFF) + ((Byte4 << 8) & 0xFF00));
  secondTwoBytes *= 65536; // multiply by 2 to power 16 - bit shift 24 to the left

  return (firstTwoBytes + secondTwoBytes);
}

void EEPROMWriteInt(int p_address, int p_value)
      {
      byte lowByte = ((p_value >> 0) & 0xFF);
      byte highByte = ((p_value >> 8) & 0xFF);

      EEPROM.write(p_address, lowByte);
      EEPROM.write(p_address + 1, highByte);
      }


unsigned int EEPROMReadInt(int p_address)
      {
      byte lowByte = EEPROM.read(p_address);
      byte highByte = EEPROM.read(p_address + 1);
      return ((lowByte << 0) & 0xFF) + ((highByte << 8) & 0xFF00);
      }

En la imágen podemos ver la información mostrada en la pantalla del lcd : en la primera línea tenemos el último número primo encontrado, en la segunda línea tenemos durante cuantos minutos ha estado ejecutándose el programa y el número de primos encontrados.

Una mejora al programa sería lograr que se grabaran los primos encontrados en una memoria eeprom externa o en una tarjeta de memoria. Ya veremos ...

Actualización del 09/03/2015
En la imágen podemos ver que después de 9768 minutos ( 163 horas) se han encontrado 731.696 números primos.

Experimentos con Arduino

Páginas

domingo, 1 de marzo de 2015

Números primos grandes en Arduino