Implementierung vom AKS-Test

[Guest]

Hallo,

da ich für mich als neues Hobby entdeckt habe alle möglichen Algorithmen zu implementieren und diese dann in einem Wiki zu sammeln hänge ich im moment beim AKS-Test[1] für Primzahlen.
Allerdings scheitere ich bereits in der ersten Zeile den dort steht.
if N=a^b mit b > 1
return "zusammengesetzt";

Nun stehe ich allerdings vor dem Problem das mir nicht so ganz klar ist, wie ich dieses N=a^b implementieren sollen wenn N die zu prüfende Zahl ist.
Meine Überlegungen gehen im moment dahin das ich einfach eine For-Schleife schreibe und überprüfe ob es eine b-te Wurzel gibt die ein ganz Zahliges Ergebnis von N liefert. Also das ich aus N die b-te Wurzel ziehe und eine ganze Zahl zurück bekomme. Allerdings bin mir da nicht sicher ob das wirklich so gemeint ist?
Ich habe nun die Hoffnung das evtl. hier im Forum jemand mit liest der diesen Algorithmus auch bereits implementiert hat und mir darum sagen kann wie die das wohl meinen?
Ich wäre über einen entsprechenden HInweis auf jeden Fall sehr dankbar.
Viele Grüsse
Dan




[1] http://de.wikipedia.org/wiki/AKS-Primzahltest

 

[Marco13]

Erster, spontaner Gedanke: Den Logarithmus von N zur Basis a ausrechnen, und schauen, ob das "genau" b ergibt. Aber warte lieber erstmal ab, was Andrey dazu sagt

 

[Landei]

Es muss bessere Verfahren geben. Z.B nehmen wir an, du hättest N % 7 == 3. Wenn du die zweiten bis vierten Potenzen einer Zahl betrachtest, haben diese durch 7 niemals den Rest 3, man kann sich also sparen, nach einer Quardat-, Kubik- oder vierten Wurzel zu suchen.

 

[0x7F800000]

uff. Marco13, du treibst mich hier noch in die totale Blamage^^ :lol:

Ich muss leider zugeben, dass mir absolut nichts brauchbares eingefallen ist, ob wohl ich jetzt rund ne stunde drüber nachgedacht hab. :cry: Leider ist die billige newton-iteration momentan das einzige, was ich da vorschlagen kann. Allerdings brauchst du nicht alle b-ten wurzeln zu testen, das ist unnötig. Du solltest bei 2 anfangen, und nur N^(1/primzahl) ausrechnen, d.h. du solltest dir nebenbei in einem sieb des eratosthenes vielfache von bereits getesteten primzahlen wegstreichen: wenn N != x^7 dann ist definitiv N != x^42 usw...
Da es etwa x/ln(x) primzahlen <x gibt, und du alle b aus [2,sqrt(N)] testen musst, und das newton-verfahren auch irgendwie in O(ln(b)) läuft, müsstest du wohl irgendwas in richtung O(2*ln(b)sqrt(N)/ln(N)) erhalten, das sieht irgendwie durch polynome abschätzbar aus, würde also die
"Primes in P" tatsache nicht kaputtmachen, wenn ich mich jetzt nicht zu sehr verhaun hab.

[edit]
Das ist jetzt irgendwie nicht besonders ausgefeilt bisher (muss drüber nochma schlafen^^) aber ich erzähl's mal trotzdem, nur um mich morgen wieder dran erinnern zu können:
Angenommen N=a^b. Suche dir irgendeine Primzahl p>N s.d. (p-1) in viele kleine faktoren zerbröselt.
Ordnung der Gruppe (Z/pZ)* ist ja genau (p-1), und es muss nach Satz von Lagrange gelten:
ord(a^b)=ord(N)|ord(G)=(p-1)
Dann geht man los, und versucht aus den faktorisierten primfaktoren irgendwas einigermaßen übersichtliches P=p1*p2*p3... zusammenzubauen sodass ord(N)|P d.h. N^P=1 gilt. Und dann weiß ich nicht was ich eigentlich machen wollte. Irgendwie hab ich so ein komisches gefühl in der magengegend, dass die übrig bleibenden faktoren und kombinationen daraus ganz gute kandidaten für das richtige b wären, wobei jetzt b die potenz der gezogenen wurzel sein soll. ???:L Aber so wirklich begründen kann ich das grad nicht.
Wenn es aber klappen würde, dann könnte man das mit mehreren threads für mehrere p durchprobieren, und so evtl schneller zum ergebnis kommen, als durch brutales durchprobieren aller b's.

Och sorry leuts, ich hab grad echt keinen Schimmer von gar nichts, ich tue nur so^^

 

