Ein np.array vorne mit 1en befüllen, hinten zuschneiden

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Es ist nicht so, daß ich das nicht hinkriege, aber die Performance ist lausig.
Es geht um ein Primzahlensieb.

1, 1, 1, 1, 1, 1, [0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1]
Ergebnis sollte so aussehen:
1, 1, 1, 1, 1, [1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1], 1
oder so:
1, 1, 1, 1, [1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1], 1, 1

Das Sieb ist hier kurz dargestellt (Size = 4 x 7 = 28). In wirklichkeit ist es so um die 20 Mio´s lang.

Jemand ne Idee ? Gruß Kalli
__deets__
User
Beiträge: 14539
Registriert: Mittwoch 14. Oktober 2015, 14:29

Ich verstehe nicht, was du da erreichen willst. numpy ist nicht fuer eine solche Datenstruktur geeignet, bei der ein Element nicht eine Zahl, sondern wieder eine Liste ist. Und wozu genau soll das gut sein?
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Oben ist das Sieb für die 7. Mein np.Array hat aber noch mehr Dimensionen (11, 13, 17 .... usw). Zum Schluß werden die Dimensionen mit fast_mul multipliziert und erzeugt so das Endsieb.
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1] (Die Dimension für die 11 zu filtern).
_______________<---------________________________________________________________ <---------
_________________n______________________________________________________________n
Concenate, stack, kopieren und append ist alles zu langsam.
Benutzeravatar
__blackjack__
User
Beiträge: 13100
Registriert: Samstag 2. Juni 2018, 10:21
Wohnort: 127.0.0.1
Kontaktdaten:

@Karl-Heinz Hofmann: Ich verstehe auch nicht was das werden soll.
„All religions are the same: religion is basically guilt, with different holidays.” — Cathy Ladman
Sirius3
User
Beiträge: 17749
Registriert: Sonntag 21. Oktober 2012, 17:20

@Karl-Heinz Hofmann: man hat nur einmal ein Sieb, das man nach und nach füllt. Sonst weiß man ja gar nicht, welche Primzahl als nächstes ausgesiebt werden soll.

Zeig doch mal Deinen Code, dann versteht man vielleicht besser, was Du machst.
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

for y in range (0,len(LP)):
z = 0
while XList[z] % LP[y] != 0:
z += 1
nkArray *= UrArray [y] [LP[y] - z : LP[y] - z + Sieblaenge]

Es geht eigentlich nur um zwei sehl lange (10 Mio (Sieblaenge) lange) Arrays mit 0 en und 1en die multipliziert werden. D.h. wenn eine 0 aus UrArray auf eine 1 aus nkArray trifft, so soll aus der 1 eine 0 in nkArray werden.
Eigentlich funzt der Code korrekt, nur es geht mir zu langsam.
Sirius3
User
Beiträge: 17749
Registriert: Sonntag 21. Oktober 2012, 17:20

Bei mir „funzt” da gar nichts, weil LP, XList, nkArray und UrArray oder Sieblaenge nicht definiert sind.
Ich verstehe auch immer noch nicht, wie Du daraus ein Sieb bauen willst. Nach welchem Algorithmus arbeitest Du?
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Code: Alles auswählen

# -*- coding: utf-8 -*-
"""
Created on Wed May  8 18:33:37 2019

@author: Kalli Hofmann
"""

from timeit import default_timer as timer
import numpy

def array_attributes(a):
   for attr in ('ndim', 'size', 'itemsize', 'dtype', 'shape', 'strides'):
        print('{:8s}: {}'.format(attr, getattr(a, attr)))


def printArray(aArray):
  print(str(aArray[0,0:15])[0:-1] + "  ...  " + str(aArray[0,-15:])[1:] + " Size: " + str(aArray.size))
 

a = numpy.full((1, 100000000), 1, numpy.float)
for x in range (0,a.size):
    if x % 3 == 0 :
        a[0,x] = 0
        
b = numpy.full((1, 100000000), 1, numpy.float)
for x in range (0,b.size):
    if x % 5 == 0 :
        b[0,x] = 0

c = numpy.full((1, 100000000), 1, numpy.float)


