A propos de nombres premiers

Partagez votre expérience de PureBasic avec les autres utilisateurs.
Avatar de l’utilisateur
djes
Messages : 4252
Inscription : ven. 11/févr./2005 17:34
Localisation : Arras, France

Re: A propos de nombres premiers

Message par djes »

Un petit HS sur un sujet qui vous intéressera sans doute : http://centenaire-shannon.cnrs.fr/
Ollivier
Messages : 4190
Inscription : ven. 29/juin/2007 17:50
Localisation : Encore ?
Contact :

Re: A propos de nombres premiers

Message par Ollivier »

Bonjour à vous,

je pense avoir trouvé la manière de coder pour expliquer la virtualisation de façon à ce qu'un débutant puisse comprendre.

C'est évidemment revenir en arrière sur les avancées du sujet, mais j'avais bien dit que je tâcherai de faire un code "tuto" prochainement.

Ça n'aidera pas fweil malheureusement, mais ça permettra peut-être à d'autres personnes de trouver un intérêt ludique de la recherche des nombres premiers.

Pour le reste que je lis au fur et à mesure, j'apprécie grandement. La technique de Miller-Rabin proposée par fweil est une nouveauté pour moi : il me faudra beaucoup de temps pour la comprendre.
Ollivier
Messages : 4190
Inscription : ven. 29/juin/2007 17:50
Localisation : Encore ?
Contact :

Re: A propos de nombres premiers

Message par Ollivier »

@FWeil

L'explication de SPH en début de sujet ne se rapproche-t-elle pas de cette technique de Miller-Rabin?
Avatar de l’utilisateur
SPH
Messages : 4721
Inscription : mer. 09/nov./2005 9:53

Re: A propos de nombres premiers

Message par SPH »

djes a écrit :Un petit HS sur un sujet qui vous intéressera sans doute : http://centenaire-shannon.cnrs.fr/
Je n'ai pas compris et vue ce qu'il y a d'interessant dans ce site. A part 4 pages a lire sur la vie de shanon, je n'ai rien vu (je ne sais pas ou cliquer)...

Merci pour tes lumieres 8)
http://HexaScrabble.com/
!i!i!i!i!i!i!i!i!i!
!i!i!i!i!i!i!
!i!i!i!
//// Informations ////
Intel Core i7 4770 64 bits - GTX 650 Ti
Version de PB : 6.00 - 64 bits
Avatar de l’utilisateur
djes
Messages : 4252
Inscription : ven. 11/févr./2005 17:34
Localisation : Arras, France

Re: A propos de nombres premiers

Message par djes »

SPH a écrit :
djes a écrit :Un petit HS sur un sujet qui vous intéressera sans doute : http://centenaire-shannon.cnrs.fr/
Je n'ai pas compris et vue ce qu'il y a d'interessant dans ce site. A part 4 pages a lire sur la vie de shanon, je n'ai rien vu (je ne sais pas ou cliquer)...

Merci pour tes lumieres 8)
Il faut utiliser les flèches du clavier, la molette de la souris, cliquer sur les liens...
Ollivier
Messages : 4190
Inscription : ven. 29/juin/2007 17:50
Localisation : Encore ?
Contact :

Re: A propos de nombres premiers

Message par Ollivier »

@DJes

Quand j'ai vu que j'allais avoir du mal, je suis allé voir la biographie sur Wikipedia.
PAPIPP
Messages : 534
Inscription : sam. 23/févr./2008 17:58

Re: A propos de nombres premiers

Message par PAPIPP »

Bonjour à tous

Voici une optimisation du prg précédent dans la recherche par zone des nombres premiers.
Travailler aux limites de la machine pose quelques difficultés en ASM.
Je me suis imposé une contrainte supplémentaire c’est de faire fonctionner le prg en 32 ou 64 bits sans avoir à recompiler ou en utilisant les mêmes instructions.

La toute première difficulté est de pouvoir réaliser certaines opérations addition division etc..
En effet si vous rechercher les diviseurs possibles de tout nombre proche de pow(2,63)-1 c'est-à-dire
$7FFFFFFFFFFFFFFF =9223372036854774807 il n’est pas nécessaire de travailler avec tout diviseur supérieur à sqr($7FFFFFFFFFFFFFFF)= 3037000499,97
valeur qui dépasse théoriquement la possibilité de la capacité d’un registre 32 bits=2147483647 en valeur signée.
Mais comme certaines instructions ASM permettent de travailler en valeur non signée ce qui permet d’atteindre pow(2,32)=
4294967296 Pb résolu en faisant attention aux instructions utilisées.

La deuxième difficulté ce trouve dans la division d’un nombre très proche de pow(2,63)-1 par un diviseur très petit on peut avoir une erreur de dépassement de capacité avec la division classique. ATTENTION le PRG suivant plante pour dépassement de la capacité.

Code : Tout sélectionner

 Macro _q_t_
"
EndMacro
Macro _n (__n)
_q_t_#__n#=_q_t_+Str(__n)+" "
EndMacro
Macro _Q (__Q)
_q_t_#__Q#=_q_t_+Str(__Q)+" "
EndMacro
nb.q=$123456789ABCDEF0
debug _q(nb)
; DIVIS.l=$FEDC
DIVIS.l=$FF 
; ; ; l'opération si dessous plante  dépassement de capacité sur type entier (dividende trop grand pour diviseur trop petit)
enableasm
mov eax, dword [v_nb] 
mov edx, dword [v_nb+4] 
mov ecx, dword [v_DIVIS] 
div ecx ;;; eici ça plante dividente trop grand pour diviseur trop petit

Il faut décomposer l’opération en deux parties comme cela.

Code : Tout sélectionner

Macro _q_t_
 		"
 EndMacro
 Macro _n (__n)
 _q_t_#__n#=_q_t_+Str(__n)+" "
 EndMacro
 Macro _Q (__Q)
 _q_t_#__Q#=_q_t_+Str(__Q)+" "
 EndMacro
nb.q=$123456789ABCDEF0
debug _q(nb)
; DIVIS.l=$FEDC
DIVIS.l=$F
debug "**************** début des 2 sortes de divisions ****************************"
define quotientp.q,DIVMAX.q,restd.l
pas=2
FOR IJ=0 TO 10
  debug "****"
  np.q=nb
  debug _q(nb)+_q(nb/DIVIS)+""+_q(nb%DIVIS)
  EnableASM
  MOV ecx,dword[v_DIVIS]
  XOR EDX,EDX
  MOV eax,dword[v_nb+4]
  DIV ecx
  MOV dword[v_quotientp],eax
  MOV eax,dword[v_nb]
  DIV ecx
  MOV dword[v_restd],edx
  MOV edx,dword[v_quotientp]
  MOV dword[v_nb+4],edx
  MOV dword[v_nb],eax  ;************************************************ Recherche de la racine carré ***************************************************************
  DisableASM
  debug _q(np)+"___"+_n(nb)+"__"+_n(restd)
 next 
 debug "************* FIN ***********"

Le troisième problème est de considérer que les nombres très grands sont dans une zone ou la racine carrée reste constante. Or le prg doit travailler dans toute zone de 1 à $7FFFFFFFFFFFFFFF avec une variation de la racine carré de 1 à 3037000499,97 ce qui fait un nombre considérable de calculs de la racine carrée. Pour éviter d’avoir à calculer la racine carrée à l’incrémentation du dividende il faut prévoir la valeur du dividende qui incrémentera la racine carrée de 1
En effet on peut prévoir à partir d’un nombre X et de sa racine RAC quant la racine carrée va changer car (A+1)^2=A^2+1+2A
si A=SQR(x)=RAC alors (SQR(X)+1)^2= X+1+2*A= X+1+2*RAC
On peut encore plus optimiser sachant que le diviseur max doit appartenir à la colonne susceptible de contenir des nb premiers
(cette optimisation n’est pas encore faite c’est pour la prochaine optimisation)


La quatrième difficulté c’est le temps de préparation de la table des différences des p_fact(x)
; exemple des p_fact(x)
; p_fact(x)--modulo--------------------Nb_occurrences des différences------------------------------Taux----------------Gain
; 3__________6__________________________2________________________________0,3333_________66,66%
; 5__________30_________________________8________________________________0,2666_________73,33%
; 7__________210________________________48_______________________________0,2285_________77,14%
; 11_________2310_______________________480______________________________0,2078_________79,22%
; 13_________30030______________________5760_____________________________0,1918_________80,81%
; 17_________510510 ____________________92160____________________________0,18052________81,94%
; 19_________9699690___________________1658880___________________________0,17102________82,89%
; 23_________223092870_________________36495360__________________________0,16358________83,64%

Or plus p_fact(x) est élevé plus le temps de préparation de la table est important donc le temps de préparation devient très important pour un p_fact(23) on avait plus de 50 secondes de préparation.

Pour diminuer ce temps il faut mémoriser la table des différences sous forme binaire au moment de la conception.
Il suffira de charger cette table en quelques milli-secondes.

Pour utiliser le prg de recherche des Nb premiers dans une zone, il faut commencer par constituer le fichier binaire de la table des différences avec le prg suivant.
Ce prg n’est à utiliser qu’une seule fois. Par contre si vous désirer voir les différences de temps, créez autant de fichiers binaires que vous désirez
avec p_fact(x) x [2 3 5 7 11 13 17 19 23] au delà de 23 les tables générées sont trop importantes (voir tableau ci-dessus).

Code : Tout sélectionner

; exemple des p_fact(x)
; p_fact(x)--modulo-----------------Nb_occurrences des différences--------Taux-----------Gain
; 3__________6__________________________2_________________________________0,3333_________66,66%
; 5__________30_________________________8_________________________________0,2666_________73,33%
; 7__________210________________________48________________________________0,2285_________77,14%
; 11_________2310_______________________480_______________________________0,2078_________79,22%
; 13_________30030______________________5760______________________________0,1918_________80,81%
; 17_________510510 ____________________92160_____________________________0,18052________81,94%
; 19_________9699690____________________1658880___________________________0,17102________82,89%
; 23_________223092870__________________36495360__________________________0,16358________83,64%

 Macro _q_t_
 		"
 EndMacro
 Macro _n (__n)
 _q_t_#__n#=_q_t_+Str(__n)+" "
 EndMacro
 Macro _Q (__Q)
 _q_t_#__Q#=_q_t_+Str(__Q)+" "
 EndMacro
EnableExplicit

OpenConsole("Résultats partiel")
Dim T_MODULO.l(12)
Structure colonne
;   nb.q
	prem.a
; 	dif_prec.l
	dif_prec.a

; 	Dif_act.l
EndStructure
Structure DIVIS
  NBPFACT.a
  PREMDIV.l
  NBSEQ.l
  MODULO.q
  ;   Array TDIF.a(40000000)  ;; *** ATTENTION   A utiliser avec DATA.A et read.A  mais très lent avec p_fact(23) P_fact(19) est préférable remplacer 40 par 2
  Array TDIF.a(40000000)  
