Python rechnet uneffizient?

[snafu]

Stutzig macht mich ja, dass der GCC beim kompakten Code `INC` zum Erhöhen der Register-Wertes in der Schleife verwendet, dies aber in der performanten Variante über ein explizites `ADD` erledigt. Die GCC-Entwickler werden sicherlich ihre Gründe haben, warum sie das so implementiert haben, aber eigenartig finde ich es trotzdem.

[BlackJack]

@snafu: Der Grund ist ganz einfach: ``INC reg`` ist kompakter und ``ADD reg, 1`` ist auf den meisten x64-Prozessoren mindestens gleich schnell. Eine Ebene tiefer werden beide ziemlich wahrscheinlich zum gleichen, oder zumindest sehr ähnlichen Mikroinstruktionen. Ganz äquivalent sind die beiden von der Funktion her nicht, weil INC nicht alle Flags beeinflusst die von ADD beeinflusst werden.

Für Intel-Prozessoren generiert mir der GCC für den gezeigten Quelltext immer die ADD-Variante. Für AMD-Prozessoren dagegen INC. Das ist ja auch so ein Punkt den man bei handoptimiertem Assembler beachten muss: Verschiedene Prozessoren behandeln die gleichen Opcodes unterschiedlich. Und dann kommen noch Sachen wie Pipelining dazu, die dafür sorgen können, dass neben der Wahl des Opcodes auch die Reihenfolge der Abarbeitung eine Rolle spielen. Wenn man also Assembler von Hand schreibt, weil man glaubt man ist grundsätzlich besser als der Compiler, dann sollte man sich auch wirklich sicher sein, dass der handgeschriebene Code auch dann noch besser ist, wenn das C- oder C++-Gegenstück auf einem beliebigen Prozessor (der gleichen Architektur) speziell für *den* Prozessor optimiert wurde.

[DerTürke]

Hallo Jungs,

ich kanns einfach nicht sein lassen aber mein Algo ist nun mal schneller.
Diesmal habe ich wie von euch beanstandet das Sieb in seiner ursprünglichen Form in reinem C geschrieben und dennoch ist meine IsPrime Funktion schneller, tut mir echt leid aber das ist nun mal Tatsache ob ihr das hinnehmen wollt oder nicht

Code: Alles auswählen

#include<iostream>
#include<vector>
#include<stdlib.h>

using namespace std;

bool IsPrime(int pValue);

int main()
{	
	int n_Sieb[200000];
	int n_Primes[33860];
	int n_PrimesIndex = 1;	

	n_Sieb[0] = 2;
	n_Primes[0] = 2;
	for (int i = 1; i < 200000; i++)
	{
		n_Sieb[i] = i * 2 + 1;
	}

	for (int i = 1; i < 200000; i++)
	{
		if (n_Sieb[i] > 0)
		{
			n_Primes[n_PrimesIndex++] = n_Sieb[i];
			for (int j = i + 1; j < 200000; j++)
			{
				if (n_Sieb[j] % n_Sieb[i] == 0)
				{
					n_Sieb[j] = 0;
				}
			}
		}
	}
	
	/// Jetzt kommt mein Algorithmus
	n_PrimesIndex = 0;
	for (int i = 2; i < 400000; i++)
	{
		bool bIsPrime = IsPrime(i);				
		if (bIsPrime)
		{
			n_Primes[n_PrimesIndex++] = i;
		}		
	}
	return 0;
}


bool IsPrime(int pValue)
{
	bool retValue = true;
	int m_Divisor = 2;
	int m_Divident = pValue;
	div_t div_Result;
	div_Result.quot = m_Divident;
	while (m_Divisor < div_Result.quot)
	{
		div_Result = div(m_Divident, m_Divisor);
		if (div_Result.rem == 0)
		{
			retValue = false;
			break;
		}
		m_Divisor++;
	}
	return retValue;
}

und auch mit Code Optimierung sieht das compilat zwar nun viel besser als vorher aus aber dennoch ist meins kompakter:

Code: Alles auswählen

using namespace std;

bool IsPrime(int pValue);

