Klasse zur Primfactor-zerlegung und zum primzahlen errechnen

Stellt hier eure Projekte vor.
Internetseiten, Skripte, und alles andere bzgl. Python.
nuss
User
Beiträge: 53
Registriert: Donnerstag 28. August 2008, 11:36

Beitragvon nuss » Mittwoch 29. Oktober 2008, 17:46

so hab noch'n paar Sachen verbessert:
http://paste.pocoo.org/show/89434/

@EyDu:
Hab jetzt in der Funktion get_new_primes proof und cache auf None gesetzt.
Die anderen Funktionen können ruhig die errechneten Primzahlen weiterverwenden.


Wiie kriegst du das raus, obs Tabs oder Spaces sind ?
Wenn ich python -tt ausführe, das Modul exportiere,
und alle Funktionen einmal aufrufe, gibt er mir jedenfalls keine Fehler aus.

man könnte die funktionen get_factors und isprime so vereinen:

Code: Alles auswählen

def _prime_ops(number, split, cache=[2]):
   
    proof = _get_first(cache)
    factors = []
    counter = 0

    while True:
        try:
            prime = cache[counter]
        except IndexError:

            cache, proof = get_new_primes(proof+1, proof, cache)
            continue

        if number % prime is 0:
            if split is False:
                return (False, cache)
            else:
               factors.append(prime)
               number /= prime
        elif number / prime < prime:
            if split is False:
                return (True, cache)
            else:
                if number is not 1:
                   factors.append(number)
                return (factors, cache)
        else:
           counter += 1

def get_factors(number, cache=[2]):
    return _prime_ops(number, True, cache)

def isprime(number, cache=[2]):
    return _prime_ops(number, False, cache)
Benutzeravatar
name
User
Beiträge: 254
Registriert: Dienstag 5. September 2006, 16:35
Wohnort: Wien
Kontaktdaten:

Beitragvon name » Freitag 31. Oktober 2008, 17:03

Schmerz! Need a Medic!
Du ueberpruefst Werte mit is, das ist ganz ganz boese. "is" prueft Identitaet, "==" den Wert.
PS: Nicht das ich euren Spass verderben will, aber wie waer es damit:

Code: Alles auswählen

def get_primes(max_):
    sieve = [True] * (max_ + 1)
    sieve[0] = sieve[1] = False
    for cur_number in xrange(2, max_ + 1):
        if sieve[cur_number]:
            j = 2 * cur_number
            while j <= max_:
                sieve[j] = False
                j += cur_number
            yield cur_number
Ohloh | Mein Blog | Jabber: segfaulthunter@swissjabber.eu | asynchia – asynchrone Netzwerkbibliothek

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.
Qubit
User
Beiträge: 75
Registriert: Dienstag 7. Oktober 2008, 09:07

Beitragvon Qubit » Samstag 1. November 2008, 02:27

name hat geschrieben:Schmerz! Need a Medic!
Du ueberpruefst Werte mit is, das ist ganz ganz boese. "is" prueft Identitaet, "==" den Wert.
PS: Nicht das ich euren Spass verderben will, aber wie waer es damit:


Ein Thema: Primzahlen finden.
Anderes Thema: Primfaktoren einer Zahl finden

Das erste ist simpel und dafür hast du eine mögliche Lösung geboten.
Das andere ist alles andere als simple. Dafür hast du keine Lösung geboten und alle anderen hier sind ziemlich "unperformant".
Um das zu sehen, nehme man mal so "einfache" Zahlen wie:
348234823482389479
oder
3482348234823894793453453455


Natürlich kann man darüber streiten, ob Python überhaupt für solche Kategorien numerischer Probleme geeignet ist. Aber mit einem halbwegs vernünftigen Algorithmus lässt sich auch hier etwas rausholen (mal mit den obigen Zahlen testen) ;-)

Code: Alles auswählen