[Guest]

guten morgen,

erst einmal vielen Dank für die vielen antworten. Ich habe gerade mal ein klein bißchen im Netz gestöbert und bin über diesen[1] Forumbeitrag gestolpert, offenbar genügt es für b das Intervall [2, log_2 N] zu testen...Allerdings bin ich mir da immer noch nicht wirklich sicher ob das wirklich effektiv ist, da das doch immer noch verteufelt viele Tests sind wenn die Primzahlen groß werden.
Viele Grüsse
Dan


[1] http://newsgroups.derkeiler.com/Archive/De/de.sci.informatik.misc/2007-01/msg00035.html

 

[Landei]

Das hier scheint zu funktionieren, müsste nur für BigInt umgeschrieben werden. Der Test gibt *mögliche* Wurzeln (Kandidaten) aus, also falls value eine Potenz ist, ist es eine der genannten:

import java.util.BitSet;
import static java.lang.Math.*;

public class Powertest {
  
  public static final int[] primes = {3,5,7,11,13,17,19,23,29,31,37,41,43,47};
  
  public static void test(long value) {
    int maxtest = (int) ceil(log(value) / log(2));
    BitSet bitset = new BitSet();
    for(int prime : primes) {
      int rem = (int)(value % (long)prime);
      if (rem == 0) {
        System.out.println("" + value + " % " + prime +" == 0");
        continue;
      }
      int array[] = new int[prime-1];
      for(int i = 1; i < prime; i++) {
         array[i-1] = i;
      }
      
      for(int power = 2; power < prime; power++) {
        for(int i = 1; i < prime; i++) {
          array[i-1] = (array[i-1]* i) % prime;
        }
        boolean notPossible = true;
        for(int i = 1; i < prime; i++) {
          if(array[i-1] == rem) {
            notPossible = false;
            break;
          }
        }
        if (notPossible) {
          //System.out.println("Found " + power + "th power mod" + prime);
          for (int b = power; b <= maxtest; b += prime - 1) {
            bitset.set(b);
          }
        }
      }
    }
    for(int i = 2; i <= maxtest; i++) {
      if (! bitset.get(i)) {
        System.out.println("Needs test for " + i + "th power");
      }
    }
  }
  
  public static void main(String... args) {
    System.out.println("Test for 418195493 = 53^5");
    test(418195493); //5,13,17,19,25
    System.out.println("\nTest for 2488651484819 = 59^7");
    test(2488651484819L); //7,13,17,19,29,31,37
    System.out.println("\nTest for 2488651484509 (keine Potenz)");
    test(2488651484509L); //13,17,19,29,31,37
  }
}

Wenn man davon ausgeht, dass die Zahl eigentlich prim sein sollte, es einem also auch reicht, einen Primfaktor zu finden, kann man bei der Berechnung von maxtest die obere Grenze mit der nächsten Primzahl, die nicht im Array ist (anstatt 2) vornehmen, um die Anzahl Kandidaten zu verkleinern: int maxtest = (int) ceil(log(value) / log(53));
Nach dieser Änderung werden für die drei Tests nur 5, 7 und gar nichts angezeigt.

 

[hawkeye78]

Hallo landei,

ich muss ehrlich sagen das ich den Code nicht verstehe. Du übergibtst eine Zahl die du überprüfen möchtest ob sie Prim ist und prüfst dann ob die zu prüfende Zahl ein vielfaches einer anderen Primzahl ist. Vielleicht fehlen mir die mathematisches Grundlagen, aber es gibt doch bestimmt auch Zahlen die das vielfache von nicht Primzahlen sind.
Möglicherweise kann mich ja in dieser Beziehung ja jemand aufklären.
viele Grüsse
Dan

PS
hawkeye78 = Gast nicht damit verwirrung gestiftet wird

 

[0x7F800000]

Allerdings bin ich mir da immer noch nicht wirklich sicher ob das wirklich effektiv ist, da das doch immer noch verteufelt viele Tests sind wenn die Primzahlen groß werden.

Dir ist schon bewusst, dass diese ganze sache als beweis für die tatsache "Primes in P" gewesen ist, und sich eigentlich kaum für den praktischen einsatz lohnt, sondern eher vom theoretischen interesse ist?

Übrigens, hab hier das da gefunden. Im Teil 4 steht da, dass dieser Test zumindest mal mit der newton-iteration erledigt werden kann, was meinen oben stehenden vorschlag etwas plausibler macht, und die hoffnung auf ein effizientes verfahren schwinden lässt.

 

[Marco13]

@Andrey: Die totale Blamage liegt eher bei mir: a ist ja garnicht gegeben