print("vor der Multiplikation:")
printArray(a)
printArray(b)
printArray(c)
print()

start = timer()
c = a * b
duration = timer() - start

print("nach der Multiplikation:")
printArray(a)
printArray(b)
printArray(c)
print()
array_attributes(c)
print()
print("Sekunden gebraucht: " + str(duration))
Zuletzt geändert von Karl-Heinz Hofmann am Dienstag 19. November 2019, 11:21, insgesamt 1-mal geändert.
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Hallo
Ich habs mal stark vereinfacht und sofort lauffähig gemacht. Jetzt zu meinem Problem: Das " c = a * b " geht mir zu langsam.
0,45 sek. zeigt er mir an. Für das Herstellen von a und b braucht er sehr viel länger. Ist aber unherheblich, da a und b nur einmal am Anfang
hergestellt werden und für den weiteren Gebrauch unverändert zur Verfügung stehen und immer wieder verwendet werden.
Gruß Kalli
Benutzeravatar
sparrow
User
Beiträge: 4193
Registriert: Freitag 17. April 2009, 10:28

Ohne zu verstehen warum du das tust (es erscheint mir wenig sinnvoll im Bezug auf deinen ersten Post) und ohne auch nur die geringste Ahnung von numpy zu haben:

Warum baust du Arrays aus Fließkommazahlen, wenn du dort nur Ganzzahlen hinein steckst? Der Unterschied ist dir bekannt?
Benutzeravatar
__blackjack__
User
Beiträge: 13100
Registriert: Samstag 2. Juni 2018, 10:21
Wohnort: 127.0.0.1
Kontaktdaten:

@Karl-Heinz Hofmann: Pro Multiplikation also nur etwas weniger als 4½ Nanosekunden. Warum/wofür ist das zu langsam?

Warum `float` für Werte die nur 1 oder 0 sein können? Mit `bool` braucht das bei mir nur 0.048534920962993056 Sekunden.

Warum ein `shape` das zwei Indizes für den Zugriff nötig macht, wo die Daten eigentlich nur eindimensional sind?

`c` wird unnötigerweise definiert weil der zuerst zugewiesene Wert nie verwendet wird.

Im Gegensatz zur Geschwindigkeit bei der Multiplikation kann man bei der Erstellung von `a` und `b` tatsächlich etwas machen, denn das ist massiv ineffizient.

Das zusammengestückel in der `print_array()`-Funktion (so sollte der Name nach Konvention geschrieben werden) ist unschön. Das ist eher BASIC als Python. Zudem ist es nicht schön/robust an der Zeichenkettenrepräsentation von Objekten rumzuschnibbeln.

Zwischenstand:

Code: Alles auswählen

#!/usr/bin/env python3
"""
Created on Wed May  8 18:33:37 2019

@author: Kalli Hofmann
"""
from timeit import default_timer as timer

import numpy


def array_attributes(array):
    for attr in ["ndim", "size", "itemsize", "dtype", "shape", "strides"]:
        print("{:8s}: {}".format(attr, getattr(array, attr)))


def print_array(array):
    with numpy.printoptions(edgeitems=15):
        print(f"{array} Size: {array.size}")


def main():
    size = 100_000_000
    a = numpy.full(size, True, numpy.bool)
    a[0::3] = False

    b = numpy.full(size, True, numpy.bool)
    b[0::5] = False

    print("vor der ”Multiplikation”:")
    print_array(a)
    print_array(b)
    print()

    start = timer()
    c = a & b
    duration = timer() - start

    print("nach der ”Multiplikation”:")
    print_array(a)
    print_array(b)
    print_array(c)
    print()
    array_attributes(c)
    print()
    print("Sekunden gebraucht: " + str(duration))


if __name__ == "__main__":
    main()
