Primzahlen

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
gast

Primzahlen

Beitragvon gast » Sonntag 2. Mai 2004, 14:31

Hallo allerseits,

ich habe vor in Python Primzahlen zu berechnen. Im allgemeinen ist das natürlich kein Problem. Ich möchte jedoch den AKS-Algorithmus der 2002 gefunden wurde verwenden. Leider habe ich den bisher nur in C/C++ gesehen. Verionen kann man sich hier anschauen:
http://fatphil.org/maths/AKS/#Implementations
Meine frage ist nun, ob jemand weiss wie man den AKS-Alogithmus in Python implementiert.

Vielen Dank schonmal im Vorraus.
Benutzeravatar
Dookie
Python-Forum Veteran
Beiträge: 2010
Registriert: Freitag 11. Oktober 2002, 18:00
Wohnort: Salzburg
Kontaktdaten:

Beitragvon Dookie » Sonntag 2. Mai 2004, 16:13

Hallo Gast,

hast Du mal auf Deinen Link geklickt? Bei mir kommt dann diese Fehlerseite:
http://fatphil.org/errors/403.html#Implementations


Gruß

Dookie
Milan
User
Beiträge: 1078
Registriert: Mittwoch 16. Oktober 2002, 20:52

Beitragvon Milan » Sonntag 2. Mai 2004, 16:53

Hi. Bei mir kommt eine Seite voller Links zu Implementationen in C und C++, ganz wie beschrieben... aber übersetzen kann ich das aus 2 Gründen nicht: a) ich kenn den AKS Algorhytmus nicht und b) da werden einige Funktionen benutzt die weder im Text definiert sind oder mir algemein bekannt sind (Bsp.: nextprime,IsPower ...).

Vielleicht wäre eine deutschsprachige Beschreibung zu ASK nicht schlecht...

Milan
Benutzeravatar
Dookie
Python-Forum Veteran
Beiträge: 2010
Registriert: Freitag 11. Oktober 2002, 18:00
Wohnort: Salzburg
Kontaktdaten:

Beitragvon Dookie » Sonntag 2. Mai 2004, 17:59

Hi nochmal,

hab jetzt mal google auf "AKS-Algorithmus Primzahlen" da kommen doch einige Treffer.
http://www.google.de/search?q=AKS-Algor ... Primzahlen


Gruß

Dookie
Benutzeravatar
reggid
User
Beiträge: 120
Registriert: Dienstag 8. Oktober 2002, 19:04
Wohnort: Dinslaken
Kontaktdaten:

Beitragvon reggid » Sonntag 2. Mai 2004, 19:18

Zum Thema umsetzung in python
Die C++ Umsetzungen benutzen nur die Standard-Libs und was die Math lib von C kann, kann die von Python schon lange, afaik!
Benutzeravatar
Dookie
Python-Forum Veteran
Beiträge: 2010
Registriert: Freitag 11. Oktober 2002, 18:00
Wohnort: Salzburg
Kontaktdaten:

Beitragvon Dookie » Sonntag 2. Mai 2004, 19:18

hmm hab mir jetzt ein paar Treffer angeschaut, aber ausser dem Link, auf den wir Ösis nicht dürfen, gibts wohl keine Codebeispiele.


Gruß

Dookie
Milan
User
Beiträge: 1078
Registriert: Mittwoch 16. Oktober 2002, 20:52

Beitragvon Milan » Montag 3. Mai 2004, 16:44

Hi. Stimmt, jetzt hab ich's auch gesehn (1. Beispiel @reggid). Mal schauen, wenn ich mal Zeit finde fange ich mit übersetzen an, muss aber erst mal Schule (jaja, es ist Klausurzeit) machen. Für alle anderen, die nicht zum Quelltext kommen, aber trotzdem was machen möchten... hier:

Code: Alles auswählen

#include <iostream>
using namespace std;
#include <math.h>

#ifdef _M_IX86

#define umulrem(z, x, y, m)   \
__asm   mov      eax, x   \
__asm   mul      y   \
__asm   div      m   \
__asm   mov      z, edx

#define umuladdrem(z, x, y, a, m)   \
__asm   mov      eax, x   \
__asm   mul      y   \
__asm   add      eax, a   \
__asm   adc      edx, 0   \
__asm   div      m   \
__asm   mov      z, edx

#else

#ifdef _MSC_VER
typedef unsigned __int64   Tu64;
#else
typedef unsigned long long   Tu64;
#endif

#define umulrem(z, x, y, m)   \
   {   \
   z = (unsigned int)(x * (Tu64)y % m);   \
   }

#define umuladdrem(z, x, y, a, m)   \
   {   \
   z = (unsigned int)((x * (Tu64)y + a) % m);   \
   }

#endif

static bool IsPrime(unsigned int n)
{
   if (n < 2) return false;
   if (n < 4) return true;
   if (n % 2 == 0) return false;

   const unsigned int iMax = (int)sqrt(n) + 1;
   unsigned int i;
   for (i = 3; i <= iMax; i += 2)
      if (n % i == 0)
         return false;

   return true;
}

static unsigned int LargestPrimeFactor(unsigned int n)
{
   if (n < 2) return 1;

   unsigned int r = n, p;
   if (r % 2 == 0)
   {
      p = 2;
      do { r /= 2; } while (r % 2 == 0);
   }
   unsigned int i;
   for (i = 3; i <= r; i += 2)
   {
      if (r % i == 0)
      {
         p = i;
         do { r /= i; } while (r % i == 0);
      }
   }
   return p;
}

static unsigned int Powm(unsigned int n, unsigned int e, unsigned int m)
{
   unsigned int r = 1;
   unsigned int t = n % m;
   unsigned int i;
   for (i = e; i != 0; i /= 2)
   {
      if (i % 2 != 0)
      {
         umulrem(r, r, t, m);
      }
      umulrem(t, t, t, m);
   }
   return r;
}

static unsigned int Inv(unsigned int n, unsigned int m)
{
   unsigned int a = n, b = m;
   int u = 1, v = 0;
   do
   {
      const unsigned int q = a / b;

      const unsigned int t1 = a - q*b;
      a = b;
      b = t1;

      const int t2 = u - (int)q*v;
      u = v;
      v = t2;
   } while (b != 0);
   if (a != 1) u = 0;
   if (u < 0) u += m;
   return u;
}

class CPolyMod
{
protected:
   // (mod x^r - 1, n)
   const unsigned int m_r;
   const unsigned int m_n;
   unsigned int m_deg;
   unsigned int * mp_a;

private:
   CPolyMod():m_r(0), m_n(0) { mp_a = NULL; };

public:
   // default value is x
   CPolyMod(unsigned int r, unsigned int n)
      : m_r(r), m_n(n)
   {
      m_deg = 1;
      mp_a = new unsigned int [2];
      mp_a[0] = 0; mp_a[1] = 1;
   }

   CPolyMod(const CPolyMod & p)
      : m_r(p.m_r), m_n(p.m_n)
   {
      m_deg = p.m_deg;
      mp_a = new unsigned int [p.m_deg + 1];
      unsigned int i;
      for (i = 0; i <= p.m_deg; ++i)
         mp_a[i] = p.mp_a[i];
   }

   virtual ~CPolyMod()
   {
      if (mp_a != NULL)
         delete [] mp_a;
   }

private:
   void _polySquare()
   {
      const unsigned int deg = m_deg;
      const unsigned int n = m_n;
      const unsigned int * const p_a = mp_a;

      const unsigned int degr = deg + deg;
      unsigned int * const p_ar = new unsigned int [degr + 1];
      unsigned int k;
      for (k = 0; k <= degr; ++k)
         p_ar[k] = 0;

      unsigned int j;
      for (j = 1; j <= deg; ++j)
      {
         const unsigned int x = p_a[j];
         if (x != 0)
         {
            unsigned int i;
            for (i = 0; i < j; ++i)
            {
               const unsigned int y = 2 * p_a[i];
               unsigned int t = p_ar[j + i];
               umuladdrem(t, x, y, t, n);
               p_ar[j + i] = t;
            }
         }
      }
      unsigned int i;
      for (i = 0; i <= deg; ++i)
      {
         const unsigned int x = p_a[i];
         unsigned int t = p_ar[2 * i];
         umuladdrem(t, x, x, t, n);
         p_ar[2 * i] = t;
      }

      m_deg = degr;
      delete [] mp_a;
      mp_a = p_ar;
   }

   void _polyMul(const CPolyMod & p)
   {
      const unsigned int deg = m_deg;
      const unsigned int n = m_n;
      const unsigned int * const p_a = mp_a;

      const unsigned int degr = deg + p.m_deg;
      unsigned int * const p_ar = new unsigned int [degr + 1];
      unsigned int k;
      for (k = 0; k <= degr; ++k)
         p_ar[k] = 0;

      unsigned int j;
      for (j = 0; j <= p.m_deg; ++j)
      {
         const unsigned int x = p.mp_a[j];
         if (x != 0)
         {
            unsigned int i;
            for (i = 0; i <= deg; ++i)
            {
               const unsigned int y = p_a[i];
               unsigned int t = p_ar[j + i];
               umuladdrem(t, x, y, t, n);
               p_ar[j + i] = t;
            }
         }
      }

      m_deg = degr;
      delete [] mp_a;
      mp_a = p_ar;
   }

