Seite 1 von 2
Verfasst: Mittwoch 29. Oktober 2008, 17:46
von nuss
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)
Verfasst: Freitag 31. Oktober 2008, 17:03
von name
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
Verfasst: Samstag 1. November 2008, 02:27
von Qubit
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)
Verfasst: Samstag 1. November 2008, 19:30
von name
Scheint zu funktionieren, schaut aber eher, naja, nicht sehr uebersichtlich aus. Koenntest du den Algorithmus bitte einfach so, bzw mit Pseudocode erklaeren?
Verfasst: Sonntag 2. November 2008, 00:38
von Qubit
name hat geschrieben:
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
>>>
Verfasst: Sonntag 2. November 2008, 02:07
von audax
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.
Verfasst: Sonntag 2. November 2008, 11:37
von name
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.
Verfasst: Sonntag 2. November 2008, 15:08
von Qubit
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
>>>
Verfasst: Sonntag 2. November 2008, 15:11
von Leonidas
Python 2.6 bringt einen Bruch-Datentyp in der Stdlib mit.
Verfasst: Sonntag 2. November 2008, 17:08
von BlackJack
@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.
Verfasst: Sonntag 2. November 2008, 17:44
von HWK
@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: Alles auswählen
[11L, 29576089L, 107038193L]
0.0694772659654
[413243L, 842687773253L]
23.8080616391
[29576089L, 107038193L]
0.258115486745
Ergebnis bei Deiner Variante mit denselben Zahlen:
Code: Alles auswählen
[11, 29576089, 107038193L]
('19.252', 'sec')
[413243, 842687773253L]
('1.840', 'sec')
[29576089, 107038193L]
('19.322', 'sec')
MfG
HWK
Verfasst: Sonntag 2. November 2008, 19:59
von HWK
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: Alles auswählen
34823482348238947
[107038193L, 29576089L, 11L]
0.0919432497706
348234823482389479
[842687773253L, 413243L]
0.0252710889233
3165771122567177
[107038193L, 29576089L]
0.191376557635
MfG
HWK
Verfasst: Sonntag 2. November 2008, 23:46
von Qubit
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