Aber eine Websuche liefert schnell Ergebnisse: In der LiDIA-Bibliothek ( http://www.cdc.informatik.tu-darmstadt.de/TI/LiDIA/ ) gibt es eine praktische funktion namens "is_power", die genau die gestellte Aufgabe löst. Eine genauere Erklärung gibt's vielleicht in der Doku. Aber @Andrey: Da kommt zumindest eine "newton_root" vor - scheint also als lagst du mit deiner ersten Idee schon richtig.

Ich bin mal so frei \lidia-2.2.0\src\base\simple_classes\bigint\is_power.cc und is_power.cc hier zu posten - wer es compilieren will kann sich die Bibliothek downloaden, aber für einen Java-Port müßte das (soweit ich das beim ersten Überfliegen beurteilen kann) schon reichen:

//==============================================================================================
//
//	This file is part of LiDIA --- a library for computational number theory
//
//	Copyright (c) 1994--2001 the LiDIA Group.  All rights reserved.
//
//	See [url]http://www.informatik.tu-darmstadt.de/TI/LiDIA/[/url]
//
//----------------------------------------------------------------------------------------------
//
//	$Id: is_power.cc,v 2.2 2004/06/15 10:19:25 lidiaadm Exp $
//
//	Author	: Thomas Papanikolaou (TP)
//	Changes	: See CVS log
//
//==============================================================================================


#ifdef HAVE_CONFIG_H
# include	"config.h"
#endif
#include	"LiDIA/bigint.h"



#ifdef LIDIA_NAMESPACE
namespace LiDIA {
#endif



#define prime_table_length 100



static int prime_table[prime_table_length] = {
	2, 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, 257, 263, 269, 271, 277, 281, 283,
	293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,
	383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461,
	463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541};



long
bigint::is_power(bigint & b) const
{
	bigint c;
	int i = 0;

	if (!is_negative()) {
		int l = bit_length();

		for (; i < prime_table_length && prime_table[i] <= l; i++) {
			newton_root(b, *this, prime_table[i]);

			power(c, b, prime_table[i]);
			if (compare(c) == 0)
				return prime_table[i];
		}
		for (i = 545; i <= l; i += 4) {
			// testing 6k-1 and 6k+1 for  k >= 91
			newton_root(b, *this, i);
			power(c, b, i);
			if (compare(c) == 0)
				return i;
			i += 2; // increase from 6k-1 to 6k+1. Do again.
			newton_root(b, *this, i);
			power(c, b, i);
			if (compare(c) == 0)
				return i;
		}
	}
	else {
		i++; // negative numbers can't be even powers
		bigint a2(*this);
		a2.abs();

		int l = a2.bit_length();

		for (; i < prime_table_length && prime_table[i] <= l; i++) {
			newton_root(b, a2, prime_table[i]);
			power(c, b, prime_table[i]);
			if (a2.compare(c) == 0) {
				b.negate();
				return prime_table[i];
			}
		}
		for (i = 545; i <= l; i += 4) {
			// testing 6k-1 and 6k+1 for  k >= 91
			newton_root(b, a2, i);
			power(c, b, i);
			if (a2.compare(c) == 0) {
				b.negate();
				return i;
			}
			i += 2; // increase from 6k-1 to 6k+1. Do again.
			newton_root(b, a2, i);
			power(c, b, i);
			if (a2.compare(c) == 0) {
				b.negate();
				return i;
			}
		}
	}
	return 0;
}



#ifdef LIDIA_NAMESPACE
}	// end of namespace LiDIA
#endif

//==============================================================================================
//
//	This file is part of LiDIA --- a library for computational number theory
//
//	Copyright (c) 1994--2001 the LiDIA Group.  All rights reserved.
//
//	See [url]http://www.informatik.tu-darmstadt.de/TI/LiDIA/[/url]
//
//----------------------------------------------------------------------------------------------
//
//	$Id: newton_root.cc,v 2.5 2004/06/15 10:19:26 lidiaadm Exp $
//
//	Author	: Thomas Papanikolaou (TP)
//	Changes	: See CVS log
//
//==============================================================================================


#ifdef HAVE_CONFIG_H
# include	"config.h"
#endif
#include	"LiDIA/bigint.h"



#ifdef LIDIA_NAMESPACE
namespace LiDIA {
#endif



void
newton_root(bigint & b, const bigint & a, int n)
{
	bigint c, d;

	if (n < 0)
		lidia_error_handler("bigint", "newton_root::negative n.");

	if (a.is_negative())
		lidia_error_handler("bigint", "newton_root::negative a.");

	switch (n) {
	case 0:
		b.assign_one();
		break;
	case 1:
		b.assign(a);
		break;
	case 2:
		sqrt(b, a);
		break;
	default:
		b.assign_one();
		shift_left(b, b, static_cast<unsigned int>((a.bit_length() + n - 1) / n));
		do {
			power(c, b, n - 1);
			div_rem(c, d, a, c);
			subtract(c, b, c);
			divide(c, c, static_cast<long>(n));
			subtract(b, b, c);
		} while (c.sign() > 0);
		power(c, b, n);
		if (c.compare(a) > 0)
			dec(b);
		if (b.compare(3UL) == 0) {
			power(c, b, n);
			if (c.compare(a) > 0)
				dec(b);
		}
		break;
	}
}



#ifdef LIDIA_NAMESPACE
}	// end of namespace LiDIA
#endif

 

[0x7F800000]

Also, ich hab heute beinahe rein zufällig ein Buch "Elementare und Algebraische Zahlentheorie" vom Herrn Prof. Müller-Stach in die Hände gekrigt, dort war dieser Algorithmus beschrieben. Ich hab nur ganz kurz reingeschaut, aber da stand tatsächlich, dass man für alle b aus [2,log_2(N)] numerisch die nullstellen von x^b-1=N ausrechnet, und dann eben prüft, ob das ergebnis ganzzahlig ist, bzw potenziert es wieder hoch b um sicherzustellen, dass a^b=N ergibt.
(korrektur: srtt(N) war da oben natürlich nicht richtig, das war die schranke für die Basis a nicht für den exponenten b, sorry)

Und ich möchte wiederholt darauf hinweisen, dass das einzige tolle an diesem test ist, dass er garantiert in polynomieller laufzeit terminiert. Der ist verglichen zu den üblichen randomisierten Algos relativ langsam, dafür können die randomisierte Algos eben nicht garantieren, dass die für irgendwelche ausnahmefälle schnell eine Lösung finden, das kommt aber extrem selten vor (besonders wenn man mehrere randomisierte verfahren mit verschiedenen startwerten nebeneinander laufen lässt, kommt man meistens schneller zu ziel)

 

[Landei]

Hallo landei,
ich muss ehrlich sagen das ich den Code nicht verstehe. Du übergibtst eine Zahl die du überprüfen möchtest ob sie Prim ist und prüfst dann ob die zu prüfende Zahl ein vielfaches einer anderen Primzahl ist.

Ich will nicht testen, ob die Zahl Prim ist, sondern, ob es sich eine Potenz handelt - was man erst einmal ausschließen muss, wenn mit dem "eigentlichen" AKS-Test anfangen will

Vielleicht fehlen mir die mathematisches Grundlagen, aber es gibt doch bestimmt auch Zahlen die das vielfache von nicht Primzahlen sind.

Jede natürliche Zahl hat eine eindeutige Zerlegung in Primfaktoren. Das spielt an dieser Stelle aber keine Rolle.

Der Grundgedanke von obigem Code ist eigentlich ganz einfach. Nehmen wir mal die Quadratzahlen:
1,4,9,16,25,36...
Jetzt bilden wir die Reste bei Division durch 3:
1,1,0,1,1,0,1...
Wir vermuten mal ganz scharf, das bei der Division einer Quadratzahl durch 3 niemals der Rest 2 herauskommt (der Beweis ist trivial). Wenn wir also die Zahl 1234567891234567891234567891234567891234567892 haben, die den Rest 2 läßt, wissen wir haargenau, dass es keine Quadratzahl sein kann - ohne dass wir die Wurzel wirklich ausgerechnet hätten. Wenn eine Zahl keine Quadratzahl ist, ist sie natürlich erst recht keine vierte, sechste, achte, zehnte Potenz u.s.w. - die können wir uns auch sparen zuberechnen.
Mein Code stellt solche Regeln nicht nur für die Division durch 3, sondern auch durch andere Primzahlen auf (zusammengesetzte Zahlen sind hier uninteressant, da sie die gleichen Regeln liefern wie ihre Primfaktoren). Und dann bestimme ich die Regeln nicht nur für Quadratzahlen, sondern für dritte, vierte, fünfte... Potenzen. Damit siebe ich einfach aus, um welche Potenzen es sich bei unserer Zahl mit Sicherheit nicht handelt, und reduziere damit die Zahl der Wurzeln, die ich tatsächlich berechnen muss.

 

[0x7F800000]

Und dann bestimme ich die Regeln nicht nur für Quadratzahlen, sondern für dritte, vierte, fünfte... Potenzen.

könntest du vielleicht nochmal in worten erklären, wie du das "erstellen von tests" bewerkstelligst? ???:L aus deinem code werde ich auch nicht wirklich schlau :bahnhof: