witness.c

Denna kod är public domain. Om ni hittar fel eller vill ändra något i koden blir jag jätteglad om ni skickar dessa ändringar till jesper [at] fantasi [punkt] se.

/*
 * Witness är en algoritm som kontrollerar om ett tal inte kan vara
 * ett primtal. Om algoritmen svarar 'sant' är talet garanterat inte
 * ett primtal, om resultatet är 'falskt' kan det vara ett primtal.
 *
 * Algoritmen anropas med två tal, slump och prima. prima är talet vi
 * vill testa om det kan vara ett primtal och slump är ett slumptal
 * som är lägre än prima. Eftersom slump ska vara lägre än prima och
 * det inte fungerar så bra med 0 och 1 så funkar denna implementation
 * av algoritmen inte jättebra på de lägsta talen (1 kommer att
 * klassas som "inte primtal").
 *
 * Sannorlikheten för att witness ska svara fel är mindre än 0,5.
 * Detta innebär att om algoritmen anropas ett antal gånger i rad med
 * olika slumptal och den svarar falskt varje gång så kan man till
 * slut vara ganska säker på att det är ett primtal.  Om algoritmen
 * anropas n gånger är sannorlikheten för att prima är ett slumptal
 * 1-2^(-n).
 *
 * Teorin bakom algoritmen härstammar från G. Miller (1975) och
 * M. Rabin (1980) och finns med största sannorlikhet att läsa i
 * närmaste bok om kryptologi.
 */

#include<stdio.h>
#include<time.h>
#include<stdlib.h>

static int witness(unsigned int slump, unsigned int prima)
{
   unsigned int bits = prima - 1;
   unsigned int d = 1;
   int i;

   for (i = 31; i >= 0; i--)
   {
      unsigned int x = d;
      d = (d * d) % prima;
      if (d == 1 && x != 1 && x != prima - 1)
         return 1;
      if ((bits >> i) % 2 == 1)
         d = (d * slump) % prima;
   }

   if (d != 1)
      return 1;

   return 0;
}

int main(int argc, char *argv[])
{
   unsigned int slump;
   unsigned int prima;

   if (argc != 2) 
   {
      fprintf(stderr, "Anropa programmet med det tal som ska testas.\n");
      exit(EXIT_FAILURE);
   }

   prima = (unsigned int)atoi(argv[1]);
   srand((unsigned int)time(NULL));
   slump = (unsigned int)((rand() / (double)RAND_MAX) * (prima - 3)) + 2;

   if (witness(slump,prima) == 0)
      printf("%u är troligen ett primtal\n", prima);
   else
      printf("%u är garanterat inte ett primtal\n", prima);

   return EXIT_SUCCESS;
}