Und jetzt stelle ich mir natürlich die Frage ob/warum man hier überhaupt die Multiplikation oder die Und-Verknüpfung von zwei Arrays braucht die ein drittes Array erstellt‽ Wenn ich mit Numpy das Sieb des Erathostenes implementieren wollte, bräuchte ich genau *ein* Array, eine Schleife und die Zuweisung mit „slice“-Notation wie beim setzen der Werte mit durch 3 oder 5 teilbaren Indizes auf `False` im Beispiel oben.
„All religions are the same: religion is basically guilt, with different holidays.” — Cathy Ladman
Sirius3
User
Beiträge: 17749
Registriert: Sonntag 21. Oktober 2012, 17:20

Für was brauchst Du denn bei einem "Sieb" a und b mehrfach?
`c` wird in Zeile 30 definiert aber nie benutzt, kann also weg.
Statt mit for-Schleifen solltest Du die Fähigkeiten von numpy zur Vektorrechnung nutzen:

Code: Alles auswählen

a = (numpy.arange(100000000) % 3) == 0
Wenn Du statt mit Floats mit Bools rechnest, geht die Multiplikation auch deutlich schneller.

Trotzdem hast Du immer noch kein Sieb. Beim Sieb braucht man nämlich KEINE Division oder Multiplikation. Deshalb ist es so schnell.
Mal als Ansatz:

Code: Alles auswählen

a = numpy.ones(100000000, dtype='u1')
a[2*2::2] = 0
Benutzeravatar
__blackjack__
User
Beiträge: 13100
Registriert: Samstag 2. Juni 2018, 10:21
Wohnort: 127.0.0.1
Kontaktdaten:

@Karl-Heinz Hofmann: Nur mal um zu zeigen das Python/Numpy da gar nicht so schlecht ist, hier mal die Multiplikation der 100 Millionen 64-Bit-Gleitkommazahlen in C:

Code: Alles auswählen

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

#define SIZE 	100000000
#define timespec_to_seconds(ts)	(ts.tv_sec + ts.tv_nsec / 1e9)

void print_array(double *a, uint32_t size) {
	uint32_t i;
	
	putchar('[');
	for (i = 0; i < 15; ++i) printf("%.0f ", a[i]);
	printf("... ");
	for (i = 0; i < 15; ++i) printf("%.0f ", a[size - i - 1]);
	printf("] Size: %"PRIu32"\n", size);
}

int main(void)
{
	double *a, *b, *c;
	uint32_t i;
	struct timespec start, end;
	
	a = malloc(SIZE * sizeof(*a));
	b = malloc(SIZE * sizeof(*b));
	c = malloc(SIZE * sizeof(*c));
	
	for (i=0; i < SIZE; ++i) a[i] = b[i] = 1;
	for (i=0; i < SIZE; i += 3) a[i] = 0;
	for (i=0; i < SIZE; i += 5) b[i] = 0;
	
	clock_gettime(CLOCK_MONOTONIC, &start);
	for (i=0; i < SIZE; ++i) c[i] = a[i] * b[i];
	clock_gettime(CLOCK_MONOTONIC, &end);
	
	print_array(c, SIZE);
	printf("Sekunden für die Multiplikation: %f\n",
		timespec_to_seconds(end) - timespec_to_seconds(start));
	
	free(c);
	free(b);
	free(a);
	return 0;
}
Braucht 0.321240 Sekunden bei mir, gegenüber den 0.418104 Sekunden von Deinem Python-Code.
„All religions are the same: religion is basically guilt, with different holidays.” — Cathy Ladman
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Hallo Leute
Vielen, vielen Dank. Da sind ein Haufen guter Ideen drin. Das mit "Bool" bin ich heute Nachmittag selbst draufgekommen. Zuerst auf "int8" und dann mit "Bool" probiert. :-) Ihr fragt euch warum ich das Sieb nicht von vornherein konstruiere?
Also es geht um Mersenne Primzahlen, speziell um kleine Faktoren zu suchen. Dabei untersuche ich viele M´s hintereinander. Wenn ich z.B https://www.mersenne.ca/exponent/2917107953 = 2^2917107953 - 1 untersuche, kommen nur
folgende Faktoren in Frage: 2k * Exponent + 1 Für 2^2917107953 - 1 ist das so mit 5834215907. Der nächste mögliche 11668431813 . Diese Folge kann man auch sieben und folgt einer gewissen Regelmäßigkeit. Dennoch muß ich das Sieb bei jedem neuen Exponenten neu ausrichten. Ich hab unten nochmal das Ergebnis von heute Nachmittag gepostet und das ganze um zwei "Faktor Arrays" erweitert. Entschuldigt meinen schlechten Programmierstil. Python bringe ich mir erst seit einem halben Jahr selbst bei. Eigentlich bin ich Pascal-/Delphianer.
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Code: Alles auswählen

# -*- coding: utf-8 -*-
"""
Created on Wed May  8 18:33:37 2019

@author: Kalli Hofmann
"""

from timeit import default_timer as timer
import numpy


def array_attributes(a):
   for attr in ('ndim', 'size', 'itemsize', 'dtype', 'shape', 'strides'):
        print('{:8s}: {}'.format(attr, getattr(a, attr)))


def printArray(aArray):
    Vorne = str(aArray[0:20])[0:-1]
    Hinten = str(aArray[-20:])[1:]
    Vorne = Vorne.replace(" True","1").replace("False","0").replace("\n","")
    Hinten = Hinten.replace(" True","1").replace("False","0").replace("\n","")
    print(Vorne + "  ...  " + Hinten + " Size: " + str(aArray.size))
 

Sieblaenge = 20000000

my_array = numpy.full((4, Sieblaenge), 1, numpy.bool)
for x in range (0,Sieblaenge):
    if x %  3 == 0:  my_array[0,x] = 0
    if x %  5 == 0:  my_array[1,x] = 0
    if x %  7 == 0:  my_array[2,x] = 0 
    if x % 11 == 0:  my_array[3,x] = 0   
res = numpy.full(( Sieblaenge), 1, numpy.bool)


print("vor der Multiplikation:")
printArray(my_array[0,:])
printArray(my_array[1,:])
printArray(my_array[2,:])
printArray(my_array[3,:])
printArray(res)
print()

start = timer()
res = my_array[0,:] * my_array[1,:] * my_array[2,:] * my_array[3,:]
duration = timer() - start

print("nach der Multiplikation:")
printArray(my_array[0,:])
printArray(my_array[1,:])
printArray(my_array[2,:])
printArray(my_array[3,:])
printArray(res)
print()
array_attributes(res)
print()
print("Sekunden gebraucht: " + str(duration))

Zuletzt geändert von Karl-Heinz Hofmann am Dienstag 19. November 2019, 16:28, insgesamt 1-mal geändert.
Sirius3
User
Beiträge: 17749
Registriert: Sonntag 21. Oktober 2012, 17:20

Egal wie man das dreht, Dein Code ist um Größenordnungen langsamer als ein richtig implementiertes Sieb. Also implementiert erst ein normales Sieb bevor du hier Spezialfälle anschaust.
Benutzeravatar
__blackjack__
User
Beiträge: 13100
Registriert: Samstag 2. Juni 2018, 10:21
Wohnort: 127.0.0.1
Kontaktdaten:

@Karl-Heinz Hofmann: Zu Deinem letzten Quelltext: `my_*` ist ein unsinnige Vorsilbe wenn es nicht auch `our_*` oder `their_*` gibt.

Irgendwie hast Du vergessen die durch 2 teilbaren Zahlen auszusieben‽

Du hast immer noch die Schleife für das sieben in Python statt Zuweisung an einen ”slice”. Theoretisch mag die *eine* Schleife schneller sein als eine interne Schleife vom Numpy-Array, aber in der Praxis ist Python einfach *deutlich* langsamer.

Was vorher `c` war ist jetzt `res` — die erste Zuweisung ist einfach nur Verschwendung. Du legst da ein Array mit 20 Millionen Werten an die nie (sinnvoll) genutzt werden. Da wird Speicher mit 20 Millionen True-Werten angelegt, um dann die ersten und letzten 20 davon auszugeben. Und dann wird dieses Objekt solange unnütz im Speicher gehalten bis das Ergebnis in einem anderen Array berechnet wurde und es davon ersetzt wird.

Die Indexwerte 0 bis 3 sollte man nicht immer explizit hinschreiben. Dafür gibt es Schleifen.