int main()
{	
00007FF6D6001000  mov         qword ptr [rsp+8],rbx  
00007FF6D6001005  mov         qword ptr [rsp+10h],rsi  
00007FF6D600100A  push        rdi  
00007FF6D600100B  mov         eax,0E4630h  
00007FF6D6001010  call        __chkstk (07FF6D6001EA0h)  
00007FF6D6001015  sub         rsp,rax  
00007FF6D6001018  movdqa      xmm2,xmmword ptr [__xmm@00000003000000020000000100000000 (07FF6D6002220h)]  
	int n_Sieb[200000];
	int n_Primes[33860];
	int n_PrimesIndex = 1;	

	n_Sieb[0] = 2;
00007FF6D6001020  lea         rcx,[rsp+21134h]  
00007FF6D6001028  movdqa      xmm3,xmmword ptr [__xmm@00000001000000010000000100000001 (07FF6D6002210h)]  
	n_Primes[0] = 2;
	for (int i = 1; i < 200000; i++)
00007FF6D6001030  mov         edx,1  
00007FF6D6001035  mov         dword ptr [n_Sieb],2  
00007FF6D6001040  lea         eax,[rdx+4]  
00007FF6D6001043  movd        xmm0,edx  
00007FF6D6001047  pshufd      xmm0,xmm0,0  
00007FF6D600104C  lea         rcx,[rcx+20h]  
00007FF6D6001050  paddd       xmm0,xmm2  
00007FF6D6001054  movd        xmm1,eax  
00007FF6D6001058  pshufd      xmm1,xmm1,0  
00007FF6D600105D  paddd       xmm0,xmm0  
00007FF6D6001061  paddd       xmm1,xmm2  
	{
		n_Sieb[i] = i * 2 + 1;
00007FF6D6001065  paddd       xmm0,xmm3  
00007FF6D6001069  movdqu      xmmword ptr [rcx-20h],xmm0  
00007FF6D600106E  add         edx,8  
00007FF6D6001071  paddd       xmm1,xmm1  
00007FF6D6001075  paddd       xmm1,xmm3  
	{
		n_Sieb[i] = i * 2 + 1;
00007FF6D6001079  movdqu      xmmword ptr [rcx-10h],xmm1  
00007FF6D600107E  cmp         edx,30D39h  
00007FF6D6001084  jl          main+40h (07FF6D6001040h)  
	n_Primes[0] = 2;
	for (int i = 1; i < 200000; i++)
00007FF6D6001086  cmp         edx,30D40h  
00007FF6D600108C  jge         main+0C0h (07FF6D60010C0h)  
00007FF6D600108E  movsxd      rax,edx  
00007FF6D6001091  lea         rcx,[n_Sieb]  
00007FF6D6001099  lea         rcx,[rcx+rax*4]  
00007FF6D600109D  lea         eax,[rdx*2+1]  
00007FF6D60010A4  nop         dword ptr [rax]  
00007FF6D60010A8  nop         dword ptr [rax+rax]  
	{
		n_Sieb[i] = i * 2 + 1;
00007FF6D60010B0  mov         dword ptr [rcx],eax  
00007FF6D60010B2  lea         rcx,[rcx+4]  
00007FF6D60010B6  add         eax,2  
00007FF6D60010B9  cmp         eax,61A81h  
00007FF6D60010BE  jl          main+0B0h (07FF6D60010B0h)  
	}

	for (int i = 1; i < 200000; i++)
00007FF6D60010C0  mov         r9d,2  
00007FF6D60010C6  lea         r8,[rsp+21134h]  
00007FF6D60010CE  lea         r10,[rsp+24h]  
00007FF6D60010D3  xor         ebx,ebx  
00007FF6D60010D5  mov         r11d,30D3Fh  
00007FF6D60010DB  nop         dword ptr [rax+rax]  
	{
		if (n_Sieb[i] > 0)
00007FF6D60010E0  mov         eax,dword ptr [r8]  
	{
		if (n_Sieb[i] > 0)
00007FF6D60010E3  test        eax,eax  
00007FF6D60010E5  jle         main+122h (07FF6D6001122h)  
		{
			n_Primes[n_PrimesIndex++] = n_Sieb[i];
00007FF6D60010E7  add         r10,4  
			for (int j = i + 1; j < 200000; j++)
00007FF6D60010EB  mov         rcx,r9  
00007FF6D60010EE  cmp         r9,30D40h  
00007FF6D60010F5  jge         main+122h (07FF6D6001122h)  
00007FF6D60010F7  nop         word ptr [rax+rax]  
			{
				if (n_Sieb[j] % n_Sieb[i] == 0)
00007FF6D6001100  mov         eax,dword ptr n_Sieb[rcx*4]  
00007FF6D6001107  cdq  
00007FF6D6001108  idiv        eax,dword ptr [r8]  
00007FF6D600110B  test        edx,edx  
00007FF6D600110D  jne         main+116h (07FF6D6001116h)  
				{
					n_Sieb[j] = 0;
00007FF6D600110F  mov         dword ptr n_Sieb[rcx*4],ebx  
			for (int j = i + 1; j < 200000; j++)
00007FF6D6001116  inc         rcx  
00007FF6D6001119  cmp         rcx,30D40h  
00007FF6D6001120  jl          main+100h (07FF6D6001100h)  
	}

	for (int i = 1; i < 200000; i++)
00007FF6D6001122  inc         r9  
00007FF6D6001125  add         r8,4  
00007FF6D6001129  sub         r11,1  
00007FF6D600112D  jne         main+0E0h (07FF6D60010E0h)  
				}
			}
		}
	}

[Sirius3]

@DerTürke: Du hast immer noch kein Sieb implementiert. Der Trick ist ja, dass man gar keine Divisionen braucht. Die innere Schleife sollte ungefähr so aussehen:
[codebox=c file=Unbenannt.c] n_Primes[n_PrimesIndex++] = n_Sieb;
for (int j = i + i; j < 200000; j+=i)
{
n_Sieb[j] = 0;
}
[/code]

Damit wäre bewiesen, wer absolut keine Ahnung hat .

[Sirius3]

Eine weitere simple Optimierung betrifft das Auslassen von durch 3 teilbaren Zahlen. Dadurch erhöht sich die Geschwindigkeit folglich um 1/3 .
Bei den Primzahlen bis 5000 sieht das Rennen nun so aus:
- is_prime: 107.791 Takte
- Sieb: 12.869 Takte

@DerTürke: Absolut gesehen holt damit natürlich is_prime auf .

[BlackJack]

Okay, dann nehmen wir doch mal ein ganz einfaches, wirklich „straight forward“ implementiertes Sieb in Python, ohne irgendwas speziell zu optimieren:

Code: Alles auswählen

#!/usr/bin/env python
# coding: utf8
from __future__ import absolute_import, division, print_function
from math import sqrt


def sieve(limit):
    is_prime = [True] * limit
    is_prime[0] = is_prime[1] = False
    for i in xrange(2, int(sqrt(limit)) + 1):
        if is_prime[i]:
            for j in xrange(i * i, limit, i):
                is_prime[j] = False
    return (i for i, t in enumerate(is_prime) if t)


def main():
    print(sum(sieve(2000000)))


if __name__ == '__main__':
    main()

Und lassen das gegen das folgende C-Programm antreten, welches den Algorithmus von DerTürke implementiert:
[codebox=c file=Unbenannt.c]#include <inttypes.h>
#include <stdbool.h>
#include <stdio.h>


static _Bool is_prime(uint_fast32_t candidate)
{
uint_fast32_t i;

for (i = 2; i < candidate; ++i) {
if (candidate % i == 0) return false;
if (candidate / i < i) break;
}
return candidate > 1;
}


int main(void)
{
uint_fast64_t sum = 0;
uint_fast32_t i;

for (i = 1; i < 2000000; ++i) {
if (is_prime(i)) {
sum += i;
}
}

printf("%"PRIuFAST64"\n", sum);
return 0;
}[/code]

Kompiliert mit GCC und Optimierungen wird die Testfunktion direkt in die Hauptfunktion eingefügt, so dass dort der Funktionsaufruf eingespart wird. Ich habe den generierten Quelltext mal mit sprechenderen Label-Namen versehen und ein wenig kommentiert:
[codebox=asm file=Unbenannt.asm] .intel_syntax noprefix
.section .rodata.str1.1,"aMS",@progbits,1
.format_string:
.string "%lu\n"

.section .text.startup,"ax",@progbits
.p2align 4
.globl main
main:
; candidate = 1
mov esi, 1
; sum = 0
xor edi, edi
; Align loop start to 16 byte boundary if max 10 bytes away
; or to 8 byte boundary
.p2align 4,,10
.p2align 3
.outer_loop:
; if candidate <= 2
cmp rsi, 2
jbe .L8

; if candidate is even but not <=2 it is not prime
test sil, 1
je .is_not_prime

; if candidate is three then it is prime,
; else i = 2 and run inner loop
cmp rsi, 3
mov ecx, 2
jne .inner_loop_inc
jmp .is_prime
.p2align 4,,10
.p2align 3
.inner_loop:
; div/mod candidate and i
xor edx, edx
mov rax, rsi
div rcx
; if remainder then not prime
test rdx, rdx
je .is_not_prime
; if result of division < i then candidate is prime
cmp rax, rcx
jb .is_prime
.inner_loop_inc:
; increase i and if not reached candidate then test i
add rcx, 1
cmp rcx, rsi
jne .inner_loop
.is_prime:
; sum += candidate
add rdi, rsi
.is_not_prime:
; increase candidate and test if upper limit is hit
add rsi, 1
cmp rsi, 2000000
jne .outer_loop

; print sum
mov rdx, rdi
mov esi, OFFSET FLAT:.format_string
mov edi, 1
xor eax, eax
call __printf_chk

; return 0
xor eax, eax
ret
.L8:
; if candidate is 1 then candidate = 2
; otherwise candidate stays 2
cmp rsi, 1
mov eax, 2
cmove rsi, rax
jmp .is_prime[/code]

Wie man sieht entspricht das ziemlich dem was DerTürke an Assembler gezeigt hat. Etwas schräg ist der erste Test und der Code bei `.L8`. Das hätte man als Mensch wahrscheinlich anders geschrieben. Letztlich fällt das bei der Laufzeit aber kaum ins Gewicht, denn es betrifft nur einen von zwei Millionen Schleifendurchläufen.

Interessant ist, dass der Compiler selbst gemerkt hat, das ungerade Kandidaten ausser der 2 nicht zu einer Addition auf die Summe führen und das deswegen mit einem billigen Bit-Test gesondert behandelt, um die teurere Division für jede zweite Zahl zu sparen.

Und das Ergebnis: Python ist fast doppelt so schnell.

Um mal einen Fusballtrainer zu zitieren: Schön spielen ist nichts, wenn das Resultat nicht stimmt. Man kann sich viel Mühe geben handgeklöppelten Assembler zu schreiben, und dafür sicher ein Fleissbienchen verdienen, aber wenn man dann *deutlich* von einer dynamisch interpetierten Programmiersprache überholt wird, wo sogar noch der Compilerlauf von Quelltext nach Bytecode für die VM mit in der Laufzeit enthalten ist, dann war das wohl nix. Aber dicke Töne spucken, andere hätten keine Ahnung. Süss.

[kbr]

Für alle, die sich auf die schnelle weder C noch Assembler einlesen können, habe ich ein einfaches Sieb in Python umgesetzt und dem erwähnten 'algo' gegenübergestellt. Diesen habe ich zudem etwas modifiziert, damit die gleiche Datenstruktur zurück gegeben wird. Ansonsten habe ich den 'algo' unverändert gelassen. Es fällt sofort auf, dass dieser nicht schneller sein kann. Das Sieb berechnet die ersten 2000000 Primzahlen in unter 1 sec. Bei dem 'algo' wage ich so große Zahlen gar nicht erst.

Code: Alles auswählen

def sieve(n):
    primes = [True] * n
    for i in range(2, n):
        if primes[i]:
            for j in range(i+i, n, i):
                primes[j] = False
    return primes

def algo(n):
    primes = [True] * n
    for i in range(2, n):
        for j in range(2, i):
            if not i % j:
                primes[i] = False
                break
    return primes

[BlackJack]

@kbr: Du hast `algo` falsch umgesetzt. Du müsstest die innere Schleife abbrechen sobald ``i // j < j`` gilt. Das hat DerTürke gemacht um nur bis j nur bis sqrt(i) laufen zu lassen. Wobei sich das in Assembler anbietet weil man den Modulo-Wert und den Quotienten ”kostenlos” von der Division bekommt die man sowieso durchführt.

[DerTürke]

@Sirius:
Du hasst dir wohl meinen Code nicht angeschaut, ansonsten würdest du nicht behaupten ich hätte das Sieb ( was ja auch vom Schwierigkeitsgrad nicht sehr fordernd ist ) nicht implementiert, vielmehr habe ich auch ein 2 Array of int (also wirklich keine List kein vector oder der Gleichen) benutzt um auch die im Sieb ermittelten Primzahlen zu sammeln, das mit der Ahnung war von meiner Seite nicht so gemeint sorry falls ich dir damit zu Nahe getreten bin. Ich glaube ich bin hier sowieso fehl am Platz ich sollte mich am besten in Foren aufhalten wo man keine so große Antipathie gegen Assembler pflegt. Nochmal auch an BlackJack sorry wenn ich mit meinen Äußerungen euch zu Nahe getreten bin ich wünsche euch allen gutes gelingen und eine schöne Zeit.

[kbr]

@BlackJack: Stimmt, das hatte ich übersehen und mich schon gewundert. Diese Optimierung hatte ich ursprünglich auch angewandt, was am Ergebnis aber nichts ändert.

[Sirius3]

@DerTürke: ich hab Dir sogar gesagt, was Du falsch gemacht hast. Also zu behaupten, das wäre so einfach, dann eine falsche Implementierung zu zeigen und dann noch sich zu freuen, dass der »eigene Algo« schneller sei, ...
Aber hier zu behaupten, wir würden ein gutes Assemblerprogramm nicht wertschätzen, das ist eine unverschämte Unterstellung .

[DerTürke]

@Sirius

@DerTürke: Du hast immer noch kein Sieb implementiert. Der Trick ist ja, dass man gar keine Divisionen braucht. Die innere Schleife sollte ungefähr so aussehen:

Code: Alles auswählen

1.         n_Primes[n_PrimesIndex++] = n_Sieb[i];
2.         for (int j = i + i; j < 200000; j+=i)
3.         {
4.            n_Sieb[j] = 0;
5.         }

Damit wäre bewiesen, wer absolut keine Ahnung hat .

Ich wollte eigentlich dieses Thema mal so langsam abhaken, daß muss ich aber dann doch noch mal fragen Siri
Frage : "Wo in der Implementierung habe ich denn eine divison ?" , und dein Vorschlag wie ich es machen müsste
Zitat: "Der Trick ist ja, dass man gar keine Divisionen braucht. Die innere Schleife sollte ungefähr so aussehen:"
würde doch das ganze Array mit 0 füllen oder ? sehe ich da was nicht ?

Trotzdem vielen Dank für die Mühe

[DerTürke]

@Sirius:

Nein ok deine Schleife habe ich im nachhinein nochmal analysiert, dadurch das der Inkrement von j mit dem vielfachen von i durchgeführt wird, werden natürlich nur die vielfachen des zuvor emrittelten Primes auf 0 gesetzt.
sorry hasst komplett recht und dadurch sparst du dir den Modulo
ok habe kapiert, tatsächlich eine sehr gut Implementierung

Danke und Gruß
DT

[BlackJack]

@DerTürke: Wer hat denn hier Antipathien gegen Assembler? Ich mag Assembler! Habe hier im Thema doch schon eine Menge davon gepostet. Und in meinem Leben auch schon so einiges geschrieben und gelesen. Auch Sirius3 hat sich nicht gescheut auf diese hardwarenahe Ebene zu gehen und mal zu schauen und zu zeigen was ein optimierender Compiler für Code erzeugt und das der an der Stelle per Hand kaum zu verbessern ist.

Unser Problem mit Dir ist nicht, das Du eine Lösung in Assembler geschrieben hast, sondern das die unglaublich ineffizient ist. So ineffizient das eine in keiner Weise optimierte Python-Lösung *deutlich* schneller ist. Trotzdem verteidigst Du Deinen Assembler-Code mit der Argumentation das nichts an handgeschriebenen Assembler heran kommt. Selbst wenn das in diesem Fall stimmen *würde*, sind solche Mikrooptimierungen völlig unsinnig solange Du diesen grottenschlechten Algorithmus verwendest. Selbst bei dem hat der GCC geringfügig besseren Code produziert, in dem er in der innersten Schleife ein ``TEST reg, reg`` statt eines ``CMP reg, 0`` verwendet hat, wie ich weiter oben gezeigt habe.

So sieht das Sieb in C komplett & korrekt aus:
[codebox=c file=Unbenannt.c]#include <inttypes.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>

#define LIMIT 2000000

static _Bool is_prime[LIMIT];


int main(void)
{
uint_fast64_t sum = 0;
uint_fast32_t i, j;

memset(is_prime, true, sizeof(is_prime));
is_prime[0] = is_prime[1] = false;
for (i = 2; i < ((uint_fast32_t) sqrt(LIMIT)) + 1; ++i) {
if (is_prime) {
for (j = i * i; j < LIMIT; j += i) {
is_prime[j] = false;
}
}
}

for (i = 2; i < LIMIT; ++i) {
if (is_prime) {
sum += i;
}
}

printf("%"PRIuFAST64"\n", sum);
return 0;
}[/code]
Das ist 57× schneller als Dein Algorithmus!

Und der C-Compiler spuckt auch hier wieder ziemlich guten Assembler aus (von mir kommentiert):
[codebox=asm file=Unbenannt.asm] .intel_syntax noprefix
.section .rodata.str1.1,"aMS",@progbits,1
.format_string:
.string "%lu\n"

.section .text.startup,"ax",@progbits
.p2align 4
.globl main
main:
; set all entries in is_prime to true
mov edx, 2000000
mov esi, 1
mov edi, OFFSET FLAT:is_prime
call memset

; i = 2
mov edx, 2
.p2align 4,,10
.p2align 3
.outer_loop:
; test if is_prime is prime
cmp BYTE PTR is_prime[rdx], 0
je .not_prime
; j = i * i
mov rax, rdx
imul rax, rdx
.p2align 4,,10
.p2align 3
.inner_loop:
; is_prime[j] = false
mov BYTE PTR is_prime[rax], 0
; j += i
add rax, rdx
; is j < limit?
cmp rax, 1999999
jbe .inner_loop
.not_prime:
; i += 1
add rdx, 1
; is i == sqrt(LIMIT)?
cmp rdx, 1415
jne .outer_loop

; i = 2
mov eax, 2
; sum = 0
xor dx, dx
.p2align 4,,10
.p2align 3
.sum_loop:
; is is_prime == false?
cmp BYTE PTR is_prime[rax], 0
; tmp = sum + i
lea rcx, [rax+rdx]
; conditionally move tmp to sum if i is prime
cmovne rdx, rcx
; i += 1
add rax, 1
; if i <= LIMIT goto loop start
cmp rax, 2000000
jne .sum_loop

; print the sum
mov esi, OFFSET FLAT:.format_string
mov edi, 1
xor eax, eax
call __printf_chk
xor eax, eax
ret[/code]

Die Addition mit ``LEA`` und den bedingte ``MOV``-Befehl finde ich clever. Auf so etwas muss man als Mensch erst einmal kommen. Also entweder selbst, oder weil man diese Vorgehensweise schon mal gesehen und sich gemerkt hat.

Das die Register nicht immer als 64-Bit register behandelt werden, wenn sie nicht müssen, ist auch etwas was ein Compiler wunderbar hinbekommt, wo man als Mensch gerne mal Fehler macht wenn man Code entwickelt oder verändert. Beispielsweise das die 64-Bit-Summe die im rdx-Register aufsummiert wird, mit einem ``xor dx,dx`` auf 0 gesetzt wird, ist fehleranfällig bei Veränderungen von Quelltext davor, weil man an der Stelle ja sicher sein muss, das der Inhalt in dem Register zu dem Zeitpunkt kleiner als 65536 war. Das mag in diesem Fall recht offensichtlich sein, aber ein Compiler macht das halt immer richtig, egal von wievielen Punkten aus so ein Code angesprungen werden kann und wie weit der letzte Zugriff auf das Register auch entfernt sein mag.

Geschwindigkeitsvergleich ”grafisch”:
[codebox=text file=Unbenannt.txt]Laufzeit (< besser):

DerTürke-Algo *********************************************************
Python-Sieb ******************************
C-Sieb *[/code]

[DerTürke]

Hallo Jungs,

habe mich erkundigt und muss leider bestätigen das mein alter Alg.o grotten schlecht war.
Der neue ist von Miller Rabin super schnell ermittelt die Primzahlen von 1-1.000.000.000 in weinger als eine sek.

hier isser:
[codebox=cpp file=Unbenannt.cpp]#include<iostream>
#include<stdlib.h>

using namespace std;

//bool IsPrime(int pValue);
bool IsPrimeExt(uint64_t pZahl, int pBasis);

int main()
{
int m_Value;
for (int i = 10; i < 1000000001; i++)
{
bool bIsPrime = IsPrimeExt(i, 2) & IsPrimeExt(i, 7) & IsPrimeExt(i, 61);
/*
if (i % 100000000 == 0)
{
cout << "Zähler gerade bei " << i << endl;
}
*/
}
cout << "Fertig" << endl;
cin >> m_Value;
return 0;
}

bool IsPrimeExt(uint64_t pZahl, int pBasis)
{
bool retValue = false;
uint64_t m_ZahlMinus1 = pZahl - 1;
uint64_t m_Ungerade = m_ZahlMinus1;
int m_2T = 0;
while ((m_Ungerade & 1) == 0)
{
m_Ungerade >>= 1;
m_2T += 1;
}

uint64_t m_Power = pBasis;
uint64_t m_T = pBasis;
while (m_Ungerade >>= 1)
{
m_Power = m_Power*m_Power % pZahl;
if (m_Ungerade & 1)
{
m_T = m_T*m_Power % pZahl;
}
}

if ((m_T == m_ZahlMinus1) || (m_T == 1))
{
retValue = true;
}

for (int k = 1; k < m_2T; k++)
{
m_T = (m_T*m_T) % pZahl;
if(m_T == m_ZahlMinus1)
{
retValue = true;
break;
}
if (m_T <= 1)
{
retValue = false;
break;
}
}
return retValue;
}[/code]

Vielleicht kann man diesen auch in Python implementieren
Gruß
DerTürke

[BlackJack]

Kann man, aber der ist ein probalistischer Test, das heisst der kann nur sicher sagen was *keine* Primzahl ist, aber nicht was ganz sicher eine Primzahl ist.

[DerTürke]

Ja habe ich gerade erfahren,

bei den Primzahlen bis 100 fehlt die 13, 37 usw...
Schade war total begeistert von diesem nun ja weiter schauen.

[DerTürke]

Ich hab's nochmal getestet aber nur zu Basis 2 da klappt es seltsamer Weise sollte man lt. Rabin Miller aber mit 2 Basen probieren und zwar mit 2 und 3, Naja werde dann wohl doch bei meinem Grotten schlechten bleiben es sei denn ihr kennt einen Algo. ( bitte nicht das Sieb ) der effiezienter ist.

VG
DT

[BlackJack]

@DerTürke: Was hast Du gegen das Sieb? Wobei es ja nicht nur *das* Sieb gibt. Sieve of Atkins wäre auch eine Möglichkeit.

Oder verbessere Deinen Algorithmus. Der hat den Vorteil parallelisierbar zu sein. Damit bekomme ich den auf meinem alten Laptop an die Python Implementierung heran. Zwei Kerne + Hyperthreading. Mit mehr und oder echten Prozessoren/Kernen, kann der dann auch gegen Python gewinnen. Wenn man das nicht in Assembler macht, geht das sogar sehr einfach. Bei dem C-Programm einfach ein OpenMP-Pragma um die äussere Schleife zu parallelisieren:
[codebox=c file=Unbenannt.c]int main(void)
{
uint_fast64_t sum = 0;
uint_fast32_t i;

#pragma omp parallel for reduction(+:sum)
for (i = 1; i < 2000000; ++i) {
if (is_prime(i)) {
sum += i;
}
}

printf("%"PRIuFAST64"\n", sum);
return 0;
}[/code]

[DerTürke]

Ja danke
gute Idee mit OpenMP.
Im Moment studiere ich gerade den Algorithmus von Magenheimer "Integer multiplication and Division on the HP Architecture" ich glaube wenn ich den verstehe dann kann ich die ineffiziente Division ( bis zu 200 Taktzyklen ) durch eine Multiplikation ( irgendwie mit Umkehrwert ) ersetzen, und dass ganze verpacke ich dann in VPMULQD (SIMD) 2 x 64 Bit (pseudo divs) naja ich habs auf jeden Fall vor, mal schaun wann ich die Lust dran verliere

VG
DT