import math
def primtest(z):
        """
        Primfaktorenzerlegung
        (c)2008 Qubit
        hannover.post@gmail.com
        """
       
        primfak=[]
        m = 0
        found = 0
        limit = int((math.sqrt(z)-1)/2)
        while z % 2 == 0:
                z = z/2
                primfak.append((2,z))
                print "found:",primfak[-1]
        i = (z-1)/2
        while (i-1) % 3 == 0:
                i = (i - 1)/3
                primfak.append((3,2*i+1))
                print "found:",primfak[-1]
        if (i-2) % 3 == 0:
                u = (i-2)/3
                s = 2
        if (i-3) % 3 == 0:
                u = (i-3)/3
                s = 3
        while m < u and m < limit:
                if s == 2:
                        if not found and (u-m) % (6*m+5) == 0:
                                n = (u-m) / (6*m+5)
                                m0 = 3*m + 2
                                n0 = 3*n
                                primfak.append((2*m0+1,2*n0+1))
                                print "found:",primfak[-1]
                                found = 1
                        if not found and (u-5*m-5) % (6*m+7) == 0:
                                n = (u-5*m-5) / (6*m+7)
                                m0 = 3*m + 3
                                n0 = 3*n + 2
                                primfak.append((2*m0+1,2*n0+1))
                                print "found:",primfak[-1]
                                found = 1
                if s == 3:
                        if not found and (u-5*m-3) % (6*m+5) == 0:
                                n = (u-5*m-3) / (6*m+5)
                                m0 = 3*m + 2
                                n0 = 3*n + 2
                                primfak.append((2*m0+1,2*n0+1))
                                print "found:",primfak[-1]
                                found = 1
                        if not found and (u-m) % (6*m+7) == 0:
                                n = (u-m) / (6*m+7)
                                m0 = 3*m + 3
                                n0 = 3*n
                                primfak.append((2*m0+1,2*n0+1))
                                print "found:",primfak[-1]
                                found = 1
                if found == 1:
                        i = n0
                        while (i-1) % 3 == 0:
                                i = (i - 1)/3
                                primfak.append((3,2*i+1))
                                print "found:",primfak[-1]
                        if (i-2) % 3 == 0:
                                u = (i-2)/3
                                s = 2
                        if (i-3) % 3 == 0:
                                u = (i-3)/3
                                s = 3
                        limit = int((math.sqrt(2*i+1)-1)/2)
                        found = 0
                        m -= 1
                m += 1

        if len(primfak) > 0:
                x = [x[0] for x in primfak]+[primfak[-1][1]]
                if x[-1] == 1: x.pop()
        else: x=[z]
        print "Primfaktoren: ",sorted(x)
Benutzeravatar
name
User
Beiträge: 254
Registriert: Dienstag 5. September 2006, 16:35
Wohnort: Wien
Kontaktdaten:

Beitragvon name » Samstag 1. November 2008, 19:30

Qubit hat geschrieben:

Code: Alles auswählen

....
....

Scheint zu funktionieren, schaut aber eher, naja, nicht sehr uebersichtlich aus. Koenntest du den Algorithmus bitte einfach so, bzw mit Pseudocode erklaeren?
Ohloh | Mein Blog | Jabber: segfaulthunter@swissjabber.eu | asynchia – asynchrone Netzwerkbibliothek

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.
Qubit
User
Beiträge: 75
Registriert: Dienstag 7. Oktober 2008, 09:07

Beitragvon Qubit » Sonntag 2. November 2008, 00:38

name hat geschrieben:
Qubit hat geschrieben:

Code: Alles auswählen

....
....

Scheint zu funktionieren, schaut aber eher, naja, nicht sehr uebersichtlich aus. Koenntest du den Algorithmus bitte einfach so, bzw mit Pseudocode erklaeren?

Nun, die mathematischen Hintergründe halte ich in diesem Forum für "beyond the scope", nur soviel, es ist etwas "Selbstgestricktes" und soll leicht portierbar sein ;-)

Habe es noch mal als Klasse implementiert..

Code: Alles auswählen

import math,time

class primetest(object):
        """
        Primfaktorenzerlegung, Primzahlentest, Primzahlensuche
        (c)2008 Qubit
        hannover.post@gmail.com
        """
        __time = {
                'fak': 0,
                'is_prime': 0,
                'primes': 0
                }

        def __getattr__(self,attr):
                attrs = {
                        't_fak': self.__time['fak'],
                        't_is_prime': self.__time['is_prime'],
                        't_primes': self.__time['primes']
                        }
                get = attrs.get(attr,None)
                if not get == None:
                        return tuple(str('%.3f,sec' %(get)).split(','))
                else:
                        object.__getattribute__(self,attr)
               
        def __fak(self,z,mode):
                primfak=[]
                m = 0
                found = 0
                limit = int((math.sqrt(z)-1)/2)
                while z % 2 == 0:
                        if mode == 1 and not z == 2: return [z,False]
                        z = z/2
                        primfak.append((2,z))
                i = (z-1)/2
                while (i-1) % 3 == 0:
                        if mode == 1 and not i == 1: return [z,False]
                        i = (i - 1)/3
                        primfak.append((3,2*i+1))
                if (i-2) % 3 == 0:
                        u = (i-2)/3
                        s = 2
                if (i-3) % 3 == 0:
                        u = (i-3)/3
                        s = 3
                while m < u and m < limit:
                  if s == 2:
                          if (u-m) % (6*m+5) == 0:
                                n = (u-m) / (6*m+5)
                                m0 = 3*m + 2
                                n0 = 3*n
                                primfak.append((2*m0+1,2*n0+1))
                                found = 1
                          if not found and (u-5*m-5) % (6*m+7) == 0:
                                n = (u-5*m-5) / (6*m+7)
                                m0 = 3*m + 3
                                n0 = 3*n + 2
                                primfak.append((2*m0+1,2*n0+1))
                                found = 1
                  if s == 3:
                          if (u-5*m-3) % (6*m+5) == 0:
                                n = (u-5*m-3) / (6*m+5)
                                m0 = 3*m + 2
                                n0 = 3*n + 2
                                primfak.append((2*m0+1,2*n0+1))
                                found = 1
                          if not found and (u-m) % (6*m+7) == 0:
                                n = (u-m) / (6*m+7)
                                m0 = 3*m + 3
                                n0 = 3*n
                                primfak.append((2*m0+1,2*n0+1))
                                found = 1
                  if found == 1:
                          if mode == 1 and not 2*m0+1 == z:
                                  return [z,False]
                          i = n0
                          while (i-1) % 3 == 0:
                                i = (i - 1)/3
                                primfak.append((3,2*i+1))
                          if (i-2) % 3 == 0:
                                u = (i-2)/3
                                s = 2
                          if (i-3) % 3 == 0:
                                u = (i-3)/3
                                s = 3
                          limit = int((math.sqrt(2*i+1)-1)/2)
                          found = 0
                          m -= 1
                  m += 1

                if len(primfak) > 0:
                        x = [x[0] for x in primfak]+[primfak[-1][1]]
                        if x[-1] == 1: x.pop()
                else: x=[z]
                return x

        def __is_prime(self,z): return len(self.__fak(z,mode=1)) == 1
        def __primes(self,bis,von): return filter(self.__is_prime, range(von,bis+1))

        def fak(self,z=2):
                t_start = time.clock()
                res = self.__fak(z,mode=0)
                t_elapsed = time.clock() - t_start
                self.__time['fak'] = t_elapsed
                return res
        def is_prime(self,z=2):
                t_start = time.clock()
                res = self.__is_prime(z)
                t_elapsed = time.clock() - t_start
                self.__time['is_prime'] = t_elapsed
                return res
        def primes(self,bis=2,von=2):
                t_start = time.clock()
                res = self.__primes(bis,von)
                t_elapsed = time.clock() - t_start
                self.__time['primes'] = t_elapsed
                return res