   void _Mod()
   {
      unsigned int deg = m_deg;
      unsigned int * const p_a = mp_a;
      while (deg >= m_r)
      {
         p_a[deg - m_r] += p_a[deg];
         if (p_a[deg - m_r] >= m_n) p_a[deg - m_r] -= m_n;
         --deg;

         while (p_a[deg] == 0) --deg;
      }
      m_deg = deg;
   }

   void _Norm()
   {
      const unsigned int deg = m_deg;
      const unsigned int n = m_n;
      unsigned int * const p_a = mp_a;
      if (p_a[deg] != 1)
      {
         const unsigned int y = Inv(p_a[deg], m_n);
         unsigned int i;
         for (i = 0; i <= deg; ++i)
         {
            unsigned int t = p_a[i];
            umulrem(t, t, y, n);
            p_a[i] = t;
         }
      }
   }

public:
   CPolyMod & operator = (const CPolyMod & p)
   {
      if (&p == this) return *this;
      m_deg = p.m_deg;
      delete [] mp_a;
      mp_a = new unsigned int [p.m_deg + 1];
      unsigned int i;
      for (i = 0; i <= p.m_deg; ++i)
         mp_a[i] = p.mp_a[i];
      return *this;
   }

   int operator != (const CPolyMod & p) const
   {
      if (m_deg != p.m_deg)
         return true;
      unsigned int i;
      for (i = 0; i <= m_deg; ++i)
         if (mp_a[i] != p.mp_a[i])
            return true;
      return false;
   }

   CPolyMod & operator += (unsigned int i)
   {
      const unsigned int t = i % m_n;
      mp_a[0] += t;
      if (mp_a[0] >= m_n) mp_a[0] -= m_n;
      return *this;
   }

   CPolyMod & operator -= (unsigned int i)
   {
      const unsigned int t = m_n - i % m_n;
      mp_a[0] += t;
      if (mp_a[0] >= m_n) mp_a[0] -= m_n;
      return *this;
   }

   CPolyMod Pow(unsigned int e) const
   {
      unsigned int er = 1;
      unsigned int j;
      for (j = e; j != 1; j /= 2)
      {
         er = 2 * er + (j % 2);
      }

      CPolyMod t(*this);
      unsigned int i;
      for (i = er; i != 1; i /= 2)
      {
         t._polySquare();
         t._Mod();
         if (i % 2 != 0)
         {
            t._polyMul(*this);
            t._Mod();
         }
      }
      t._Norm();
      return t;
   }
};

int main(int argc, char * argv[])
{
// AKS
// n=11701, r=11699, q=5849, s=2923
// n=1000000007, r=57287, q=28643, s=14311
//
// Bernstein
// n=349, r=347, q=173, s=140
// n=1000000007, r=3623, q=1811, s=1785

   unsigned int n0;
   cout << "n ? ";
   cin >> n0;

   unsigned int n;
   for (n = n0; true; n += 2)
   {
      bool b = false;
      unsigned int q, r, s;
      for (r = 3; r < n; r += 2)
      {
         if (IsPrime(r))
         {
            const unsigned int m = n % r;
            if (m == 0) break;
            q = LargestPrimeFactor(r - 1);
            if (Powm(m, (r - 1) / q, r) <= 1) continue;
/*
            // AKS
            if (q >= 4 * sqrt(r) * n.Log()/log(2))
            {
               s = (unsigned int)(2 * sqrt(r) * n.Log()/log(2));
               b = true;
               break;
            }
*/
            // Bernstein
            const double cMin = 2 * floor(sqrt(r)) * log(n);
            double c = 0;
            for (s = 1; s <= q; s++)
            {
               c += log(q + s - 1) - log(s);
               if (c >= cMin)
               {
                  b = true;
                  break;
               }
            }
            if (b) break;

         }
      }
      if (b)
      {
         cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << "\r" << flush;;

         bool b = true;

         CPolyMod rPoly(r, n);            // x
         try
         {
            rPoly = rPoly.Pow(n);         // x^n
         }
         catch (...)
         {
            b = false;
         }
         if (b)
         {
            unsigned int a;
            for (a = 1; a <= s; ++a)
            {
               cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << " " << a << "\r" << flush;;

               CPolyMod lPoly(r, n);      // x
               lPoly -= a;               // x - a
               try
               {
                  lPoly = lPoly.Pow(n);   // (x - a)^n
               }
               catch (...)
               {
                  b = false;
                  break;
               }
               lPoly += a;               // (x - a)^n + a
               if (lPoly != rPoly)
               {
                  b = false;
                  break;
               }
            }
         }
         if (b)
            cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << ", power of a prime." << endl;
         else
            cout << "n=" << n << ", r=" << r << ", q=" << q << ", s=" << s << ", composite." << endl;
      }
   }

   return 0;
}


mfg, Milan

Wer ist online?

Mitglieder in diesem Forum: Bing [Bot]