Und das mit dem Multiplizieren sollte eigentlich ein logisches ”und” sein.

Code: Alles auswählen

#!/usr/bin/env python3
"""
Created on Wed May  8 18:33:37 2019

@author: Kalli Hofmann
"""
from timeit import default_timer as timer

import numpy


def print_array(array):
    with numpy.printoptions(edgeitems=15):
        print(f"{array} Size: {array.size}")


def main():
    size = 20_000_000

    primes = [2, 3, 5, 7, 11]
    array = numpy.full((len(primes), size), True, bool)
    for i, prime in enumerate(primes):
        array[i, ::prime] = False

    print("vor der Multiplikation:")
    for row in array:
        print_array(row)
    print()

    start = timer()
    result = numpy.logical_and.reduce(array, 0)
    duration = timer() - start

    print("nach der Multiplikation:")
    for row in array:
        print_array(row)
    print_array(result)
    print()
    print()
    print("Sekunden gebraucht: " + str(duration))


if __name__ == "__main__":
    main()
Wobei das hier das gleiche Problem hat wie Dein Code: es werden auch die fest im Code verankerten Primzahlen selbst mit ausgesiebt. Das möchte man vielleicht korrigieren.

Um den Punkt von Sirius3 noch mal als Code deutlich zu machen: Ein ganz normales, simples Sieb findet mit den folgenden einfachen Zeilen auf meinem Rechner die 5.761.455 Primzahlen zwischen 0 und 100.000.000 und braucht dafür etwa 6 Sekunden:

Code: Alles auswählen

#!/usr/bin/env python3
import numpy


def print_array(array):
    with numpy.printoptions(edgeitems=15):
        print(f"{array} Size: {array.size}")


def main():
    size = 100_000_000
    is_prime = numpy.full(size, True, bool)
    is_prime[0:2] = False
    for i in range(2, int(numpy.sqrt(size)) + 1):
        is_prime[i * 2 :: i] = False

    primes = numpy.nonzero(is_prime)[0]
    print_array(primes)


if __name__ == "__main__":
    main()
Ausgabe:

Code: Alles auswählen

[       2        3        5        7       11       13       17       19
       23       29       31       37       41       43       47 ...
 99999643 99999677 99999703 99999721 99999773 99999787 99999821 99999827
 99999839 99999847 99999931 99999941 99999959 99999971 99999989] Size: 5761455
„All religions are the same: religion is basically guilt, with different holidays.” — Cathy Ladman
Benutzeravatar
__blackjack__
User
Beiträge: 13100
Registriert: Samstag 2. Juni 2018, 10:21
Wohnort: 127.0.0.1
Kontaktdaten:

Hier mal ein Pascal-Programm das zusätzlich zum Sieb für den gleichen Wertebereich auch noch alle Mersenne-Zahlen aus dem Bereich ausgibt und ob sie prim sind oder nicht:

Code: Alles auswählen

program Sieve;

const MaxNumber = 100000000;

var
  isPrime: array[0..MaxNumber] of Boolean;
  primes: array of LongInt;

procedure DoSieve;
var
  i, j: LongInt;
begin
  for i := Low(isPrime) to High(isPrime) do isPrime[i] := TRUE;
  isPrime[0] := FALSE;
  isPrime[1] := FALSE;
  for i := Low(isPrime) to Round(Sqrt(High(isPrime))) do begin
    if isPrime[i] then begin
      j := i + i;
      while j <= High(isPrime) do begin
        isPrime[j] := FALSE;
        Inc(j, i);
      end;
    end;
  end;
end;

function CountPrimes: LongInt;
var
  n, i: LongInt;
begin
  n := 0;
  for i := Low(isPrime) to High(isPrime) do begin
    if isPrime[i] then Inc(n);
  end;
  CountPrimes := n;
end;

procedure CollectPrimes;
var
  i, j: LongInt;
begin
  SetLength(primes, CountPrimes());
  j := 0;
  for i := 0 to High(isPrime) do begin
    if isPrime[i] then begin
      primes[j] := i;
      Inc(j);
    end;
  end;
end;

procedure PrintPrimes;
var
  i: LongInt;