>>> prim=primetest()
>>> prim.is_prime(348234823482389479)
False
>>> prim.t_is_prime
('0.350', 'sec')
>>> prim.fak(348234823482389479)
[413243, 842687773253L]
>>> prim.t_fak
('2.057', 'sec')
>>> prim.is_prime(123456789)
False
>>> prim.t_is_prime
('0.000', 'sec')
>>> prim.fak(123456789)
[3, 3, 3607, 3803]
>>> prim.t_fak
('0.001', 'sec')
>>> res = prim.primes(10000)
>>> prim.t_primes
('0.142', 'sec')
>>> res[-10:]
[9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973]
>>> len(res)
1229
>>>
audax
User
Beiträge: 830
Registriert: Mittwoch 19. Dezember 2007, 10:38

Beitragvon audax » Sonntag 2. November 2008, 02:07

Ich erzittere vor deinen Mathkenntnissen. Trotzdem ist der Code hässlich wie dar Tag und die Klasse mehr als überflüssig.

Vor allem ist in dem Code irre viel Wiederholung und keine einzige Erklärung.
Benutzeravatar
name
User
Beiträge: 254
Registriert: Dienstag 5. September 2006, 16:35
Wohnort: Wien
Kontaktdaten:

Beitragvon name » Sonntag 2. November 2008, 11:37

audax hat geschrieben:Ich erzittere vor deinen Mathkenntnissen. Trotzdem ist der Code hässlich wie dar Tag und die Klasse mehr als überflüssig.

Vor allem ist in dem Code irre viel Wiederholung und keine einzige Erklärung.
Amen. Wuerde mich echt voll intressieren wie das so schnell geht.
Ohloh | Mein Blog | Jabber: segfaulthunter@swissjabber.eu | asynchia – asynchrone Netzwerkbibliothek

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.
Qubit
User
Beiträge: 75
Registriert: Dienstag 7. Oktober 2008, 09:07

Beitragvon Qubit » Sonntag 2. November 2008, 15:08

audax hat geschrieben:Ich erzittere vor deinen Mathkenntnissen. Trotzdem ist der Code hässlich wie dar Tag und die Klasse mehr als überflüssig.

Warte erst mal ab, bis die richtigen Zahlentheoretiker kommen ;-)
Ja, der Code ist nicht hübsch, aber er ist portable und erfüllt seinen Zweck.
Die Klasse ist auch eher "standalone", mir fallen nicht allzu viele (praktische) Anwendungen ein, wo man Primzahlen bräuchte.
Aber zB bei Brüchen zum Kürzen ;-):

Code: Alles auswählen

class ratio(object):
    def __init__(self,zaehler,nenner):
        if nenner == 0:
            raise ZeroDivisionError
        else:
            abs_n = abs(int(nenner))
            sign_n = int(nenner) /abs_n           
        if zaehler == 0:
            sign_z = 1
            r =(0,abs_n)
        else:
            abs_z = abs(int(zaehler))
            sign_z = int(zaehler)/abs_z
            r = self.__reduce(abs_z,abs_n)
        self.z = sign_z*divint(r[0])
        self.n = sign_n*divint(r[1])
    def __add__(x,y):
        return divint(x.z*y.n+y.z*x.n)/divint(x.n*y.n)
    def __sub__(x,y):
        return divint(x.z*y.n-y.z*x.n)/divint(x.n*y.n)
    def __mul__(x,y):
        return divint(x.z*y.z)/divint(x.n*y.n)
    def __rmul__(x,y):
        return divint(x.z*y)/divint(x.n)
    def __pow__(x,p):
        return divint(x.z**p)/divint(x.n**p)
    def __div__(x,y):
        return divint(x.z*y.n)/divint(x.n*y.z)
    def __repr__(self):
        return str(self.z) + '/' + str(self.n)
    def __mini(self,x,y):
        if len(x) <= len(y):
            return x[:],x,y
        else: return y[:],y,x
    def __reduce(self,a,b):
        a = primetest().fak(a)
        b = primetest().fak(b)
        x=self.__mini(a,b)
        for i in x[0]:
            if i in x[2]:
                x[1].remove(i)
                x[2].remove(i)
        x[1].append(1)
        x[2].append(1)
        return reduce(int.__mul__,a),reduce(int.__mul__,b)