EndStructure
define SEQD.DIVIS,mess$,nb$,nbs,rest.l,i,modulo,elem,NBprem,Tdep.q,inb,inbp,prem_p,PREM_NB_PREM_AV,NB_PREM,SOM_DIF,vecteur$,NBSEQ,lvecteur,rapportq.d,RAPPORT.D
define ADRDSEQ.q,AdresdebF.q,AdresFinF.q,Adresdeb.q,delta

; **************** Réalisation du vecteur modulo npp à partir d'une table des nombres premiers des 100 premiers nombres *********
;***** Choix du premier Nb premier pour lequel le vecteur modulo sera appliqué.
;***** pour éviter la génération d'un vesteur trop important nous limiterons ce choix aux 23 premiers nombres premiers soit :2 3 5 7 11 13 17 19 23 29 31
;***** Recherche du modulo < 2*3*5*7*11=2310 éléments du vecteur
SAISIE0:
mess$+"Filtre de réduction des diviseurs"
If Len(mess$)>120
	MessageRequester("Erreur","Colonne Div max 2 3 5 7 11 13 17 19 23 29 31"+#CR$+"Relancez le prg") ;
	End
EndIf
nb$=InputRequester(mess$,"Colonne Div max 2 3 5 7 11 13 17 19 23 29 31","5") ;
If Len(nb$)>Len("31")
EndIf
nbs=Val(nb$)
If nbs<1 Or nbs>31
	Goto SAISIE0
EndIf
rest.l=nbs%30
If rest=2 Or rest=3 Or rest=5 Or rest=7 Or rest=11 Or rest=13 Or rest=17 Or rest=19 Or rest=23 Or rest=29
Else
  Goto SAISIE0
EndIf 
; NBPFACT=NBS
; FINSAISIE0:
Restore lab_pnbp
i=0
Modulo=1
Repeat
	Read.l ELEM
	If ELEM<>0
		T_modulo(I)=elem
		Modulo*elem
		i+1
	EndIf
Until ELEM=0 Or ELEM=>nbs
redim SEQD\TDIF(Modulo+10)
SEQD\MODULO=MODULO 
SEQD\NBPFACT=nbs
NBprem=i-1
nbs=ELEM
Dim tcol.colonne(SEQD\MODULO+2)
;*** Recherche des colonnes ayant des nb premiers ***
tcol(0)\prem=2
tcol(1)\prem=2
Tdep.q=ElapsedMilliseconds()
For inb=2 To modulo+1
	For inbp=0 To nbprem
		If inb%t_modulo(inbp)=0
			tcol(inb)\prem=2
		EndIf
	Next
Next
prem_p=0
PREM_NB_PREM_AV=0
NB_PREM=0
SOM_DIF=0
vecteur$="DIV_"+Str(nbs)+":"+#CRLF$
;****** Recherche des différence entre NB premiers *****
If CreateFile(0,"DIVA_"+Str(nbs)+".PB")         ; création d'un nouveau fichier texte...
  NBSEQ=0
  For inb=2 To modulo+1
;     tcol(inb)\nb=inb
		If tcol(inb)\prem=0 And inb>nbs
			If prem_p=0
				prem_p=inb
				PREM_NB_PREM_AV=inb
				vecteur$+"DATA.A  "+Str(inb)+","
				SEQD\PREMDIV=inb
				NBSEQ+1
			Else
				NB_PREM+1
				If nb_prem%101=0
				  If nb_prem%201=0
				    PrintN(_n(NB_prem)+_n(prem_p)+_n(tcol(inb)\dif_prec)+_n(som_dif)+_n(ElapsedMilliseconds()-Tdep)+_n(lvecteur))
					EndIf
					If nb_prem%5=0  ; ici tous les =101*x  c'est avec 5 que l'on obtient le meilleur résultat
					  vecteur$+Str(inb-prem_p)+#CRLF$
; 					  SEQD\TDIF(NBSEQ-1)=inb-prem_p
					  SEQD\TDIF(NBSEQ)=inb-prem_p

					  NBSEQ+1
						lvecteur=Len(vecteur$)
						WriteString(0,vecteur$) ;
						vecteur$="Data.A "
					Else
					  vecteur$+Str(inb-prem_p)+#CRLF$+"Data.A "
; 					  SEQD\TDIF(NBSEQ-1)=inb-prem_p
					  SEQD\TDIF(NBSEQ)=inb-prem_p

					  NBSEQ+1
					EndIf
					
					tcol(inb)\dif_prec=inb-prem_p
					SOM_DIF+(inb-prem_p)
; 					tcol(prem_p)\Dif_act=inb-prem_p
					prem_p=inb
				Else
					tcol(inb)\dif_prec=inb-prem_p
					vecteur$+Str(inb-prem_p)+","
; 					SEQD\TDIF(NBSEQ-1)=inb-prem_p
					SEQD\TDIF(NBSEQ)=inb-prem_p
					NBSEQ+1
					SOM_DIF+(inb-prem_p)
; 					tcol(prem_p)\Dif_act=inb-prem_p
					prem_p=inb
				EndIf
			EndIf
		EndIf
	Next
	;  tcol(modulo+2)\dif_prec=nb_prem+modulo-(modulo+1)
	NB_PREM+1
	tcol(modulo+2)\dif_prec=PREM_NB_PREM_AV-1
	SEQD\TDIF(NBSEQ)=PREM_NB_PREM_AV-1
	SEQD\TDIF(0)=PREM_NB_PREM_AV-1
	SEQD\NBSEQ=NBSEQ
	redim SEQD\TDIF(NBSEQ+2)

; 	tcol(1)\Dif_act=PREM_NB_PREM_AV-1
; 	tcol(1)\nb=1
	SOM_DIF+(PREM_NB_PREM_AV-1)
	rapportq.d=100.0*(modulo-Nb_prem)/modulo
	vecteur$+Str(PREM_NB_PREM_AV-1)+",0 ;"+#CRLF$+";; Modulo="+Str(modulo)+" Nb Elements dans vecteur="+Str(NB_prem)+"  NBPFACT="+str(SEQD\NBPFACT)+" Nb sequences="+str(NBSEQ)+"  GAin="+Str(rapportq)+"%"
	
	WriteString(0,vecteur$) ;
	CloseFile(0)            ; ferme le fichier précédemment ouvert et enregistre les données
	PrintN(_n(ElapsedMilliseconds()-Tdep))
	RAPPORT.D=modulo/NB_PREM
	MessageRequester("C'est tout bon","Somme des nb du vecteur="+_n(SOM_DIF)+#CRLF$+_n(modulo)+#CRLF$+"voir le fichier pour le détail du vecteur"+#CRLF$+"Fichier: "+GetCurrentDirectory()+"DIVA_"+Str(nbs)+".PB"+#CRLF$+"Nomnre d'éléments dans le vecteur="+Str(NB_PREM)+" Rapport="+StrD(RAPPORT))
Else
	MessageRequester("Information","Impossible de créer le fichier! DIVA_"+Str(nbs)+".PB")
EndIf
; Adresdeb.l=@SEQD\NBPFACTStructure DIVIS
;   NBPFACT.a
;   PREMDIV.l
;   NBSEQ.l
;   MODULO.q
;   ;   Array TDIF.a(40000000)  ;; *** ATTENTION   A utiliser avec DATA.A et read.A  mais très lent avec p_fact(23) P_fact(19) est préférable remplacer 40 par 2
;   Array TDIF.a(40000000)  
ADRDSEQ.q=@SEQD
AdresdebF.q=@SEQD\NBPFACT
AdresFinF.q=@SEQD\MODULO+sizeof(quad)
Adresdeb.q=@SEQD\TDIF(0)
delta=AdresFinF-AdresdebF
If CreateFile(1, "MEMA_"+Str(nbs)+".bin")
  WriteData(1,ADRDSEQ ,delta )
;   WriteData(1,Adresdebf , 17)
  WriteData(1,Adresdeb , SEQD\NBSEQ+1)
  closefile(1)
Else
  messagerequester("ATTENTION", "Fichier "+ "MEMA_"+Str(nbs)+".bin non créer")
endif

Input()
CloseConsole()

DataSection
	lab_pnbp:
	Data.l 2,3,5,7,11,13,17,19,23,29,31,0
EndDataSection

DataSection
  lab_pnbp2:
  Data.L   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
  Data.L  101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199
  Data.L  211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293
  Data.L  307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397
  Data.L  401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499
  Data.L  503,509,521,523,541,547,557,563,569,571,577,587,593,599
  Data.L  601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691
  Data.L  701,709,719,727,733,739,743,751,757,761,769,773,787,797
  Data.L  809,811,821,823,827,829,839,853,857,859,863,877,881,883,887
  Data.L  907,911,919,929,937,941,947,953,967,971,977,983,991,997,0
EndDataSection
; ***************** Ci dessous les Nombres premiers pour les 1000 premiers nombres
; 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 547 557 563 569 571 577 587 593 599
; 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691
; 701 709 719 727 733 739 743 751 757 761 769 773 787 797
; 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887
; 907 911 919 929 937 941 947 953 967 971 977 983 991 997

Voici enfin le programme de recherche par zone des Nb premiers.

Code : Tout sélectionner

; Le nombre 1 n'est pas par convention un nombre premier.

; exemple des p_fact(x)
; p_fact(x)--modulo-----------------Nb_occurrences des différences--------Taux-----------Gain
; 3__________6__________________________2_________________________________0,3333_________66,66%
; 5__________30_________________________8_________________________________0,2666_________73,33%
; 7__________210________________________48________________________________0,2285_________77,14%
; 11_________2310_______________________480_______________________________0,2078_________79,22%
; 13_________30030______________________5760______________________________0,1918_________80,81%
; 17_________510510 ____________________92160_____________________________0,18052________81,94%
; 19_________9699690____________________1658880___________________________0,17102________82,89%
; 23_________223092870__________________36495360__________________________0,16358________83,64%
Macro _q_t_
				"
EndMacro
Macro _n (__n)
_q_t_#__n#=_q_t_+Str(__n)+" "
EndMacro
Macro _Q (__Q)
_q_t_#__Q#=_q_t_+Str(__Q)+" "
EndMacro

EnableExplicit
; DisableDebugger
; 2 3 5 7 11 13 17 19 23
; dim Nb_prem(10)

		Structure DIVIS
	NBPFACT.a
	PREMDIV.l
	NBSEQ.l
	MODULO.q
		Array TDIF.a(37000000)  ;; *** ATTENTION   A utiliser avec DATA.A et read.A  mais très lent avec p_fact(23) P_fact(19) est préférable remplacer 40 par 2
;   Array TDIF.a(2000000)  
EndStructure

Global result$,rest,fichier$
Define nb$,nb.q,DIVIS.q,pas,mess$,quotient.q,ind.l,SEQD.DIVIS,MIN.Q=1E12,pos, logmin, logmax, rmes,MIN_DEP.Q,ZONEP.Q,presultat$,resultat$
Define MAX.Q=MIN+1E4,nbg$,nbd$,pose.l,t_ind.l,t_col.l,t_ind_dep.l,MAXMAXDIV.q,RAC2.Q,MINRAC2.Q,nbseq.l,nbseqp,deb0.q
Define ind_dep,DIFF,FLAGFIN,IND_NB,CENTM,DIVISM,INDDD, DIVIS$,MAXDIV.q,dep_time.q,cent,nbseq_1.q,deb00.q,ADRDSEQ.q,AdresdebF.q,AdresFinF.q,Adresdeb.q,delta

Procedure.s choixdiv()
	Protected     COUR_DIR$,Filtre$,fichier$,Tfichier,fichierp$,leftfich$,RIGHTfich$
	COUR_DIR$ = GetCurrentDirectory()
	Filtre$ = "MEMA (MEMA_*.BIN)|MEMA_*.BIN;|Tous les fichiers (*.*)|*.*"
	fichier$=OpenFileRequester("Choisissez un fichier MEMA ou annulez", COUR_DIR$+"\MEMA_17.BIN", Filtre$, 0)
	fichierp$=UCase(Trim(GetFilePart(fichier$)))
	Tfichier=FileSize(Fichier$)
	leftfich$=Left(fichierp$,5)
	RIGHTfich$=Right(fichierp$,4)
	If Tfichier<1 Or Left(fichierp$,5)<>"MEMA_" Or Right(fichierp$,4) <>".BIN"
		MessageRequester( "ATTENTION", "fichier vide ou le nom n'est pas conforme "+ _n(Tfichier)+"  "+fichierp$)
		End
	EndIf
	ProcedureReturn fichier$
EndProcedure
FICHIER$=choixdiv()
deb00.q=ElapsedMilliseconds()
;************************************************************************************************************************************
ADRDSEQ=@SEQD
AdresdebF=@SEQD\NBPFACT
AdresFinF=@SEQD\MODULO+SizeOf(quad)
Adresdeb=@SEQD\TDIF(0)
delta=AdresFinF-AdresdebF

If OpenFile(2,FICHIER$);,#PB_File_SharedRead|#PB_File_NoBuffering)
	ReadData(2,ADRDSEQ ,delta)
	ReDim SEQD\TDIF(SEQD\NBSEQ+1)
	ReadData(2, Adresdeb,SEQD\NBSEQ+1)
	CloseFile(2)
EndIf
deb0.q=ElapsedMilliseconds()


mess$=""
SAISIE2:
mess$+"MIN Blanc +Zone à explorer avec MAX=MIN+1Ex  <= 1E2 +1E3 "

If Len(mess$)>250
		MessageRequester("Erreur","Donnez Min BLANC et Zone à explorer sous forme 1E3 +1E3 1E4 +1E4  1E9 +1E4 "+#CR$+"Relancez le prg") ;
		End
EndIf
nb$=UCase(LTrim(InputRequester(mess$,"MIN Blanc +Zone à explorer avec MAX=MIN+1Ex <= 1E2 +1E3 ","1E9 +1E4")));
pos=FindString(nb$," ")
nbg$=Trim(Left(nb$,pos))
nbd$=Trim(Right(NB$,Len(nb$)-pos))
POSE=FindString(nbg$,"E")
If pose>0
	MIN=ValD(nbg$)
Else 
	MIN=Val(nbg$)
EndIf 
POSE=FindString(nbd$,"E")
If pose>0
	ZONEP=ValD(nbd$)
Else 
	ZONEP=Val(nbd$)
EndIf 
MIN_DEP=MIN

MAX=MIN+ZONEP

If MAX<=MIN Or MIN<0 Or MIN=>$7FFFFFFFFFFFFFFF ;;; = 9223372036854775807 valeur max positive en .q
		Goto SAISIE2
EndIf 
logmin=Log10(MIN)
logmax=Log10(max-min)
If logmin+logmax>22
		rmes=MessageRequester(" ATTENTION Temps très long","Oui=>continuez Non=>donnez autre zone",#PB_MessageRequester_YesNo)
		If Rmes = 6        ; le bouton Oui a été choisi (Resultat = 6)
		Else               ; le bouton Non a été choisi (Resultat = 7)
				Goto SAISIE2
		EndIf
EndIf
MIN_DEP=MIN
resultat$=""
presultat$=" 2  3  5  7 11 13 17 19 23 29 31"
IF MIN < SEQD\NBPFACT 
	pose=findstring(presultat$,str(SEQD\NBPFACT))+1
	resultat$=left(presultat$,pose)
	MIN=SEQD\NBPFACT
endif   

; ******** sysnchronistion des Nombres à tester avec le vecteur des différences ***********
Define restd.q=MIN%SEQD\MODULO 
t_ind=0
t_col=1

For ind=0 To SEQD\MODULO 
	If restd<=t_col
		t_ind_dep=t_ind
		Break
	Else 
	t_col+seqd\TDIF(t_ind)
	t_ind+1
	EndIf
Next
DIFF=t_col-restd

MIN+DIFF  
;************** Fin de sysnchronisation *************************************************
OpenConsole()

divis.q=SEQD\PREMDIV
flagfin=0
; maxdiv.q=Pow(min,0.5)
maxdiv.q=Sqr(min)
MAXMAXDIV.q=Sqr(MAX)
;;;*******************************************************************************************************************************************************************
;;; On peut anticiper la prochaine racine différente de celle actuelle (A+1)^2 =A^2+1+2A si A=SQR(X)=RAC donc (RAC^2)=X alors => (RAC+1)^2= (RAC^2)+1+2*RAC=X+1+(2*RAC)
;;;*******************************************************************************************************************************************************************

RAC2.Q=(2*maxdiv)+1
MINRAC2.Q=MIN+RAC2
If MINRAC2<0
MINRAC2=9223372036854775807
EndIf  
PrintN(_n(MAXMAXDIV-MAXDIV)+_n(RAC2)+_N(MINRAC2))
printn(resultat$)

IND_NB=t_ind_dep

centm=min/100
divis$=Str(min)+" "
divism=0
nbseq_1.q=SEQD\NBSEQ-1
nbseq=SEQD\NBSEQ
nbseqp=nbseq+1
dep_time.q=ElapsedMilliseconds()
pas=SEQD\TDIF(1)
;;;***********************************************************************************************************************
;;;******************************************  Début des 2 boucles pour rechercher  les NB premiers dans une zone MIN MAX
Repeat 
INDDD=1
pas=SEQD\TDIF(INDDD)
		
		While divis<=maxdiv; or FLAGFIN>0

		;   printn(_n(divis)+_n(pas)+_n(idivis))
	;     ;Le nom de @@: signifie étiquette anonyme,vous pouvez avoir défini un grand nombre d'entre eux dans la source.
	;     ;Symbole @b(ou l'équivalent @r)le plus proche des références du label précédent anonyme
	;     ;,symbole @f références plus proche de l'étiquette suivante anonyme.
	;     ;Ces symboles spéciaux sont insensibles à la casse.
	;;*******************************  TESTS EN ASM ********************************************************************
	;*************************************** Division d'un nombre de 64 bits par un nombre de 32 bits ******************************************
	EnableASM
	MOV ecx,dword[v_DIVIS]
	XOR EDX,EDX
	MOV eax,dword[v_MIN+4]
	DIV ecx
	;   MOV dword[v_quotientp],eax   ;;;; Le reste de la division nous suffit 
	MOV eax,dword[v_MIN]
	DIV ecx
	CMP edx,0
	JNZ @f
	INC dword [v_FLAGFIN] 
	Break
	!@@:
	MOV edx,dword[v_pas]
	ADD dword[v_DIVIS],edx  ; pas de PB l'addition est réalisée en valeur absolue 
													;     donc $FFFFFFFF =>4294967295 donc pas de risque de dépassement sqr(pow(2,63)-1))=>3037000499,9760496922867524030306 <4294967295
		;********************************************************************************************************************************************************
	;********************************************************************************************************************************************************
	; ;******************  peu efficace Pour accélerer encore l' algo   et surtout incompatible avec PB64bits*******************************
	
	
	
;   INC dword [v_INDDD]
;   mov ebx,dword[v_INDDD]
;   cmp ebx,dword[v_nbseqp]
; ; ;   jne @f
;   jne _FinFin
;   mov dword[v_INDDD],1
; ; ;   !@@:
;   !_FinFin:
;   LEA    ebp,[v_SEQD]
;   MOV    ebx,dword [v_INDDD]
;   MOV    ebp,dword [ebp+17]
;   MOVZX  eax,byte [ebp+ebx]
;   mov    dword[v_pas],eax

	DisableASM
	;********************************************************************************************************************************************************
	;********************************************************************************************************************************************************
	; ;***********************************  Pour accélerer encore l' algo  *******************************
	INDDD+1
	If INDDD=nbseqp  ;;; inutil de passer en ASM  voir ci dessus c'est plus long de près du double en temps sur la zone  9223372036854774807 1000
		INDDD=1 ;;;;;9223372036854770999
	EndIf
		pas=SEQD\TDIF(INDDD)
EnableASM
	;************************************  Pour accélerer encore l' algo vous pouvez décommentez les 5 lignes PB précedentes ******************************
	;*******************************************************************************************************************************************************
			Wend
		;;;; ***********************************************  FIN DE LA PREMI7RE BOUCLE *************************************************************
		If flagfin=0
;   debug _n(MIN)
				cent=min/100
				If centm<>cent 
						PrintN(result$)
						centm=cent 
										result$=Str(min)+" " 
				Else 
						result$+Str(min)+" " 
				EndIf
		Else
				flagfin=0
		EndIf  
		
		MIN+ seqd\TDIF(ind_nb)
;;;*******************************************************************************************************************************************************************
;;; On peut anticiper la prochaine racine différente de celle actuelle (A+1)^2 =A^2+1+2A si A=SQR(X)=RAC donc (RAC^2)=X alors => (RAC+1)^2= (RAC^2)+1+2*RAC=X+1+(2*RAC)
;;;*******************************************************************************************************************************************************************

		If MIN =>MINRAC2
			maxdiv.q=Sqr(MIN)
			MINRAC2.Q=MIN+(2*maxdiv)+1
;      PrintN("Changement de MAXDIV="+_q(maxdiv)+_q(minrac2)+_q(MIN)) ;;;; Pour vérifier le changement de racine carree
		EndIf  

						;************************************************ Recherche de la racine carré ASM  peut de gain ***************************************************************
;      enableasm 
;      !FILD qword[v_MIN]
;    !FSQRT
;     ; !FISTTP qword [v_MAXDIV] ; avec arrondi si vous décommentez cette instruction commentez la précédente
;     ;                           et laissez le ! car l'instruction FISTTP n'est pas reconnue par PureBasic avec !FISTTP c'est tout bon !!!!
;    !FISTP qword[v_MAXDIV] ; sans arrondi
;     disableasm 
		divis=SEQD\PREMDIV 
		ind_nb+1
;    if ind_nb>SEQD\NBSEQ-1
		If ind_nb>NBSEQ_1
				ind_nb=0  ;;; remise à zéro si le cycle des différences est atteint pour incrémentation des nombres susceptibles d'être premiers
		EndIf  
Until min>max Or min<0 ;;Si vous n'allez pas aux confins de la machine =$7FFFFFFFFFFFFFFF ;;; = 9223372036854775807 valeur max positive en ..q. 
;;                        Vous pouvez retirer (or min<0) ou dépassement de $7FFFFFFFFFFFFFFF
PrintN(result$)
PrintN( "Nombre Minimal = "+Str(MIN_DEP)+"  Nombre Maximal = "+Str(MAX))
PrintN("Nombre d'éléments du vecteur =:"+Str(SEQD\NBSEQ)+"  temps de préparation ="+Str(deb0-deb00))
PrintN("Temps millisecondes ="+Str(ElapsedMilliseconds()-dep_time))
Input() 
CloseConsole()
; debug " "

; ***************** Ci dessous les Nombres premiers pour les 1000 premiers nombres
; 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 547 557 563 569 571 577 587 593 599
; 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691
; 701 709 719 727 733 739 743 751 757 761 769 773 787 797
; 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887
; 907 911 919 929 937 941 947 953 967 971 977 983 991 997


A+
Il est fort peu probable que les mêmes causes ne produisent pas les mêmes effets.(Einstein)
Et en logique positive cela donne.
Il est très fortement probable que les mêmes causes produisent les mêmes effets.
Ollivier
Messages : 4190
Inscription : ven. 29/juin/2007 17:50
Localisation : Encore ?
Contact :

Re: A propos de nombres premiers

Message par Ollivier »

@PAPIPP

Bonjour à toi PAPIPP, n'hésite pas à ajouter une intro, une synthèse ou une conclusion à ton travail pour situer le lecteur : par exemple, cette histoire de racine concerne la méthode par division.

Aussi ta recherche d'anticipation du choix de la racine permet de gagner en vitesse en éludant l'ensemble des tests diviseurs à l'ensemble des premiers uniquement.

Cette anticipation permet d'éviter environ un million de tests inutiles pour moins d'un millier de modifications effectives dans un extrait d'environ 2 à 3 millions de nombres testés.

Conclusion : ta technique d'anticipation est indispensable dans le long chemin de la sophistication des tests par division.

Merci donc. Je vais tâcher d'apporter un code source promis de date assez conséquente sur le travail initial de fweil et la manière de comprendre son astuce aussi subtile que rapide. Ce n'est plus qu'une question d'heures pour pouvoir publier. Après, j'attaquerai mon application sur la compression des gaps (il y en a pour quelques semaines...).
Ollivier
Messages : 4190
Inscription : ven. 29/juin/2007 17:50
Localisation : Encore ?
Contact :

Re: A propos de nombres premiers

Message par Ollivier »

Je n'arrive pas à représenter ce que je veux : si fweil avait pu m'expliquer pourquoi 2*N?

J'arrive à comprendre N * N : tout a un ordre et une logique, par contre, avec N * 2, c'est "vaseux".

J'obtiens un algo à 10^9 nombres testés par seconde, non optimisé avec N * N. Ça me désoeuvre de tenter de démontrer N * 2 qui est forcément plus lent...

@PAPIPP

Bon, désolé, mais j'ai quand même une question : quand tu veux vérifier un nombre de 18 chiffres, avec 80 millions de divisions, au lieu de juste conclure :
"tous les restants sont non nuls, donc c'est un premier"
ou
"au moins un restant est nul, ce n'est pas un premier."

Est-ce que tu pourrais rajouter toutes les décisions que chaque restant non nul apporte au lieu de les jeter à la corbeille?

Un pico-exemple :
41 indique après test par division qu'il est premier certes, mais, en plus, il informe que 42 et 45 sont non premiers, non?
PAPIPP
Messages : 534
Inscription : sam. 23/févr./2008 17:58

Re: A propos de nombres premiers

Message par PAPIPP »

Bonjour Ollivier
@PAPIPP

Bon, désolé, mais j'ai quand même une question : quand tu veux vérifier un nombre de 18 chiffres, avec 80 millions de divisions, au lieu de juste conclure :
"tous les restants sont non nuls, donc c'est un premier"
ou
"au moins un restant est nul, ce n'est pas un premier."

Est-ce que tu pourrais rajouter toutes les décisions que chaque restant non nul apporte au lieu de les jeter à la corbeille?

Le prg comporte 2 boucles.
1) la plus externe incrémente les Nb à tester suivant la table des différences.
2) La plus interne passe tous les diviseurs possibles jusqu’à trouver un reste nul et de jeter ce NB de la boucle externe à la poubelle.
3) On incrémente à nouveau le Nb à tester suivant la table des différences et on repasse au paragraphe 2) en recommençant à tester tous les diviseurs en utilisant là aussi la table des différences (Les indices d'incrémentation dans la table des différences sont propres à chaque boucle)

Comme le prg recherche les NB premiers et tout nombre qui n’est pas premier est jeté.

Par contre j’ai réalisé un prg de décomposition d’un seul nombre en facteur premier.
On peut l’étendre sur une zone ce qui n’est pas le cas actuellement.
Ce prg est basé sur les mêmes préparations que le prg de recherche des Nb premiers
Il faut donc le placer dans le même répertoire que les 2 autres prg

Voici ce prg de décomposition d’un nombre <= 9223372036854775807
résultat de la décomposition de 9223372036854770999 => 26 secondes alors que la décomposition de 9223372036854775807 => 2ms

Code : Tout sélectionner

9223372036854770999=2807010931*3285834029
9223372036854770999: contrôle
temps=26182ms
Recherche des diviseurs <=:23 + diviseurs =>23
Nombre d'éléments du vecteur =:36495360  temps de préparation =56
9223372036854775807=7*7*73*127*337*92737*649657
9223372036854775807: contrôle
temps=2ms
Recherche des diviseurs <=:23 + diviseurs =>23
Nombre d'éléments du vecteur =:36495360 temps de préparation =56


Code : Tout sélectionner

Macro _q_t_ 
		" 
EndMacro 
Macro _n (__n) 
_q_t_#__n#=_q_t_+Str(__n)+" " 
EndMacro 
Macro _Q (__Q) 
_q_t_#__Q#=_q_t_+Str(__Q)+" " 
EndMacro 
; exemple des p_fact(x)
; p_fact(x)--modulo-----------------Nb_occurrences des différences--------Taux-----------Gain
; 3__________6__________________________2_________________________________0,3333_________66,66%
; 5__________30_________________________8_________________________________0,2666_________73,33%
; 7__________210________________________48________________________________0,2285_________77,14%
; 11_________2310_______________________480_______________________________0,2078_________79,22%
; 13_________30030______________________5760______________________________0,1918_________80,81%
; 17_________510510 ____________________92160_____________________________0,18052________81,94%
; 19_________9699690____________________1658880___________________________0,17102________82,89%
; 23_________223092870__________________36495360__________________________0,16358________83,64%
; ind
EnableExplicit
Structure DIVIS
  NBPFACT.a
  PREMDIV.l
  NBSEQ.l
  MODULO.q
  ;   Array TDIF.a(40000000)  ;; *** ATTENTION   A utiliser avec DATA.A et read.A  mais très lent avec p_fact(23) P_fact(19) est préférable remplacer 40 par 2
  Array TDIF.a(40000000)  
EndStructure
Define nb$,nb.q,l2nb.d,nbl2.q,DIVMAX.q,NBREC.q,DIVIS.q,i,ipas,fact_prem$,j,deb1.q,mess$,quotientp.q,iprem,indt,Restd.q,ind.l,pas,idivis,FICHIER$
Global nbs.l,NBSEQ.l,SEQD.DIVIS
Dim Tab.q(63) ;;; max = pow(2,63)-1 donc 63 éléments au max
Define *Tab=@Tab(), ADRDSEQ.q,AdresdebF.q,AdresFinF.q,Adresdeb.q,delta
             ; 2 3 5 7 11 13 17 19 23
Dim Nbprem(10)
Nbprem(0)=2:Nbprem(1)=3:Nbprem(2)=5:Nbprem(3)=7:Nbprem(4)=11:Nbprem(5)=13:Nbprem(6)=17:Nbprem(7)=19:Nbprem(8)=23:Nbprem(9)=29
OpenConsole()
Procedure.s choixdiv()
  Protected mess$, nb$, ELEM, ipas, i,COUR_DIR$,Filtre$,fichier$,Tfichier,enreg1$,enreg2$,pos,restenreg2$,flag,nb,jj,restenreg$,fichierp$,AdredebF,Adredeb,TTAB,RTTAB,AdresDEBA,decoupe,rest
  Protected EOFF,leftfich$,RIGHTfich$
  Global NBtotal=0, NBTOTNB=0,deb00.q,deb0.q
  COUR_DIR$ = GetCurrentDirectory()
  Filtre$ = "MEMA (MEMA_*.BIN)|MEMA_*.BIN;|Tous les fichiers (*.*)|*.*"
  fichier$=OpenFileRequester("Choisissez un fichier MEMA ou annulez", COUR_DIR$+"\MEMA_17.BIN", Filtre$, 0)
  fichierp$=ucase(trim(GetFilePart(fichier$)))
  Tfichier=FileSize(Fichier$)
  leftfich$=left(fichierp$,5)
  RIGHTfich$=right(fichierp$,4)
  If Tfichier<1 or left(fichierp$,5)<>"MEMA_" or right(fichierp$,4) <>".BIN"
    MessageRequester( "ATTENTION", "fichier vide ou le nom n'est pas conforme "+ _n(Tfichier)+"  "+fichierp$)
    End
  EndIf
  ProcedureReturn fichier$
EndProcedure
FICHIER$=choixdiv()
deb00.q=ElapsedMilliseconds()
;************************************************************************************************************************************
ADRDSEQ.q=@SEQD
AdresdebF.q=@SEQD\NBPFACT
AdresFinF.q=@SEQD\MODULO+sizeof(quad)
Adresdeb.q=@SEQD\TDIF(0)
delta=AdresFinF-AdresdebF

if openfile(2,FICHIER$);,#PB_File_SharedRead|#PB_File_NoBuffering)
  readdata(2,ADRDSEQ ,delta)
  REDIM SEQD\TDIF(SEQD\NBSEQ+1)
  readdata(2, Adresdeb,SEQD\NBSEQ+1)
  closefile(2)
endif

deb0.q=ElapsedMilliseconds()

saisie:
mess$+"Décomposition en facteur premier"
If Len(mess$)>99
  MessageRequester("Erreur","Donnez un entier positif < ou  9223372036854775807"+#CR$+"Relancez le prg") ;
  End
EndIf
nb$=InputRequester(mess$,"Donnez un entier positif < ou =9223372036854770999","9223372036854775807") ; c'est la limite d'un type q  2^63-1
                                                                                                      ;;                                                                                   en puissance de 2=>2^62 donne le plus grand nombre de facteurs
If Len(nb$)>Len("9223372036854775807")                                                                ;; en puissance de 3 => 3^39
  Goto saisie
EndIf
nb.q=Val(nb$)
If nb<1
  Goto saisie
EndIf
;*********************************************************************************************************************************************************
deb1.q=ElapsedMilliseconds()

DIVMAX.Q=Sqr(NB)      ;limite de recherche des facteurs premiers
NBREC.Q=1
DIVIS=2
i=0
fact_prem$=nb$+"="
deb1.q=ElapsedMilliseconds()

;************************************** Recherche des diviseurs premiers qui se trouvent dans p_fact(x) ***********************
DIVMAX=Sqr(nb)
iprem=0
DIVIS=nbprem(iprem)
indt=0
While DIVIS<=SEQD\NBPFACT And DIVIS<=DIVMAX 
  Restd.q=nb%DIVIS
  If Restd=0
    tab(indt)=DIVIS
    indt+1
    quotientp=nb/DIVIS
    nb=quotientp
    DIVMAX=Sqr(nb)
  Else 
    iprem+1
    DIVIS=nbprem(iprem)
  EndIf  
Wend
;****************  Partie principale pour le reste de la décomposition **************************************
If DIVIS<=seqd\PREMDIV
  DIVIS=SEQD\PREMDIV
EndIf
ipas=0
;************************************************************************************************************************
pas=SEQD\TDIF(1)
idivis=1
NBSEQ=SEQD\NBSEQ
;   ind=indt*8
Repeat
;   printn(_n(divis)+_n(pas)+_n(idivis))
  ;     ;Le nom de @@: signifie étiquette anonyme,vous pouvez avoir défini un grand nombre d'entre eux dans la source.
  ;     ;Symbole @b(ou l'équivalent @r)le plus proche des références du label précédent anonyme
  ;     ;,symbole @f références plus proche de l'étiquette suivante anonyme.
  ;     ;Ces symboles spéciaux sont insensibles à la casse.
  ;;*******************************  TESTS EN ASM ********************************************************************
  ;*************************************** Division d'un nombre de 64 bits par un nombre de 32 bits ******************************************
  EnableASM
  MOV ecx,dword[v_DIVIS]  ;;;; la division euclidienne X=kQ+R avec R<Q nous oblige à utiliser une astuce
  XOR EDX,EDX             ;;;; pour éviter avec un diviseur trop petit et un dividente trop grand d'avoir un reste de type qword or EDX est de type dword 
  MOV eax,dword[v_nb+4]
  DIV ecx
  MOV dword[v_quotientp],eax
  MOV eax,dword[v_nb]
  DIV ecx
  CMP edx,0
  JNZ @f
  ;     MOV dword[v_reste],edx
  ; quotientp

  MOV edx,dword[v_quotientp]
  MOV dword[v_nb+4],edx
  MOV dword[v_nb],eax
  ;************************************************ Recherche de la racine carré ***************************************************************
  FILD qword[v_nb]
  FSQRT
  ; !FISTTP dword [v_DIVMAX] ; avec arrondi si vous décommentez cette instruction commentez la précédente
  ;                           et laissez le ! car l'instruction FISTTP n'est pas reconnue par PureBasic avec !FISTTP c'est bon !!!!
  FISTP qword[v_DIVMAX] ; sans arrondi
                        ;                          ;*********************************************** Affectation de la valeur du DIViseur dans la table *******************************************
                        ;     MOV ecx, [v_DIVIS]; à decommenter si ecx est effacée
                        ;     mov ebx, [p_Tab]
                        ;     mov edx, [v_ind]
                        ;     add ebx,edx
                        ;     mov[ebx],ecx
                        ;     add edx,8
                        ;     mov[v_ind],edx
  DisableASM
  tab(indt)=DIVIS   ;;;; Plus rapide que l'option ASM
  indt+1
  If nb<=1:Break:EndIf
  EnableASM
  ;     JMP lab#MacroExpandedCount
  JMP lab
  !@@:
  MOV edx,[v_pas]
  ADD dword[v_DIVIS],edx  ; pas de PB l'addition est réalisée en valeur absolue 
                          ;     donc $FFFFFFFF =>4294967295 donc pas de risque de dépassement sqr(pow(2,63)-1))=>3037000499,9760496922867524030306 <4294967295
  DisableASM
  ;********************************************************************************************************************************************************
  ;********************************************************************************************************************************************************
  ; ;***********************************  Pour accélerer encore l' algo  *******************************
  idivis+1
  If idivis=NBSEQ+1  ;;; inutil de passer en ASM  voir ci dessous c'est plus long de 5 à 6 secondes  sur la recherche de 9223372036854770999
    idivis=1
  EndIf
    pas=SEQD\TDIF(idivis)
EnableASM
  ;     mov EAX, dword [v_idivis]
  ;     xor edx,edx
  ;     mov ECX, dword [v_NBSEQ]
  ;     div ECX 
  ;     mov  dword [v_idivis],EDX
  ;     idivis=idivis%NBSEQ
  ;************************************  Pour accélerer encore l' algo vous pouvez décommentez les 7 lignes PB précedentes ******************************
  ;*******************************************************************************************************************************************************
  ;     !lab#MacroExpandedCount:
  !lab:
  If DIVIS>DIVMAX Or DIVIS<0:Break:EndIf
ForEver


;*************************************************Fin macro decompose2****************************************************************************************************
; indt=ind/8
If nb >1
  tab(indt)=NB
EndIf   

;************************************   EDITION ***********************************************
For j=0 To indt
  fact_prem$+tab(j)
  NBREC*tab(j) ; Vérification de la décomposition
  If j<indt
    fact_prem$+"*"
  EndIf
Next
;essai pour ces nombres 9223372036854769861, 9223372036854770317, 9223372036854770999, 9223372036854773927 et 9223372036854775309 qui ne sont pas premiers.

PrintN(fact_prem$)
PrintN(Str(NBREC)+": contrôle")
PrintN("temps="+Str(ElapsedMilliseconds()-deb1)+"ms")
PrintN("Recherche des diviseurs <=:"+Str(SEQD\NBPFACT)+" + diviseurs =>"+Str(SEQD\NBPFACT))
PrintN("Nombre d'éléments du vecteur =:"+Str(SEQD\NBSEQ)+"  temps de préparation ="+Str(deb0-deb00))
MessageRequester("Résultat",fact_prem$+#CR$+Str(NBREC)+": contrôle"+#CR$+"temps="+Str(ElapsedMilliseconds()-deb1)+"ms"+#CR$+"Recherche des diviseurs <=:"+Str(SEQD\NBPFACT)+" + diviseurs =>"+Str(SEQD\NBPFACT)+#CR$+"Nombre d'éléments du vecteur =:"+Str(SEQD\NBSEQ))
Input()
CloseConsole()
; ces nombres 9223372036854769861, 9223372036854770317, 9223372036854770999, 9223372036854773927 et 9223372036854775309 qui ne sont pas premiers; 2295911257* 4017303373=9223372036854769861;==> div=-2147483645 divm=2147483647 DIVMAX=3037000500
; ***************** Ci dessous les Nombres premiers pour les 1000 premiers nombres
; 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 547 557 563 569 571 577 587 593 599
; 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691
; 701 709 719 727 733 739 743 751 757 761 769 773 787 797
; 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887
; 907 911 919 929 937 941 947 953 967 971 977 983 991 997



A+
Il est fort peu probable que les mêmes causes ne produisent pas les mêmes effets.(Einstein)
Et en logique positive cela donne.
Il est très fortement probable que les mêmes causes produisent les mêmes effets.
Ollivier
Messages : 4190
Inscription : ven. 29/juin/2007 17:50
Localisation : Encore ?
Contact :

Re: A propos de nombres premiers

Message par Ollivier »

Un petit code pour comprendre le fonctionnement de l'astuce de fWeil:

Code : Tout sélectionner

 
;- Compiler : PB 5.50 X64
;-

Structure Graphic
        Width.I
        Height.I
        Depth.I
        Array CharFont.I(80)
        Array CharS.I(80, 255)
        Array CharW.I(80, 255)
        Array CharH.I(80, 255)
        MouseImage.I
        MouseSprite.I
EndStructure

Procedure GraphicCreate()
        Define *This.Graphic = AllocateMemory(SizeOf(Graphic) )
        InitializeStructure(*This, Graphic)
        ProcedureReturn *This
EndProcedure

Procedure GraphicOpen(*This.Graphic)
        With *This
                ExamineDesktops()
                \Width = DesktopWidth(0)
                \Height = DesktopHeight(0)
                \Depth = DesktopDepth(0)
                InitSprite()
                InitKeyboard()
                InitMouse()
                OpenScreen(\Width, \Height, \Depth, "")
        EndWith
EndProcedure

Procedure GraphicInitMouse(*This.Graphic)
        
        With *This
                \MouseImage = CreateImage(#PB_Any, 32, 32, 32)
                If StartDrawing(ImageOutput(\MouseImage) )
                        DrawingMode(#PB_2DDrawing_AllChannels)
                        Box(0, 0, 32, 32, RGBA(0, 0, 0, 0) )
                        StopDrawing()
                        If StartVectorDrawing(ImageVectorOutput(\MouseImage) )
                                ScaleCoordinates(0.3, 0.3)
                                MovePathCursor(2, 4)
                                AddPathLine(0, 60, #PB_Path_Relative)
                                AddPathLine(15, -15, #PB_Path_Relative)
                                AddPathLine(10, 20, #PB_Path_Relative)
                                AddPathLine(10, -5, #PB_Path_Relative)
                                AddPathLine(-10, -20, #PB_Path_Relative)
                                AddPathLine(20, 0, #PB_Path_Relative)
                                ClosePath()
                                VectorSourceColor(RGBA(255, 255, 255, 255))
                                FillPath(#PB_Path_Preserve)
                                VectorSourceColor(RGBA(0, 0, 0, 255))
                                StrokePath(4)
                                StopVectorDrawing()
                        EndIf
                EndIf
                
                \MouseSprite = CreateSprite(#PB_Any, 32, 32, #PB_Sprite_AlphaBlending)
                If StartDrawing(SpriteOutput(\MouseSprite) )
                        DrawingMode(#PB_2DDrawing_AllChannels)
                        Box(0, 0, 32, 32, RGBA(0, 0, 0, 0) )
                        DrawingMode(#PB_2DDrawing_AlphaBlend)
                        DrawImage(ImageID(\MouseImage), 0, 0)
                        StopDrawing()
                EndIf
        EndWith
        
EndProcedure

Procedure GraphicInitFont(*This.Graphic, FontName.S, SizeList.S)
        Define SizeSize = CountString(SizeList, ";")
        Dim Size.I(SizeSize)
        For K = 0 To SizeSize
                Size(K) = Val(StringField(SizeList, 1 + K, ";") )
        Next
        For K = 0 To SizeSize
                J = Size(K)
                *This\CharFont(J) = LoadFont(#PB_Any, FontName, J)
        Next
        
        StartDrawing(ScreenOutput() )
        For K = 0 To SizeSize
                J = Size(K)
                DrawingFont(FontID(*This\CharFont(J) ) )
                For I = 1 To 255
                        *This\CharW(J, I) = TextWidth(Chr(I) )
                        If *This\CharW(J, I) = 0
                                *This\CharW(J, I) = 1
                        EndIf
                        *This\CharH(J, I) = TextHeight(Chr(I) )
                Next
        Next
        StopDrawing()
        
        For K = 0 To SizeSize
                J = Size(K)
                For I = 1 To 255
                        *This\CharS(J, I) = CreateSprite(#PB_Any, *This\CharW(J, I), *This\CharH(J, I), #PB_Sprite_AlphaBlending)
                        StartDrawing(SpriteOutput(*This\CharS(J, I) ) )
                        DrawingFont(FontID(*This\CharFont(J) ) )
                        DrawingMode(#PB_2DDrawing_AllChannels)
                        DrawText(0, 0, Chr(I), RGBA(0, 0, 0, 255), RGBA(255, 255, 255, 0) )
                        StopDrawing()
                Next
        Next
EndProcedure

Procedure GraphicDisplay(*This.Graphic, x, y, A$, Size, ColorAlpha, Angle = 0)
        Alpha = Alpha(ColorAlpha)
        Color = RGB(Red(ColorAlpha), Green(ColorAlpha), Blue(ColorAlpha) )
        Define AngleRad.F = Angle * #PI / 180.0
        If Angle
                For I = 1 To Len(A$)
                        A = Asc(Mid(A$, I, 1) )
                        If Angle
                                RotateSprite(*This\CharS(Size, A), Angle, #PB_Absolute)
                        EndIf
                        DisplayTransparentSprite(*This\CharS(Size, A), X, Y, Alpha, Color)
                        X + (*This\CharW(Size, A) * Cos(AngleRad) )
                        Y + (*This\CharW(Size, A) * Sin(AngleRad) )
                Next
        Else
                For I = 1 To Len(A$)
                        A = Asc(Mid(A$, I, 1) )
                        RotateSprite(*This\CharS(Size, A), 0, #PB_Absolute)
                        DisplayTransparentSprite(*This\CharS(Size, A), X, Y, Alpha, Color)
                        X + *This\CharW(Size, A)
                Next
        EndIf
EndProcedure

;-

        Global *Gr.Graphic = GraphicCreate()
        GraphicOpen(*Gr)
        GraphicInitFont(*Gr, "Verdana", "10;14;16")
        GraphicInitMouse(*Gr)
        
Structure Box
        ImageNum.I[3]
        SpriteNum.I[3]
EndStructure         

Procedure BoxCreate()
        Define *This.Box = AllocateMemory(SizeOf(Box) )
        ProcedureReturn *This
EndProcedure

Procedure BoxInit(*This.Box, Index.I, Width.I, Height.I, Color.I)
        Define W, H
        With *This
                \ImageNum[Index] = CreateImage(#PB_Any, Width, Height, 32)
                If StartDrawing(ImageOutput(\ImageNum[Index] ) )
                        DrawingMode(#PB_2DDrawing_AllChannels)
                        W = OutputWidth()
                        H = OutputHeight()
                        Box(0, 0, W, H, RGBA(0, 0, 0, 255) )
                        Box(1, 1, W - 2, H - 2, Color)
                        StopDrawing()
                EndIf
                \SpriteNum[Index] = CreateSprite(#PB_Any, Width, Height, #PB_Sprite_AlphaBlending)
                If StartDrawing(SpriteOutput(\SpriteNum[Index] ) )
                        W = OutputWidth()
                        H = OutputHeight()
                        DrawingMode(#PB_2DDrawing_AllChannels)
                        Box(0, 0, W, H, RGBA(0, 0, 0, 0) )
                        DrawingMode(#PB_2DDrawing_AlphaBlend)
                        DrawImage(ImageID(\ImageNum[Index] ), 0, 0)
                        StopDrawing()
                EndIf
        EndWith
        
EndProcedure

Procedure.I HarmoDisplay(*Box.Box, X, Y, Column, Modulo, N)
        Define R
        H = Modulo - (N % Modulo)
        If H = 1
                DisplayTransparentSprite(*Box\SpriteNum[1], X + 16 * Column, Y + 16 * H)
                R = 1
        Else
                DisplayTransparentSprite(*Box\SpriteNum[2], X + 16 * Column, Y + 16 * H)
                R = 0
        EndIf
        GraphicDisplay(*Gr, X + 16 * Column, Y - 16, Str(Modulo), 10, RGBA(255, 255, 0, 127), 270)
        ProcedureReturn R
EndProcedure

        Define *Box.Box = BoxCreate()
        BoxInit(*Box, 1, 16, 16, RGBA(255, 255, 255, 255) )
        BoxInit(*Box, 2, 16, 16, RGBA(255, 0, 0, 255) )
        
        Dim P(1 << 16 - 1)
        MaxP = 0
        Repeat
                Delay(16)
                If IsScreenActive()
                        ExamineKeyboard()
                        ExamineMouse()
                        ClearScreen(RGB(63, 159, 255) )
                        MouseW = MouseWheel()
                        MouseX = MouseX()
                        MouseY = MouseY()
                        MouseBForward = MouseB
                        MouseB = MouseButton(1)
                        If MouseB <> MouseBForward
                                If MouseB
                                        Auto ! 1
                                EndIf
                        EndIf
                        GraphicDisplay(*Gr, 0, 0, "Press [Escape] key to 'kit'... Houiler la molette pour +1 / -1... Cliquez pour Play/Stop...", 14, RGBA(255, 255, 0, 127) )
                        If Auto
                                N + 1
                        Else
                                N - MouseW
                                If N < 1
                                        N = 1
                                EndIf
                        EndIf
                        GraphicDisplay(*Gr, 0, 20, Str(N + 1), 14, RGBA(255, 255, 0, 127) )
                        Etat = 0
                        For I = 1 To MaxP
                                Etat + HarmoDisplay(*Box, 100, 100, I, P(I), N)
                        Next
                        If Etat = 0
                                MaxP + 1
                                P(MaxP) = N + 1
                        EndIf
                        DisplayTransparentSprite(*Gr\MouseSprite, MouseX, MouseY)
                        FlipBuffers()
                EndIf
        Until KeyboardPushed(#PB_Key_Escape)
PAPIPP
Messages : 534
Inscription : sam. 23/févr./2008 17:58

Re: A propos de nombres premiers

Message par PAPIPP »

Bonjour à tous

Avant de m’attaquer à la recherche des nombres premiers entre 1 et (2^64)-2, j’ai préféré me frotter à la décomposition d’un nombre compris dans cette plage,
car dans ce prg je n’ai à extraire la racine d’un nombre >(2^63)-1 qu’une seule fois au départ car un diviseur du nombre à décomposer réduit le quotient à un entier positif
qui rentre dans une variable X.q alors que dans le prg de recherche il faut vérifier ou calculer la racine à chaque nombre testé.

Voici donc le programme qui décompose un nombre compris entre 1 et (2^64)-2
Le précédent prg décomposait entre 1 et (2^63)-1
Toujours avec les mêmes contraintes de fonctionnement sous 32 ou 64 bits sans avoir à recompiler ou en utilisant les mêmes instructions. Et j’ai essayé de conserver la vitesse d’exécution.

Pourquoi rechercher la décomposition des nombres entre (2^63) et (2^64] ?
Il y a autant de nombres entre 1 et (2^63] qu’entre (2^63) et (2^64)
En effet (2^64]=(2^63)*2= (2^63)+(2^63) or les nombres entre (2^-63) et (2^64)
NB_nombre=(2^64)-(2^63) = (2^63)+(2^63)-(2^63)= (2^63) CQFD
Il y a donc (2^63) nombres entre (2^63) et (2^64)
J’ai donc doublé les possibilités du prg.

Les principaux problèmes rencontrés sont les difficultés à travailler avec des nombres signés.
En effet tous les nombres introduits dans une variable de type .q (quad) sont négatifs à partir de $8000000000000000. (2^63]
Ainsi par exemple prendre la racine carrée d’un nombre négatif n’est pas très facile SQR(-1)=i .
Or nous considérons tous les nombres => $8000000000000000 comme entiers et positifs.
Il faut donc utiliser des instructions ASM qui ne tiennent pas compte du signe.
Par contre extraire la racine carrée d’un nombre négatif n’est toujours pas résolu avec les instructions en 32 bits.

Voici ce que je propose.
$FFFFFFFFFFFFFFFF= (2^64)-1 et sqr($FFFFFFFFFFFFFFFF)=$FFFFFFFF=4294967295
On voit que le diviseur max que l’on aura à manipuler =(2^32)-1 qui est positif dans une variable de structure (quad) X.q
D’autre part la valeur max possible d’être conserver positive dans un X.q est (2^63)-1
Alors si NBmax>(2^63)-1 donc > 9223372036854775807 avec LIMITQ.q=$7FFFFFFFFFFFFFFF
Nous pouvons travailler avec la différence entre NBmax et (2^63)-1
DIFF.q=NBmax-LIMITQ pour simplifier appelons DIFF=D NBmax=A et LIMITQ=L
A=L+D on peut écrire A=L+KL avec KL=D donc k=D/L
A=L+KL => L(1+k) => L(1+(D/L))
Et SQR(A) = SQR(L)*(sqr(1+k) => SQR(L)*(SQR(1+(D/L)) or cette équation conserve tous les termes réels positifs facilement convertible en entiers positifs

On voit que k varie dans une plage de -1 à +1
et D dans une plage -=$7FFFFFFFFFFFFFFF à +$7FFFFFFFFFFFFFFF
Ceci nous donne le nombre max possible que l’on peut traité =
L+D_max= $7FFFFFFFFFFFFFFF +$7FFFFFFFFFFFFFFF=$FFFFFFFFFFFFFFFE cette valeur est la valeur max que je peux traiter
Si vous introduisez $FFFFFFFFFFFFFFFF le prg refusera et vous demandera une autre valeur.
Par exemple $FFFFFFFFFFFFFFFE donne comme décomposition

$FFFFFFFFFFFFFFFE=2*7*7*73*127*337*92737*649657
$FFFFFFFFFFFFFFFE: contrôle
temps=2ms
Vous pouvez vérifier avec la calculette de Microsoft qui donne des résultats fiables

Code : Tout sélectionner

 Macro _q_t_ 
		" 
 EndMacro 
 Macro _n (__n) 
 _q_t_#__n#=_q_t_+Str(__n)+" " 
 EndMacro 
 Macro _Q (__Q) 
 _q_t_#__Q#=_q_t_+Str(__Q)+" " 
 EndMacro 
; exemple des p_fact(x)
; p_fact(x)--modulo-----------------Nb_occurrences des différences--------Taux-----------Gain
; 3__________6__________________________2_________________________________0,3333_________66,66%
; 5__________30_________________________8_________________________________0,2666_________73,33%
; 7__________210________________________48________________________________0,2285_________77,14%
; 11_________2310_______________________480_______________________________0,2078_________79,22%
; 13_________30030______________________5760______________________________0,1918_________80,81%
; 17_________510510 ____________________92160_____________________________0,18052________81,94%
; 19_________9699690____________________1658880___________________________0,17102________82,89%
; 23_________223092870__________________36495360__________________________0,16358________83,64%
; ind
EnableExplicit

Procedure.s HEX2DEC(HEXS.S)
  EnableExplicit
  Static Dim t_rest.a(1000)
  Protected DIVIDENDE$,vald.q,longt,resultat$,DIVIDP$,DIVIDPR$,QUOTP.q,RESTP.q,quotp$,j,ii,DIVIDP.Q,Irest
  DIVIDENDE$=LTrim(UCase(HEXS))
  vald.q=Val("$"+DIVIDENDE$)
  longt=Len(DIVIDENDE$)
  If vald>0 And longt<17
    resultat$=Str(vald)
  Else  
    Irest=0
    DIVIDP$=""
    DIVIDPR$=""
    quotp$=""
    Repeat 
      For ii=1 To longt
        DIVIDP$=DIVIDPR$+Mid(DIVIDENDE$,ii,1)
        DIVIDP.Q=Val("$"+DIVIDP$)
        ;   QUOTP.q=DIVIDP/10
        ;   RESTP.q=DIVIDP%10
        EnableASM
        MOV ax,word [p.v_DIVIDP]  
        MOV cl,10
        ;   idiv Cl 
        DIV cl  ; utilise moins de cycles machine que idiv
        MOV  [p.v_QUOTP],al
        MOV [p.v_RESTP],AH
        DisableASM
        
        DIVIDPR$=Hex(RESTP)
        quotp$+Hex(QUOTP)
      Next 
      t_rest(Irest)=RESTP
      Irest+1 
      DIVIDENDE$=QUOTP$
      longt=Len(DIVIDENDE$) 
      DIVIDP$=""
      DIVIDPR$=""
      quotp$=""
      
    Until Val("$"+ dividende$)=0
    For j=Irest-1 To 0 Step-1
      resultat$+Str( t_rest(j))
      t_rest(j)=0
    Next
  EndIf
  ProcedureReturn  resultat$
  DisableExplicit
EndProcedure

Procedure.s DEC2HEX(DECI.S)
  EnableExplicit
  Static Dim t_rest.a(1000)
  Protected DIVIDENDE$,vald.q,longt,resultat$,DIVIDP$,DIVIDPR$,QUOTP.q,RESTP.q,quotp$,i,j,ii,DIVIDP.Q,Irest
  DIVIDENDE$=UCase(DECI)
  longt=Len(DIVIDENDE$)
  If Val(DIVIDENDE$)=>0  And longt<20
    vald.q=Val(DIVIDENDE$)
    resultat$=Hex(vald)
  Else  
    irest=0
    DIVIDP$=""
    DIVIDPR$=""
    quotp$=""
    Repeat 
      For ii=1 To longt
        DIVIDP$=DIVIDPR$+Mid(DIVIDENDE$,ii,1)
        DIVIDP.Q=Val(DIVIDP$)
        ;   QUOTP.q=DIVIDP/16
        ;   RESTP.q=DIVIDP%16
        EnableASM
        MOV ax,word [p.v_DIVIDP]  
        MOV cl,16
        ;   div Cl 
        DIV Cl  ; utilise moins de cycles machine que idiv
        MOV  [p.v_QUOTP],al
        MOV [p.v_RESTP],AH
        DisableASM
        
        DIVIDPR$=Str(RESTP)
        quotp$+Str(QUOTP)
      Next 
      t_rest(Irest)=RESTP
      Irest+1 
      DIVIDENDE$=QUOTP$
      longt=Len(DIVIDENDE$) 
      DIVIDP$=""
      DIVIDPR$=""
      quotp$=""
      
    Until Val(dividende$)=0
    For j=Irest-1 To 0 Step-1
      resultat$+Hex( t_rest(j))
      t_rest(j)=0
    Next
  EndIf
  ProcedureReturn  resultat$
  DisableExplicit
EndProcedure


Structure DIVIS
	NBPFACT.a
	PREMDIV.l
	NBSEQ.l
	MODULO.q
	Array TDIF.a(40000000)  
EndStructure
Define nb$,nb.q,l2nb.d,nbl2.q,DIVMAX.q,NBREC.q,DIVIS.q,i,ipas,fact_prem$,j,deb1.q,mess$,quotientp.q,iprem,indt,Restd.q,ind.l,pas,idivis,FICHIER$
Global nbs.l,NBSEQ.l,SEQD.DIVIS,LIMIT$="7FFFFFFFFFFFFFFF",LIMITQ.Q=$7FFFFFFFFFFFFFFF,RACLIM.d=SQR(LIMITQ)
Dim Tab.q(65) ;;; max = pow(2,64)-1 donc 64+1 éléments au max
Define *Tab=@Tab(), ADRDSEQ.q,AdresdebF.q,AdresFinF.q,Adresdeb.q,delta,DIVMAXD.d
							; 2 3 5 7 11 13 17 19 23
Dim Nbprem(10)
Nbprem(0)=2:Nbprem(1)=3:Nbprem(2)=5:Nbprem(3)=7:Nbprem(4)=11:Nbprem(5)=13:Nbprem(6)=17:Nbprem(7)=19:Nbprem(8)=23:Nbprem(9)=29
OpenConsole()
Procedure.s choixdiv()
	Protected mess$, nb$, ELEM, ipas, i,COUR_DIR$,Filtre$,fichier$,Tfichier,enreg1$,enreg2$,pos,restenreg2$,flag,nb,jj,restenreg$,fichierp$,AdredebF,Adredeb,TTAB,RTTAB,AdresDEBA,decoupe,rest
	Protected EOFF,leftfich$,RIGHTfich$
	Global NBtotal=0, NBTOTNB=0,deb00.q,deb0.q
	COUR_DIR$ = GetCurrentDirectory()
	Filtre$ = "MEMA (MEMA_*.BIN)|MEMA_*.BIN;|Tous les fichiers (*.*)|*.*"
	fichier$=OpenFileRequester("Choisissez un fichier MEMA ou annulez", COUR_DIR$+"\MEMA_17.BIN", Filtre$, 0)
	fichierp$=ucase(trim(GetFilePart(fichier$)))
	Tfichier=FileSize(Fichier$)
	leftfich$=left(fichierp$,5)
	RIGHTfich$=right(fichierp$,4)
	If Tfichier<1 or left(fichierp$,5)<>"MEMA_" or right(fichierp$,4) <>".BIN"
		MessageRequester( "ATTENTION", "fichier vide ou le nom n'est pas conforme "+ _n(Tfichier)+"  "+fichierp$)
		End
	EndIf
	ProcedureReturn fichier$
EndProcedure
FICHIER$=choixdiv()
deb00.q=ElapsedMilliseconds()
;************************************************************************************************************************************
ADRDSEQ.q=@SEQD
AdresdebF.q=@SEQD\NBPFACT
AdresFinF.q=@SEQD\MODULO+sizeof(quad)
Adresdeb.q=@SEQD\TDIF(0)
delta=AdresFinF-AdresdebF

if openfile(2,FICHIER$);,#PB_File_SharedRead|#PB_File_NoBuffering)
	readdata(2,ADRDSEQ ,delta)
	REDIM SEQD\TDIF(SEQD\NBSEQ+1)
	readdata(2, Adresdeb,SEQD\NBSEQ+1)
	closefile(2)
endif

deb0.q=ElapsedMilliseconds()
;;;*******************************************************************************************************
  saisie:
mess$+"Décomposition en facteur premier"
If Len(mess$)>99
	MessageRequester("Erreur","entier >0 et <=18446744073709551615 "+#CR$+"Relancer le prg avec <$FFFFFFFFFFFFFFFF") ;
	End
EndIf
; nb$=InputRequester(mess$,"Donnez un entier positif < ou =9223372036854770999","$80123456789ABCDE") ; c'est la limite d'un type q  2^63-1
; nb$=InputRequester(mess$,"Donnez un entier positif < ou =9223372036854770999","$8000000000000001") ; c'est la limite d'un type q  2^63-1
nb$=InputRequester(mess$,"entier >0 et <$FFFFFFFFFFFFFFFF ","$80123456789ABCDE") ; c'est la limite d'un type q  2^63-1
nb$=trim(nb$)	
if mid(nb$,1,1)="$"
  nbhex$=right(nb$,len(nb$)-1)
else 
  nbhex$=dec2hex(nb$)   
endif 
;;;;; ******************* mise en forme du nombre à décomposer et recherche d'une racine carrée****************************************
lnbh=len(nbhex$)
if lnbh>16 OR nbhex$>"FFFFFFFFFFFFFFFE"
 	Goto saisie
EndIf
 
nbhex$=ReplaceString(space(16-lnbh)," ","0")+nbhex$
 A1.q=val("$"+mid(nbhex$,9,16))
 A2.q=val("$"+mid(nbhex$,1,8))
 NB.q=0 
 DIFF.Q=0
 EnableASM  ;; ici recherche en cours pour évaluer le plus précisément la racine carré des nombres >2^63-1
 mov eax,dword[v_A1]
 mov dword[v_nb],eax 
 mov edx,dword[v_A2]
 mov dword [v_nb+4],edx
 sub eax,dword[v_LIMITQ] 
 sbb edx,dword[v_LIMITQ+4]
 mov dword[v_DIFF],eax 
 mov dword [v_DIFF+4],edx
 ;;;;; recherche d'une racine carré de nb à partir de la différence 
 difd.d=diff 
 limitd.d=limitq
 K.d=DIFd/LIMITd ;;;; k est défini sur la plage -1 +1 pour une plage div de +$7FFFFFFFFFFFFFFF à -$7FFFFFFFFFFFFFFF
 RACK.D=SQR(1+K)
 DIVMAXD.d=RACLIM *RACK
  DIVMAX=divmaxd+1  ;;limite de recherche des facteurs premiers
;   printn(_n(nb)+_n(DIVMAX)+_hq(DIVMAX)+_n(divmax)+_d(divmaxd))

;*********************************************************************************************************************************************************
deb1.q=ElapsedMilliseconds()

NBREC.Q=1
DIVIS=2
i=0
fact_prem$=nb$+"="
deb1.q=ElapsedMilliseconds()
;************************************** Recherche des diviseurs premiers qui se trouvent dans p_fact(x) ***********************
iprem=0
DIVIS=nbprem(iprem)
indt=0
While DIVIS<=SEQD\NBPFACT And DIVIS<=DIVMAX 
  enableasm
	MOV ecx,dword[v_DIVIS]  ;;;; la division euclidienne X=kQ+R avec R<Q nous oblige à utiliser une astuce
	XOR EDX,EDX             ;;;; pour éviter avec un diviseur trop petit et un dividente trop grand d'avoir un reste de type qword or EDX est de type dword 
	MOV eax,dword[v_nb+4]
	DIV ecx
	MOV dword[v_quotientp],eax
	MOV eax,dword[v_nb]
	DIV ecx
	CMP edx,0
	JNZ @f
	MOV edx,dword[v_quotientp]
	MOV dword[v_nb+4],edx
	MOV dword[v_nb],eax
	;************************************************ Recherche de la racine carré ***************************************************************
	FILD qword[v_nb]
	FSQRT
	; !FISTTP dword [v_DIVMAX] ; avec arrondi si vous décommentez cette instruction commentez la précédente
	;                           et laissez le ! car l'instruction FISTTP n'est pas reconnue par PureBasic avec !FISTTP c'est bon !!!!
	FISTP qword[v_DIVMAX] ; sans arrondi
		tab(indt)=DIVIS
		indt+1
		DIVMAX=Sqr(nb)
	JMP lab0
	!@@:
    iprem+1
    DIVIS=nbprem(iprem)
    !lab0: 
    disableasm 
Wend
;****************  Partie principale pour le reste de la décomposition **************************************
If DIVIS<=seqd\PREMDIV
	DIVIS=SEQD\PREMDIV
EndIf
ipas=0
;************************************************************************************************************************
pas=SEQD\TDIF(1)
idivis=1
NBSEQ=SEQD\NBSEQ
;   ind=indt*8
Repeat
;   printn(_n(divis)+_n(pas)+_n(idivis))
	;     ;Le nom de @@: signifie étiquette anonyme,vous pouvez avoir défini un grand nombre d'entre eux dans la source.
	;     ;Symbole @b(ou l'équivalent @r)le plus proche des références du label précédent anonyme
	;     ;,symbole @f références plus proche de l'étiquette suivante anonyme.
	;     ;Ces symboles spéciaux sont insensibles à la casse.
	;;*******************************  TESTS EN ASM ********************************************************************
	;*************************************** Division d'un nombre de 64 bits par un nombre de 32 bits ******************************************
	EnableASM
	MOV ecx,dword[v_DIVIS]  ;;;; la division euclidienne X=kQ+R avec R<Q nous oblige à utiliser une astuce
	XOR EDX,EDX             ;;;; pour éviter avec un diviseur trop petit et un dividente trop grand d'avoir un reste de type qword or EDX est de type dword 
	MOV eax,dword[v_nb+4]
	DIV ecx
	MOV dword[v_quotientp],eax
	MOV eax,dword[v_nb]
	DIV ecx
	CMP edx,0
	JNZ @f
	;     MOV dword[v_reste],edx
	; quotientp

	MOV edx,dword[v_quotientp]
	MOV dword[v_nb+4],edx
	MOV dword[v_nb],eax
	;*************************************************************************************************************************************
  ;**** Recherche de la racine carré  pas de pb ici le prg à trouver un diviseur qui n'est pas négatif car (2^64-2)/2 <(2^63-1)*********
	;*************************************************************************************************************************************
FILD qword[v_nb]
	FSQRT
	; !FISTTP dword [v_DIVMAX] ; avec arrondi si vous décommentez cette instruction commentez la précédente
	;                           et laissez le ! car l'instruction FISTTP n'est pas reconnue par PureBasic avec !FISTTP c'est bon !!!!
	FISTP qword[v_DIVMAX] ; sans arrondi
												;                          ;*********************************************** Affectation de la valeur du DIViseur dans la table *******************************************
												;     MOV ecx, [v_DIVIS]; à decommenter si ecx est effacée
												;     mov ebx, [p_Tab]
												;     mov edx, [v_ind]
												;     add ebx,edx
												;     mov[ebx],ecx
												;     add edx,8
												;     mov[v_ind],edx
	DisableASM
	tab(indt)=DIVIS   ;;;; Plus rapide que l'option ASM
	indt+1
	If nb<=1:Break:EndIf
	EnableASM
	;     JMP lab#MacroExpandedCount
	JMP lab
	!@@:
	MOV edx,dword[v_pas]
	ADD dword[v_DIVIS],edx  ; pas de PB l'addition est réalisée en valeur absolue 
													;     donc $FFFFFFFF =>4294967295 donc pas de risque de dépassement sqr(pow(2,63)-1))=>3037000499,9760496922867524030306 <4294967295
	DisableASM
	;********************************************************************************************************************************************************
	;********************************************************************************************************************************************************
	; ;***********************************  Pour accélerer encore l' algo  *******************************
	idivis+1
	If idivis=NBSEQ+1  ;;; inutil de passer en ASM  voir ci dessous c'est plus long de 5 à 6 secondes  sur la recherche de 9223372036854770999
		idivis=1
	EndIf
		pas=SEQD\TDIF(idivis)
EnableASM
	;     mov EAX, dword [v_idivis]
	;     xor edx,edx
	;     mov ECX, dword [v_NBSEQ]
	;     div ECX 
	;     mov  dword [v_idivis],EDX
	;     idivis=idivis%NBSEQ
	;************************************  Pour accélerer encore l' algo vous pouvez décommentez les 7 lignes PB précedentes ******************************
	;*******************************************************************************************************************************************************
	;     !lab#MacroExpandedCount:
	!lab:
	If DIVIS>DIVMAX Or DIVIS<0:Break:EndIf
ForEver


;*************************************************Fin macro decompose2****************************************************************************************************
; indt=ind/8
;If nb >1
	tab(indt)=NB
;EndIf   

;************************************   EDITION ***********************************************
For j=0 To indt
	fact_prem$+tab(j)
	NBREC *tab(j) ; Vérification de la décomposition
	If j<indt
		fact_prem$+"*"
	EndIf
Next
NBRECH$=HEX(NBREC)
NBREC$=HEX2DEC(NBRECH$)
if mid(nb$,1,1)="$"
  NBCONT$="$"+NBRECH$
else 
  NBCONT$=NBREC$
endif  
;essai pour ces nombres 9223372036854769861, 9223372036854770317, 9223372036854770999, 9223372036854773927 et 9223372036854775309 qui ne sont pas premiers.

PrintN(fact_prem$)
PrintN(NBCONT$+": contrôle")
PrintN("temps="+Str(ElapsedMilliseconds()-deb1)+"ms")
PrintN("Recherche des diviseurs <=:"+Str(SEQD\NBPFACT)+" + diviseurs =>"+Str(SEQD\NBPFACT))
PrintN("Nombre d'éléments du vecteur =:"+Str(SEQD\NBSEQ)+"  temps de préparation ="+Str(deb0-deb00))
MessageRequester("Résultat",fact_prem$+#CR$+NBCONT$+": contrôle"+#CR$+"temps="+Str(ElapsedMilliseconds()-deb1)+"ms"+#CR$+"Recherche des diviseurs <=:"+Str(SEQD\NBPFACT)+" + diviseurs =>"+Str(SEQD\NBPFACT)+#CR$+"Nombre d'éléments du vecteur =:"+Str(SEQD\NBSEQ))
Input()
CloseConsole()
; ces nombres 9223372036854769861, 9223372036854770317, 9223372036854770999, 9223372036854773927 et 9223372036854775309 qui ne sont pas premiers; 2295911257* 4017303373=9223372036854769861;==> div=-2147483645 divm=2147483647 DIVMAX=3037000500
; ***************** Ci dessous les Nombres premiers pour les 1000 premiers nombres
; 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 547 557 563 569 571 577 587 593 599
; 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691
; 701 709 719 727 733 739 743 751 757 761 769 773 787 797
; 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887
; 907 911 919 929 937 941 947 953 967 971 977 983 991 997

PS il faut placer ce prg dans le même répertoire que les précédents.

A+
Il est fort peu probable que les mêmes causes ne produisent pas les mêmes effets.(Einstein)
Et en logique positive cela donne.
Il est très fortement probable que les mêmes causes produisent les mêmes effets.
Avatar de l’utilisateur
SPH
Messages : 4721
Inscription : mer. 09/nov./2005 9:53

Re: A propos de nombres premiers

Message par SPH »

Crible d'Eratostene (a la souris) dans une spirale d'Ulam :

Code : Tout sélectionner

Procedure ulam()
  
  Protected nombre = WindowMouseX(0) *2-1; On récupére la coordonnée x de la souris ici
  LoadFont(1, "Arial", 20)  
  StartDrawing(CanvasOutput(0))
  
  x=250
  y=250
  xx=1
  yy=0
  cmb=1
  nb=0
  modulo=0
  cos=0
  la=1
  
  Repeat
    la+2
    x+xx
    y+yy
    nb+1
    If nb=cmb
      modulo+1
      If modulo%2=0
        cmb+1
      EndIf
      
      nb=0
      cos+90
      xx=Cos(Radian(cos))
      yy=Cos(Radian(cos+90))
      
    EndIf
    
    
    i=1
    u=0
        
    If la%nombre=0
      Plot(x, y, RGB(255,255,0)) ; nombre premier
    Else
      Plot(x, y, RGB(50,0,0)) ; nombre non premier
    EndIf
        
  Until cmb=499
    

  DrawingMode( #PB_2DDrawing_Default) ;#PB_2DDrawing_Transparent)
  DrawingFont(FontID(1)) 
  BackColor(0)
  DrawText(20,20,Str(nombre), RGB(255,0,0))

  StopDrawing() 
  
EndProcedure



If OpenWindow(0, 0, 0, 500, 500, "Plot Exemple", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
  
  CanvasGadget(0, 0, 0, 500, 500)
  
  BindGadgetEvent(0, @ulam(), #PB_EventType_MouseMove) ; Quand l'événement de type MouseMove est déclenché, on exécute DrawDesTrucs()
  
  
  Repeat
    Event = WaitWindowEvent()
    
  Until Event = #PB_Event_CloseWindow
  
EndIf
J'ai lu presque toutes les pages de ce theme qui m'interesse.
Je bloque sur l'utilité des ecarts entre les nombres premiers. Quel interet ?
http://HexaScrabble.com/
!i!i!i!i!i!i!i!i!i!
!i!i!i!i!i!i!
!i!i!i!
//// Informations ////
Intel Core i7 4770 64 bits - GTX 650 Ti
Version de PB : 6.00 - 64 bits
PAPIPP
Messages : 534
Inscription : sam. 23/févr./2008 17:58

Re: A propos de nombres premiers

Message par PAPIPP »

Il est fort peu probable que les mêmes causes ne produisent pas les mêmes effets.(Einstein)
Et en logique positive cela donne.
Il est très fortement probable que les mêmes causes produisent les mêmes effets.
Avatar de l’utilisateur
SPH
Messages : 4721
Inscription : mer. 09/nov./2005 9:53

Re: A propos de nombres premiers

Message par SPH »

heuuuu, hors de porté pour moi... 8O

Merci quand meme, j'attendrais que quelqu'un vulgarise tout ca 8)
http://HexaScrabble.com/
!i!i!i!i!i!i!i!i!i!
!i!i!i!i!i!i!
!i!i!i!
//// Informations ////
Intel Core i7 4770 64 bits - GTX 650 Ti
Version de PB : 6.00 - 64 bits
Répondre