begin
  for i := 0 to 14 do Write(primes[i], ' ');
  Write('... ');
  for i := High(primes) - 15 to High(primes) do Write(primes[i], ' ');
  WriteLn(' Size: ', High(primes) + 1);
end;

procedure PrintMersenneNumbers;
var
  m: Longint;
begin
  WriteLn('Mersenne no. | is prime');
  m := 1;
  while m < High(isPrime) do begin
    WriteLn(m:12, ' | ', isPrime[m]);
    m := (m Shl 1) Or 1;
  end;
end;

begin
  DoSieve;
  CollectPrimes;
  PrintPrimes;
  PrintMersenneNumbers;
end.
Läuft in 1,34 Sekunden durch bei mir und gibt folgendes aus:

Code: Alles auswählen

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 ... 99999623 99999643 99999677 99999703 99999721 99999773 99999787 99999821 99999827 99999839 99999847 99999931 99999941 99999959 99999971 99999989  Size: 5761455
Mersenne no. | is prime
           1 | FALSE
           3 | TRUE
           7 | TRUE
          15 | FALSE
          31 | TRUE
          63 | FALSE
         127 | TRUE
         255 | FALSE
         511 | FALSE
        1023 | FALSE
        2047 | FALSE
        4095 | FALSE
        8191 | TRUE
       16383 | FALSE
       32767 | FALSE
       65535 | FALSE
      131071 | TRUE
      262143 | FALSE
      524287 | TRUE
     1048575 | FALSE
     2097151 | FALSE
     4194303 | FALSE
     8388607 | FALSE
    16777215 | FALSE
    33554431 | FALSE
    67108863 | FALSE
„All religions are the same: religion is basically guilt, with different holidays.” — Cathy Ladman
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Code: Alles auswählen

# -*- coding: utf-8 -*-
"""
Created on Wed May  8 18:33:37 2019

@author: Kalli Hofmann
"""
from urllib.request import urlopen, URLError
from math import log
from numbthy import factors, isprime
import numpy
import re
import time
import json
import os, sys
import webbrowser



def extract_charset(response): #══ Für Auslesen einer Page ═══════════════════
    m = re.search('charset=(\S+)',response.getheader('content-type') or '') 
    return m and m.group(1) or 'ISO-8859-1'                                 
                                                                            
   
def array_attributes(a): #══ Anzeigen der Arrayattribute ═════════════════════
   for attr in ('ndim', 'size', 'itemsize', 'dtype', 'shape', 'strides'):
        print('{:8s}: {}'.format(attr, getattr(a, attr)))

        
def print_array(aArray): #══ Anzeigen eines Arrays, so wie ich es haben will ══
    Vorne = str(aArray[0:20])[0:-1]
    Hinten = str(aArray[-20:])[1:]
    Vorne = Vorne.replace(" True","1").replace("False","0").replace("\n","")
    Hinten = Hinten.replace(" True","1").replace("False","0").replace("\n","")
    print(Vorne + "  ...  " + Hinten + " Size: " + str(aArray.size))
    

def kontrolle(p): #══ Kontrolle eines Wertes ═════════════════════════════════
    if (p % 3 and p % 5 and p % 7 and p % 13) == 0:
                 print(factors(p))                                   


def zusatz_aussiebung(Liste,pmin,pmax): #══ Zusatzaussiebung einer Liste ═════
    for LP in pList[pMin:pMax]: # 5 = 13
       triple = 3 * LP
       doubleminus = Sieblaenge - 2 * LP
       x = 0
       while XList[x] % LP != 0:  x += 1
       while XList[x] % 3 != 0:
           Liste[x] = 0
           x += LP
       while x < doubleminus:
              Liste[x + LP] = 0
              Liste[x + LP + LP] = 0
              x += triple


def mach_die_pList(p): #══ Macht eine kleine Primzahlenliste ════════════════════
    print("Stelle pList her.")
    pList = [p]
    for i in range (3,150001):
        if isprime(i): pList.append(i)
    print("pList hergestellt. pList hat " + str(len(pList)) +
          " Elemente. Größtes Element ist: " + str(pList[len(pList)-1] ))
    print("Wert von pMin = " + str(pList[pMin]) +
          "  Wert von pMax = " + str(pList[pMax]))
    return pList
              

#█████████████ Hauptparameter ██████████████████████████████████████████████████
                                                                             #██
Anfang       = 6007000000                                                    #██
WiederAnfang = 6007000000                                                    #██
Ende         = 6008000000                                                    #██
FaktorAnfang = 57    # in 2 ↑ er Potenz                                      #██
FaktorEnde   = 58    # in 2 ↑ er Potenz                                      #██
pMin         = 5     # 5.tes Element (=13) in der pList (nicht verändern !!) #██
pMax         = 1000  # Wert zwischen 900 und 1100 ist gut                    #██
                                                                             #██
#███████████████████████████████████████████████████████████████████████████████

kA           = (pow(2,FaktorAnfang))//Ende - 2  # - 2 ist Headroom
if kA % 2   != 0: kA -= 1
if kA        < 2: kA = 2
kA_Halb      = kA // 2
kE           = pow(2,FaktorEnde) // Anfang + 2  # + 2 ist Headroom
if kE % 2   != 0: kE += 1
Sieblaenge   = 2
while Sieblaenge % 4620 != 0:
         kE += 2
         Sieblaenge  = (kE - kA) // 2
Siebmultiplikator = Sieblaenge // 4620
print("Sieblaenge: " + str(Sieblaenge))
print("Siebmultiplikator: " + str(Siebmultiplikator))
ErgebnisSpeicherPfad = (os.path.basename(sys.argv[0])[:-3] + "_" + 
                       str(Anfang) + "_" + str(Ende) + "_" +  
                       str(FaktorAnfang) + "_" + str(FaktorEnde) + ".txt")
print("Speichere Ergebnisse unter " + str(ErgebnisSpeicherPfad) + " ab.")
content              = ''
pList                = mach_die_pList(2)


for Exponent in range (WiederAnfang,Ende): #══ Hier gehts los ══════════════════
 if isprime(Exponent):
  ZeitStart   = time.time()
  DExponent   = Exponent * 2
  FactorCount = 0
 
  XList = [] #══ X, Y Listen herstellen ════════════════════════════════════════
  for x in range(0,pList[pMax * 3]):
       XList.append(((kA_Halb + x) * DExponent) + 1)
  YList = [0]
  for y in range(1,Siebmultiplikator):
       YList.append(y * 4620 * DExponent)

  #══ Grundsieb herstellen ═════════════════════════════════════════════════════

  Class  = numpy.full(4620, True, numpy.bool)
  kArray = numpy.full(Sieblaenge, True, numpy.bool)
  if (kA + 4)  % 8 == 0: kArray[0::4] = False 
  if (kA + 6)  % 8 == 0: kArray[1::4] = False 
  if (kA + 8)  % 8 == 0: kArray[2::4] = False
  if (kA + 10) % 8 == 0: kArray[3::4] = False
  for pf in [3,5,7,11]:
    for x in range(0,pf):
       if (((Exponent * kA) + 1) + (DExponent * x)) % pf == 0:
           kArray[x::pf] = False
  Class[:] = kArray[0:4620]
  Zeit1   = time.time()
  
  #══ Zusätzliches Aussieben ═══════════════════════════════════════════════════

  #zusatz_aussiebung(kArray,5,pMax)
  for pf in pList[5:pMax]:
    for x in range(0,pf):
       if (((Exponent * kA) + 1) + (DExponent * x)) % pf == 0:
           kArray[x::pf] = False
             
  Zeit2   = time.time()
  kArray = numpy.reshape(kArray,(Siebmultiplikator,4620))

  #══ Das eigentliche Testen ═══════════════════════════════════════════════════
  
  for x in range(0,4620):
    if Class.item(x) == True:
      for y in range(0,Siebmultiplikator):
        if kArray.item(y,x) == True:  
          p = XList[x] + YList[y]
          if pow(2,Exponent,p) == 1:      # Hier wird letztendlich getestet
             if isprime(p):                           # Wird nachgeschaltet um zusammengestzte p auszuschließen ( Ist langsamer als pow)
                 k = (p - 1) // Exponent
                 FactorCount += 1
                 try:           # nachguggen ob der Faktor schon gefunden wurde und notiert ist
                   with urlopen('http://www.mersenne.ca/exponent/'+str(Exponent)) as response: data = response.read()
                   content = data.decode(extract_charset(response))
                 except URLError:
                            content = ''
                            print('URLError bei M' + str(Exponent))
                            time.sleep(1200)
                            with urlopen('http://www.mersenne.ca/exponent/'+str(Exponent)) as response: data = response.read()
                            content = data.decode(extract_charset(response))
                            exit
                 if content != '':
                    if not str(p) in content :
                      print ("Exponent: M" + str(Exponent) + "   Faktor No.: " + str(k) + " = " + str(p) + " (" + str(round(log(p,2),4)) + " bit) " + "ergibt Rest  " + "0")
                      f = open (ErgebnisSpeicherPfad, "a")
                      f.write("UID: Karl-Heinz_Hofmann/ALDI, M" + str(Exponent) + " has a factor: " + str(p) + "\n")
                      f.close()
                      
  Zeit3   = time.time()
  f = open (ErgebnisSpeicherPfad, "a") #══ Ergebnis in Datei schreiben ════════════                   
  if FactorCount == 0:
     f.write("UID: Karl-Heinz_Hofmann/ALDI, no factor for M" + str(Exponent) + " from 2^" + str(FaktorAnfang) + " to 2^" + str(FaktorEnde) + " [mfaktc 0.21 75bit_mul32_gs]" + "\n")
  if FactorCount == 1:
     f.write("UID: Karl-Heinz_Hofmann/ALDI, found 1 factor for M" + str(Exponent) + " from 2^" + str(FaktorAnfang) + " to 2^" + str(FaktorEnde) + " [mfaktc 0.21 75bit_mul32_gs]" + "\n")
  if FactorCount > 1:
     f.write("UID: Karl-Heinz_Hofmann/ALDI, found " + str(FactorCount) + " factors for M" + str(Exponent) + " from 2^" + str(FaktorAnfang) + " to 2^" + str(FaktorEnde) + " [mfaktc 0.21 75bit_mul32_gs]" + "\n")
  f.close()   
                 
  ZeitStop = time.time() #══ Ausgabe ══════════════════════════════════════════════
  
  print("Habe für M" + str(Exponent) + "%10.2f Sekunden gebraucht.  " % (ZeitStop - ZeitStart) + "( " + str(FactorCount) + " Faktoren gefunden )   k_min = " + str(kA_Halb) + "   k_max = " + str(kE // 2))
  #print(counter)
  #print("%13.4f  " % (Zeit1 - ZeitStart) )
  #print("%13.4f  " % (Zeit2 - Zeit1) )
  #print("%13.4f  " % (Zeit3 - Zeit2) )
 
Hab mal hier den verbesserten Gesamtcode gepostet. Unten kann man die Zeiten für die Abschnitte mal einschalten. Spielen kann man mal mit dem pMax Wert.
Ich glaube, wenn man das Sieben auf die Grafikkarte (Hab ne Nvidia 1070 und Cuda installiert) verlagern könnte, dann ging noch was. Auch die pow Abfrage ist
ne Engstelle. Gruß Kalli
Karl-Heinz Hofmann
User
Beiträge: 20
Registriert: Dienstag 12. November 2019, 15:34

Ach ja, 2 braucht man nicht zu sieben, da automatisch alle p´s wegen p = (2k * Exponent + 1) automatisch alle ungerade sind. Aber die 4 kann und
muß man sieben (hat mathematischen Grund), da (k+2) % 4 nie 0 sein darf. Also die k´s =[ 2, 6, 10 , 14, 18 ......] sind nicht möglich. Die Zahl 4620 ensteht durch
die Multiplikation von 3, 4, 5, 7, 11 Wenn man das Sieb zum Schluß mit reshape wieder faltet, kann man zum Schluß schön durch die Lücke schießen. Gruß Kalli
Antworten