class divint(int):
    def __init__(self,x):
        self.z = x
        self.n = 1
    def __div__(x,y):
            return ratio(x,y)

>>> r1 = ratio(2,4)
>>> r2 = ratio(8,12)
>>> r3 = divint(12)/divint(15)
>>> r1,r2,r3
(1/2, 2/3, 4/5)
>>> r1+r2+r3
59/30
>>> r1-r2+r3
19/30
>>> r2+r1*r3
16/15
>>> r1*r2**2
2/9
>>>
Zuletzt geändert von Qubit am Sonntag 2. November 2008, 23:21, insgesamt 1-mal geändert.
Benutzeravatar
Leonidas
Administrator
Beiträge: 16023
Registriert: Freitag 20. Juni 2003, 16:30
Kontaktdaten:

Beitragvon Leonidas » Sonntag 2. November 2008, 15:11

Python 2.6 bringt einen Bruch-Datentyp in der Stdlib mit.
My god, it's full of CARs! | Leonidasvoice vs Modvoice
BlackJack

Beitragvon BlackJack » Sonntag 2. November 2008, 17:08

@Qubit: Beim Anwenden sieht man ganz gut, das die Klasse überflüssig ist. Wenn man ein Exemplar nur erzeugt, um darauf eine Methode auf zu rufen und es dann einfach wieder verwirft, ist das ein deutlicher Hinweis, dass die Klasse überflüssig ist.
Benutzeravatar
HWK
User
Beiträge: 1295
Registriert: Mittwoch 7. Juni 2006, 20:44

Beitragvon HWK » Sonntag 2. November 2008, 17:44

@Qubit: Die Zeiten sind ja erst einmal beeindruckend. Sie hängen aber doch ganz erheblich von den Ausgangszahlen ab. Bei z.B. '3165771122567177' sieht das Ergebnis deutlich schlechter aus.
Ich habe als Vergleich einmal die Pollard-Rho-Methode implementiert. Die Ergebnisse sind natürlich auch hier von den Ausgangszahlen abhängig, insgesamt aber sicher nicht schlechter. Die Implementierung ist aber wohl verständlicher.

Code: Alles auswählen

from math import sqrt
from random import randint

def ggt(a, b):
    a, b = abs(a), abs(b)
    while b:
        a, b = b, a % b
    return a

def rho(n, x0=None):
    if x0 is None:
        x0 = randint(2, int(sqrt(n)))
    while True:
        d = randint(1, 1000)
        x = y = x0
        nn = n
        result = []
        while nn > 1:
            while True:
                x = (x * x + d) % nn
                t = (y * y + d) % nn
                y = (t * t + d) % nn
                g = ggt(abs(y - x), nn)
                if g > 1:
                    result.append(g)
                    nn /= g
                    break
        if result[0] != n:
            result.sort()
            return result

if __name__ == '__main__':
    from time import clock
    for n in (34823482348238947, 348234823482389479, 3165771122567177):
        t = clock()
        print rho(n)
        print clock() - t
Ergebnis:[code=][11L, 29576089L, 107038193L]
0.0694772659654
[413243L, 842687773253L]
23.8080616391
[29576089L, 107038193L]
0.258115486745[/code]Ergebnis bei Deiner Variante mit denselben Zahlen:[code=][11, 29576089, 107038193L]
('19.252', 'sec')
[413243, 842687773253L]
('1.840', 'sec')
[29576089, 107038193L]
('19.322', 'sec')[/code]MfG
HWK
Benutzeravatar
HWK
User
Beiträge: 1295
Registriert: Mittwoch 7. Juni 2006, 20:44

Beitragvon HWK » Sonntag 2. November 2008, 19:59

Das war wohl etwas voreilig. Die Implementierung funktioniert so leider nicht. Ich habe es jetzt etwas angepasst. Dazu war allerdings eine Überprüfung auf Primzahleigenschaften nötig. Dazu habe ich den Miller-Rabin-Test verwendet, auch wenn er eine gewisse Fehlerquote lässt. Mit 9.0949470177292824e-013 ist diese aber sicher zu verschmerzen. Die Geschwindigkeit wurde insgesamt durch die Änderungen erheblich gesteigert.
Hier der Code:

Code: Alles auswählen

from math import sqrt
from random import randint

TIMES = 20

def ggt(a, b):
    a, b = abs(a), abs(b)
    while b:
        a, b = b, a % b
    return a

def is_prime(n):
    if n in [1, 2, 3]:
        return True
    s = 0
    t = n - 1
    while not (t % 2):
        s += 1
        t //= 2
    for _ in xrange(TIMES):
        b = randint(2, n-2)
        z = pow(b, t, n)
        if z != 1 and z != n-1:
            for _ in xrange(TIMES):
                z = pow(z, 2, n)
                if z == n-1:
                    break
            else:
                return False
    return True

def rho(n, x0=None):
    if n == 1:
        return []
    if n == 4:
        return [2, 2]
    if is_prime(n):
        return [n]
    if x0 is None:
        x0 = randint(2, int(sqrt(n)))
    while True:
        x = y = x0
        d = randint(1, 1000)
        while True:
            x = (x * x + d) % n
            t = (y * y + d) % n
            y = (t * t + d) % n
            g = ggt(abs(y - x), n)
            if g > 1:
                if g != n:
                    return rho(n / g) + rho(g)
                break

if __name__ == '__main__':
    from time import clock
    for n in (34823482348238947, 348234823482389479, 3165771122567177):
##     for _ in xrange(5):
##         n = randint(10000000, 1000000000000000000)
        print n
        t = clock()
        print rho(n)
        print clock() - t
Und das Ergebnis:[code=]34823482348238947
[107038193L, 29576089L, 11L]
0.0919432497706
348234823482389479
[842687773253L, 413243L]
0.0252710889233
3165771122567177
[107038193L, 29576089L]
0.191376557635[/code]MfG
HWK
Qubit
User
Beiträge: 75
Registriert: Dienstag 7. Oktober 2008, 09:07

Beitragvon Qubit » Sonntag 2. November 2008, 23:46

HWK hat geschrieben:@Qubit: Die Zeiten sind ja erst einmal beeindruckend. Sie hängen aber doch ganz erheblich von den Ausgangszahlen ab. Bei z.B. '3165771122567177' sieht das Ergebnis deutlich schlechter aus.

Hallo HWK,
bin mir schon bewusst, dass der von mir gezeigte Algorithmus i.a. gegen probabilistische Tests, AKS, etc. schlechter aussieht. Zumindestens ist mein Algorithmus aber 100% sicher ;-)
Beschäftige mich mit Primzahlen nur nebensächlich, bin da kein Profi noch ist das ein Hobby von mir. Habe dabei aber ein paar Zusammenhänge erkannt, die mE so noch nicht publik sind (aber wie gesagt, bin auf dem Gebiet kein Profi) und im Algorithmus umgesetzt wurden. Dieser implementiert allerdings noch nicht alle diese Zusammenhänge, so dass ich ihn in absehbarer Zeit in dieser Richtung noch verbessern werde. Dann schauen wir weiter :wink:

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder