PureBasic

Bonjour à tous.

C'est un plaisir d'écrire quelques mots ici.

Je travaille, entre autres choses, sur un projet consistant à éditer des résultats de calculs sur les «écarts entre nombres premiers». Les écarts entre les nombres premiers ont une importance particulière pour certains scientifiques de ma connaissance, et on m'a demandé de l'aide pour concevoir un processus plus rapide pour intégrer des résultats dans leurs bases de données.

Le projet est limité dans l'immédiat à une arithmétique 64 bits, ce qui signifie que nous ne calculons sur des nombres que jusqu'à 1e18, ce qui est suffisant por le moment, pour les résultats finaux que nous voulons obtenir.

De cette façon, j'ai écrit un bout de code pour permettre de calculer les nombres premiers, compter les écarts entre les nombres premiers, et éditer des statistiques sur les quantités des écarts, les résultats normalisés (% vs quantités), et ainsi de suite.

J'ai aussi imaginé un moyen de diviser un gros calcul en tranches, car si on veut par exemple couvrir toute la plage de 0 à 1e18 c'est inaccessible aux humains travaillant avec des PC.

En fait, j'ai limité la première étape de résultats à 1e15 comme borne supérieure de calcul.

Je l'ai testé sur un périmètre 0 à 1e12, par exemple, dans les premières versions de code, et il n'a pas été possible d'obtenir une fin traitement en moins d'un mois. Maintenant 1e12 est traité en une heure sur mon PC fixe.

Le programme utilise le multithreading pour profiter des avantages des processeurs multicœurs.

Pour traiter de très grandes étendues, comme 1e15, il faudrait environ 200 jours. Ceci est la raison principale qui m'a fait choisir la fragmentation du travail en parties plus raisonnables, la reconstitution des résultats consolidés étant faite à la fin de toutes les tranches de calcul.

Le programme dont le code figure ci-dessous contient ce qui est utile pour démarrer, demandant à l'utilisateur d'entrer une plage de calcul (par défaut il propose une plage raisonnable 1-1E9).

Ils y a quelques options supplémentaires, mais c'est un autre débat pas vraiment utile dans la section Trucs et astuces.

Si vous avez un quelconque intérêt sur le sujet et souhaitez plus d'informations, ou si vous voulez des explications plus spécifiques, sauf sur le code qui peut être commenté ici, vous pouvez me contacter en PM. Et si vous pensez que vous avez une bonne tranche de temps disponible sur vos PC à partager avec moi, je peux fournir une partie de traitement dans ce projet (juste à vous rappeler qu'il n'y a ni gloire, ni argent à gagner

).

Voilà, je ne sais pas si ce préambule est bien clair et compréhensible, mais voici le code! Je metrais peut-être à jour, ici, à l'avenir pour les codeurs PureBasic.

Code : Tout sélectionner


;
; Primes Gaps project was started to help some scientists I know to compute primes, and gaps between primes,
; giving the users the possibility to process large ranges, up to 1e12, in the shortest possible time on a single
; PC.
;
; There is no competition with large computers or distributed computing. The first requirement in this project
; is to count primes and gaps between prime numbers with a single usual PC.
;
; A second requirement was to make possible to use multithreading to accelerate results rendering. As most of
; PCs have at least 2 or 4 cores in a single chip, it makes sense to propose a 64 bits, multithreaded processing.
;
; Also the scope of primes to process was intended to reach 1e12 ... in the prime numbers. Two reasons about
; this : the first is that primality tests are heavy enough to make this limit far in rendering counts for such scopes
; in reasonable delay. The other is that many students and researchers use hardware and software with 
; such performances that they do not often walk around larger bounds.
;
; By the way when writing first tests, I could not imagine to solve 1e12 in less than a couple of days on my 4cores
; 2,8 GHz PC. As it is, this program has a 1e14 limitation. But all arithmetics make possible to process easily (but
; not at a glance) up to 1e18.
;
; I will not detail the whole wlak I did, but now this program shows, when using a good PC, that 1e12 full range is
; processed in about one hour and half.
;
; I made this program a bit more detailed in possible usage to let users see the differences in using either regular,
; optimized, with or without multithreading, 32 or 64 bits releases, etc ... this is for fun. But interesting indeed.
;
; I dedicate all this work especially to 2 persons, my father, Gérard Weil and a relative of our family, Gilbert Stork,
; because if I love to write programs, I cannot forget we carry genes of science.  I would like also to thank
; much much much Patrice Brouers for asking me some help on his research, and François Brouers to support
; us in some ways.
;

;
; The sieve of Erathostenes is not a recent discovery. It is also not exactly the way we can imagine in making
; light and fast primality tests when parsing large quantities of integers. I will not explain / describe it here, but
; if you want to check a large integer to be a prime or not, you will probably not start with this method.
;
; Fast explained, you have to dress a sheet with your integers lists in the first column, and fill rows with all
; matches with multiples of 2, 3, 4, 5, 6, 7, 8, 9 ... (simple counters indeed with as many rows as the square root
; of the highest integer listed in the sheet.
;
; Thinking furthermore you will use 2, 3, 5, 7, 11, 13, ... prime numbers only in your sheet header. Well if your
; list of numbers to test is very large, the number of columns will finally be large also.
;
; The number of primes between 0 and n is known to be close to n / (ln(n) - 1).
;
; If you want to check all cells of your sheet (array) having n rows and n / (ln(n) - 1) columns, you will never
; get results probably for centuries. Using n = 1e12 would make necessary to use a memory humankind does not
; have build.
;
; So the only possible way to go on is to split your problem in smaller pieces of work.
;
; The big trick here, to process a large range of integers, is the following :
;
; - split the full range in slices
; - split each slice in pages
; - collect results from pages to slice
; - collect results from slices to range.
; - ... you got it.
;
;
;
;

Structure NOTHING
EndStructure

;
; As the application program is there to render primes' gaps counts, #maxGaps constant
; helps to dim arrays to store gaps counts. 2000 preset is because considering the limit
; of eventually 1e18 (64 bits int numbers) will not generate gaps over 2000 (less in fact).
;
#maxGaps = 2000
;
; 2 constants to identify if the user works with a regular or optimized algorithm.
;
#SliceSoE2 = 1
#SliceSoE2_opt = 2

Structure SLICE ; SLICE structure formats a record containing all needed information for a slice.
  nmin.q                                        ; nmin / nmax are lower and upper bounds of the slice
  nmax.q                                       ;
  n.q                                            ; 
  nPrimes.q                                    ; the number of primes found in the slice
  FirstPrime.q                                 ; the first prime integer found in the slice
  LastPrime.q                                 ; ... and the last
  Gaps.q[#maxGaps]                        ; a static array to collcet gaps' counts
  FirstPrimeGaps.q[#maxGaps]           ; a static array remembering the first prime found having a given gap
  StartTime.q                                 ; statistic information about slice processing duration
  EndTime.q                                   ;
  ID.q
  Status.q
EndStructure

Structure Q2                                 ; small structure to facilitate final checks of results
  p.q
  q.q
  check.s
EndStructure

Define.NOTHING

;
; A set of paragraphs here are placed to understand and report about running PC characteristics.
; This is of no use except for interested users.
;
  Select #PB_Compiler_OS
    Case #PB_OS_Windows
      OSType.s = "Windows"
    Case #PB_OS_Linux
      OSType.s = "Linux"
    Case #PB_OS_MacOS
      OSType.s = "Mac OS X"
    Default
      OSType.s = "Unknown OS type"
  EndSelect
  Select #PB_Compiler_Processor
    Case #PB_Processor_x86
      ProcessorType.s = "x86"
    Case #PB_Processor_x64
      ProcessorType.s = "x64"
    Default
      ProcessorType.s = "Unknown processor type"
  EndSelect
  
  
  Select OSVersion()
    Case #PB_OS_Windows_NT3_51
      sOSVersion.s = "Windows NT 3.51"
    Case #PB_OS_Windows_95
      sOSVersion.s = "Windows 95"
    Case #PB_OS_Windows_NT_4
      sOSVersion.s = "Windows NT 4"
    Case #PB_OS_Windows_98
      sOSVersion.s = "Windows 98"
    Case #PB_OS_Windows_ME
      sOSVersion.s = "Windows ME"
    Case #PB_OS_Windows_2000
      sOSVersion.s = "Windows 2000"
    Case #PB_OS_Windows_XP
      sOSVersion.s = "Windows XP"
    Case #PB_OS_Windows_Server_2003
      sOSVersion.s = "Windows Server 2003"
    Case #PB_OS_Windows_Vista
      sOSVersion.s = "Windows Vista"
    Case #PB_OS_Windows_Server_2008
      sOSVersion.s = "Windows Server 2008"
    Case #PB_OS_Windows_7
      sOSVersion.s = "Windows 7"
    Case #PB_OS_Windows_Server_2008_R2
      sOSVersion.s = "Windows Server 2008 R2"
    Case #PB_OS_Windows_8
      sOSVersion.s = "Windows 8"
    Case #PB_OS_Windows_Server_2012
      sOSVersion.s = "Windows Server 2012"
    Case #PB_OS_Windows_8_1
      sOSVersion.s = "Windows 8.1"
    Case #PB_OS_Windows_Server_2012_R2
      sOSVersion.s = "Windows Server 2012 R2"
    Case #PB_OS_Windows_10
      sOSVersion.s = "Windows 10"
    Case #PB_OS_Windows_Future
      sOSVersion.s = "New Windows version (not existing when the program was written)"
    Case #PB_OS_Linux_2_2
      sOSVersion.s = "Linux 2.2"
    Case #PB_OS_Linux_2_4
      sOSVersion.s = "Linux 2.4"
    Case #PB_OS_Linux_2_6
      sOSVersion.s = "Linux 2.6"
    Case #PB_OS_Linux_Future
      sOSVersion.s = "New Linux version (not existing when the program was written)"
    Case #PB_OS_MacOSX_10_0
      sOSVersion.s = "MacOSX 10.0"
    Case #PB_OS_MacOSX_10_1
      sOSVersion.s = "MacOSX 10.1"
    Case #PB_OS_MacOSX_10_2
      sOSVersion.s = "MacOSX 10.2"
    Case #PB_OS_MacOSX_10_3
      sOSVersion.s = "MacOSX 10.3"
    Case #PB_OS_MacOSX_10_4
      sOSVersion.s = "MacOSX 10.4"
    Case #PB_OS_MacOSX_10_5
      sOSVersion.s = "MacOSX 10.5"
    Case #PB_OS_MacOSX_10_6
      sOSVersion.s = "MacOSX 10.6"
    Case #PB_OS_MacOSX_10_7
      sOSVersion.s = "MacOSX 10.7"
    Case #PB_OS_MacOSX_10_8
      sOSVersion.s = "MacOSX 10.8"
    Case #PB_OS_MacOSX_10_9
      sOSVersion.s = "MacOSX 10.9"
    Case #PB_OS_MacOSX_10_10
      sOSVersion.s = "MacOSX 10.10"
    Case #PB_OS_MacOSX_10_11
      sOSVersion.s = "MacOSX 10.11"
    Case #PB_OS_MacOSX_Future
      sOSVersion.s = "New MacOS X version (not existing when the program was written)"
    Default
      sOSVersion.s = "Unknown OS version return code"
  EndSelect
  
  PB_Compiler_Version.s = Str(#PB_Compiler_Version)
  PB_Version.s = Mid(PB_Compiler_Version, 1, 1) + "." + Mid(PB_Compiler_Version, 2)
  EnvHeader.s = #PB_Compiler_Filename + " generated " + FormatDate("%yyyy/%mm/%dd %hh:%ii:%ss", #PB_Compiler_Date) + " with PureBasic " + PB_Version + #CRLF$ +
                "Available for " + OSType + " on " + ProcessorType + " architectures" + #CRLF$ + 
                "running on : " + UserName() + "@" + ComputerName() + #CRLF$ +
                "CPU name : " + CPUName() + #TAB$ + "Processors count : " + Str(CountCPUs(#PB_System_CPUs)) + #TAB$ + "Avail. processors : " + Str(CountCPUs(#PB_System_ProcessCPUs)) + #CRLF$ +
                "Total phys. mem : " + FormatNumber(MemoryStatus(#PB_System_TotalPhysical), 0, ",", ".") + #TAB$ + "Total free. mem (@prog start) : " + FormatNumber(MemoryStatus(#PB_System_FreePhysical), 0, ",", ".") + #CRLF$ + 
                "OS version on running device : " + sOSVersion
  
;
; lower / upper bounds of range to test
;
  nmin.q = 1
  nmax.q = 1E9
  
;
; Set some initial values and get user customized parameters
;
  ProgramParameters.s
  CountProgramParameters.q = CountProgramParameters()
  If CountProgramParameters
      i.q = 1
      Repeat
        ProgramParameters + " " + ProgramParameter()
        i + 1
      Until i > CountProgramParameters
    Else
      DirectoryNumber.q = ExamineDirectory(#PB_Any, GetCurrentDirectory(), "*.*")
      If DirectoryNumber
        While NextDirectoryEntry(DirectoryNumber)
          If DirectoryEntryAttributes(DirectoryNumber) & #PB_FileSystem_Hidden = 0
              ProgramParameters_FileName.s = DirectoryEntryName(DirectoryNumber)
              FileNumber.q = ReadFile(#PB_Any, ProgramParameters_FileName)
              If FileNumber
                  Tag.s = "Primes Gaps automation "
                  lTag.q = Len(Tag)
                  If Left(ProgramParameters_FileName, lTag) = Tag
                      ProgramParameters = ReadString(FileNumber, #PB_Ascii)
                      CloseFile(FileNumber)
                      SetFileAttributes(ProgramParameters_FileName, #PB_FileSystem_Hidden)
                      Break
                  EndIf
                  CloseFile(FileNumber)
              EndIf
          EndIf
        Wend
        FinishDirectory(DirectoryNumber)
      EndIf
      
  EndIf
  
  If ProgramParameters <> ""
      q.s = ProgramParameters
    
    Else
      q.s = InputRequester("upper bound value, or h for quick help", "Enter the min-max values of range", Str(nmin) + "-" + Str(nmax))
      If q = "?" Or LCase(q) = "h" Or LCase(q) = "help"
        MessageRequester("Some help ?", "Just give a min-max value to parse for primes and primes' gaps" + #LF$ +
                                        "If the program is launched from console, arguments are expected in the command line." + #LF$ + #LF$ +
                                        "You can provide some more parameters after min-max values, using /name=value :" + #LF$ +
                                        "If possible argument value is a choice, defaults are displayed here first, and alternatives in second." + #LF$ +#LF$ +
                                        "file=[yes/no] to generate a file with results." + #LF$ +
                                        "file=open  optionally you can ask the file to open at program end." + #LF$ +
                                        "algo=[opt/regular] to specify which algorithm you prefer." + #LF$ +
                                        "thread=[yes/no] to enable / disable multithreading." + #LF$ +
                                        "wait=[yes/no] wait / don't wait for user close console at end.", #PB_MessageRequester_Ok)
          End
      EndIf
  EndIf
  
  If q <> ""
      If Trim(LCase(q)) = "quit"
          End
      EndIf
      If FindString(q, "-", 1)
          LHS.s = Trim(StringField(q, 1, "-"))
          RHS.s = Trim(StringField(q, 2, "-"))
          If FindString(RHS, " ", 1)
              RHS = StringField(RHS, 1, " ")
          EndIf
          NewList Args.s()
          iField.q = 2
          While StringField(q, iField, "/") <> ""
            AddElement(Args())
            Args() = StringField(q, iField, "/")
            iField + 1
          Wend
          nmin = Val(LHS)
          nmax = Val(RHS)
        Else
          nmax = Val(q)
      EndIf
  EndIf
  
  maxStoredPrimes.q = Sqr(nmax)
  ThreadedFlag.q = #True
  Algorithm.q = #SliceSoE2_opt
  WriteDocument.q = #True
  OpenDocument.q = #False
  Wait_At_End.q = #True
  
  ForEach Args()
    LHS.s = LCase(Trim(StringField(Args(), 1, "=")))
    RHS.s = LCase(Trim(StringField(Args(), 2, "=")))
    Select LHS
      Case "file"
        If RHS = "yes"
          WriteDocument.q = #True
        EndIf
        If RHS = "no"
          WriteDocument.q = #False
        EndIf
        If RHS = "open"
          OpenDocument.q = #True
        EndIf
      Case "algo"
        If RHS = "regular"
          Algorithm.q = #SliceSoE2
        EndIf
        If RHS = "opt"
          Algorithm.q = #SliceSoE2_opt
        EndIf
      Case "thread"
        If RHS = "no"
          ThreadedFlag.q = #False
        EndIf
        If RHS = "yes"
          ThreadedFlag.q = #True
        EndIf
      Case "wait"
        If RHS = "no"
          Wait_At_End.q = #False
        EndIf
        If RHS = "yes"
          Wait_At_End.q = #True
        EndIf
      Default
    EndSelect
  Next
  
  If Algorithm = #SliceSoE2_opt
      sAlgorithm.s = "#SliceSoE2_opt"
    ElseIf Algorithm = #SliceSoE2
      sAlgorithm.s = "#SliceSoE2"
  EndIf
  
  If Not WriteDocument
    OpenDocument = #False
  EndIf
  
;
; Inform optimistics that tomorrow is another day
;
  If nmax > 1E15
      MessageRequester("Sorry !", "This app. can't process over " + FormatNumber(1E15, 0, ",", ".") + " upper bound", #PB_MessageRequester_Ok)
      End
  EndIf
  
;
; Some variables and memory space to run process
;
; Threads() is a linked list that will be used to get and manage threads
; 
Global NewList Threads.q()
Global DoublePrimeTestCheck.q = #False
Global PageSizeLimit.q = 500000 ; The page size limit (a page is a part of a slice here).
    nSlicesMaxLimit.q = 1024 * 32 ; The limit number of slices in which the range is splitted.
    ThDelay1.q = 0 ; In case of multithreading use, ThDelay1/2 allows to adapt theory to reality.
    ThDelay2.q = 0 ; These delays should not be necessary, but in case of frequent crash, just place
                         ; a short delay like 1 or 5. Here are milliseconds.
    maxConcurrentThreads.q = 16 ; This parameter indicates how many concurrent threads should
                                             ; be sent to CPU at the same moment. The CPU will process  only
                                             ; only a thread by core / CPU at a time but each core / CPU can
                                             ; manage a queue. Depending on hardware, this parameter may
                                             ; help to tune performances.

PI_Primes.q = 1.2 * maxStoredPrimes / (Log(maxStoredPrimes) - 1)
Global Dim Primes.q(PI_Primes)
Global nStoredPrimes.q
Global nSlicesStarted.q
Global nSlicesDone.q

    If nmax < 1E6
        nSlices.q = 1
      ElseIf nmax < 1E7
        nslices = 4
      ElseIf nmax < 1E8
        nslices = 8
      ElseIf nmax < 1E9
        nslices = 16
      ElseIf nmax < 1E10
        nslices = 32
      ElseIf nmax < 1E11
        nslices = 64
      ElseIf nmax < 1E12
        nslices = 128
      ElseIf nmax < 1E13
        nslices = 256
      ElseIf nmax < 1E14
        nslices = 512
      Else
        nslices = 1024
    EndIf
    
;     nSlices = 256
;     maxConcurrentThreads = 256
;     PageSizeLimit.q = 0
;     
    
;
; A small piece of code to get reference data. This to check results if possible.
; PI_n_ref will collect PI(n) for a limited set of values and also the corresponding gaps counts.
;
NewList PI_n_ref.Q2()
Repeat
  Read.q range.q
  Read.q value.q
  If range = 999999999 And value = 999999999
      Break
  EndIf
  AddElement(PI_n_ref())
  PI_n_ref()\p = range
  PI_n_ref()\q = value
ForEver
ForEach PI_n_ref()
  Read.q i.q
  If i = 999999999
      Break
  EndIf
  If i = PI_n_ref()\p
      Read.s PI_n_ref()\check
  EndIf
Next

;
; Here is the trick ! A sieve of Erathostenes algorithm to parse ranges of integers with counters.
; The most important trick is to parse any range with a n upper bound using sqr(n) counters.
; Also each counters must be initialized to the appropriate value starting with lower bound.
;
; To perform this as fast as possible, the program uses a stored primes table, containing a set
; of precalculated primes up to the value of sqr(max), where max is the highest possible upper bound.
;
; That's it and nothing more. But the way to split a large range in slices, and each slice in pages makes
; the difference.
;
; SliceSoE2 receives a data record, SLICE structure, containing bounds of slice, and some space for
; results.
; The record will be fullfilled with primes count, gaps counts, first and last primes of slice.
;
; Collecting all the slices, the main program code will process a global merge of all data.
;
; The following procedure is the regular code one, available for both x86 and x64 CPUs, and processing
; 64 bits arithmetic. It means that performances are poor on x86.
;
Procedure SliceSoE2(*Slice.SLICE)
  Shared ID.q
  nSlicesStarted + 1
  *Slice\StartTime = ElapsedMilliseconds()
  *Slice\ID = ID
  *Slice\Status = 1
  StartDate.q = Date()
  ID + 1
  Gnmin.q = *Slice\nmin
  Gnmax.q = *Slice\nmax
  nPrimes.q = 0
  PageSize.q = Sqr(Gnmax)
  If PageSizeLimit
      If PageSize > PageSizeLimit
          PageSize = PageSizeLimit
      EndIf
  EndIf
  FirstPrime.q = 0
  LastPrime.q = 0
  nSoE.q = 0
  Gn.q = Gnmin
  Repeat
    nmin.q = Gn
    nmax.q = Gn + PageSize - 1
    If nmax > Gnmax
        nmax = Gnmax
    EndIf
    maxSoE.q = Sqr(nmax)
    While nSoE <= nStoredPrimes
      If Primes(nSoE) > maxSoE
          Break
      EndIf
      nSoE + 1
    Wend
    nSoE - 1
    If nmin <= 1
        nmin = 1
        nPrimes + 1
        FirstPrime = 2
        LastPrime = 2
    EndIf
    If nmin & 1 = 0
        nmin + 1
    EndIf
    iSoE.q = 0
    If nmax & 1 = 0
        nmax - 1
    EndIf
    ThisPageSize.q = nmax - nmin
    If ThisPageSize > 0
        Dim Page.q(ThisPageSize)
        If nmin = 1
            Page(0) = 1
        EndIf
        iSoE.q = nSoE
        Repeat
          thisprime.q = Primes(iSoE)
          If nmin < thisprime
              ptr.q = thisprime << 1 - nmin
            Else
              r.q = nmin % thisprime
              If r = 0
                  ptr.q = 0
                Else
                  ptr.q = thisprime - r
              EndIf
          EndIf
          While ptr <= ThisPageSize
            If ptr & 1 = 0
                Page(ptr) = 1
            EndIf
            ptr + thisprime
          Wend
          iSoE - 1
        Until iSoE < 0
        i.q = 0
        Repeat
          If Page(i) = 0
              nPrimes + 1
              If FirstPrime = 0
                  FirstPrime = i + nmin
                  LastPrime = FirstPrime
                Else
                  Gap.q = i + nmin - LastPrime
                  If *Slice\Gaps[Gap] = 0
                      *Slice\FirstPrimeGaps[Gap] = LastPrime
                  EndIf
                  *Slice\Gaps[Gap] + 1
                  LastPrime = i + nmin
              EndIf
          EndIf
          i + 2
        Until i > ThisPageSize ; Next
    EndIf
    Gn + PageSize
  Until Gn > Gnmax
  *Slice\FirstPrime = FirstPrime
  *Slice\LastPrime = LastPrime
  *Slice\nPrimes + nPrimes
  *Slice\EndTime = ElapsedMilliseconds()
  *Slice\Status = 2
  nSlicesDone + 1
EndProcedure

;
; Optimized code is proposed in a deep rework for x64 and lighter one for x86 because the whole code
; is definitely a 64 bits one.
;
; In the 64 bits part, FASM is used to make all scopes of code as fast as possible using all registers to
; have calculations as much as possible inside the core (register to register processing eveywhere it is possible).
;
; Assuming this, the optimized code does exactly the same compared to the regular part.
;
Procedure SliceSoE2_opt(*Slice.SLICE)
  Shared ID.q
  nSlicesStarted + 1
  *Slice\StartTime = ElapsedMilliseconds()
  *Slice\ID = ID
  *Slice\Status = 1
  StartDate.q = Date()
  ID + 1
  Gnmin.q = *Slice\nmin
  Gnmax.q = *Slice\nmax
  nPrimes.q = 0
  PageSize.q = Sqr(Gnmax)
  If PageSizeLimit
      If PageSize > PageSizeLimit
          PageSize = PageSizeLimit
      EndIf
  EndIf
  FirstPrime.q = 0
  LastPrime.q = 0
  nSoE.q = 0
  Gn.q = Gnmin
  Repeat
    nmin.q = Gn
    nmax.q = Gn + PageSize - 1
    If nmax > Gnmax
        nmax = Gnmax
    EndIf
    maxSoE.q = Sqr(nmax)
    
    !  MOV     r11, qword [p.v_nSoE]
    !  MOV     r12, r11
    !  SAL      r12, 3
    !  ADD      r12, qword [a_Primes]
    !  MOV     r13, qword [p.v_maxSoE]
! SliceSoE2_opt_while_nSoE_le_nStoredPrimes_while:
    !  CMP     r11, qword [v_nStoredPrimes] ; While nSoE <= nStoredPrimes
    !  JG        SliceSoE2_opt_while_nSoE_le_nStoredPrimes_endloop
      !  CMP     qword [r12], r13 ; If Primes(nSoE) > maxSoE
      !  JLE       @f
          !  JMP      SliceSoE2_opt_while_nSoE_le_nStoredPrimes_endloop ; Break
! @@: ; EndIf
      !  INC      r11 ; nSoE + 1
      !  ADD     r12, 8
    !  JMP     SliceSoE2_opt_while_nSoE_le_nStoredPrimes_while
! SliceSoE2_opt_while_nSoE_le_nStoredPrimes_endloop:
    !  DEC      r11 ; nSoE - 1
    !  MOV     qword [p.v_nSoE], r11
    
    If nmin <= 1
        nmin = 1
        nPrimes + 1
        FirstPrime = 2
        LastPrime = 2
    EndIf
    If nmin & 1 = 0
        nmin + 1
    EndIf
    iSoE.q = 0
    If nmax & 1 = 0
        nmax - 1
    EndIf
    ThisPageSize.q = nmax - nmin
    
    *pGaps = @*Slice\Gaps[0]
    *pFirstPrimeGaps = @*Slice\FirstPrimeGaps[0]
    
    If ThisPageSize > 0
        Dim Page.q(ThisPageSize)
        
; ptr => r8   thisprime => r9   iSoE => r10   nmin => r11   ThisPageSize => r12   @Primes() => r13    @Page() => r14
        !  MOV     r11, qword [p.v_nmin]
        !  MOV     r14, qword [p.a_Page]
        !  CMP     r11, 1 ; If nmin = 1
        !  JNE      @f
            !  MOV     qword [r14], 1 ; Page(0) = 1
! @@: ; EndIf
        !  MOV     r10, qword [p.v_nSoE] ; iSoE.q = nSoE
        !  MOV     r12, qword [p.v_ThisPageSize]
        !  MOV     r13, qword [a_Primes]
! SliceSoE2_opt_iSoE_nSoE_0_loop: ; Repeat
          !  MOV     rax, r10 ; thisprime.q = Primes(iSoE)
          !  SAL      rax, 3
          !  MOV     r9, qword [r13+rax]
          !  MOV     rax, r11 ; If nmin < thisprime
          !  CMP     rax, r9
          !  JGE      @f
              !  MOV     rax, r9 ; ptr.q = thisprime << 1 - nmin
              !  SAL      rax, 1
              !  SUB      rax, r11
              !  MOV     r8, rax
              !  JMP      SliceSoE2_opt_while_ptr_le_ThisPageSize_loop
! @@: ; Else
              !  MOV     rax, r11 ; r.q = nmin % thisprime
              !  MOV     rcx, r9
              !  CQO
              !  IDIV     rcx
              !  CMP     rdx, 0 ; If r = 0
              !  JNZ      @f
                  !  MOV     r8, 0 ; ptr.q = 0
                  !  JMP     SliceSoE2_opt_while_ptr_le_ThisPageSize_loop
! @@: ; Else
                  !  MOV     r8, r9 ; ptr.q = thisprime - r
                  !  SUB      r8, rdx
! SliceSoE2_opt_while_ptr_le_ThisPageSize_loop: ; While ptr <= ThisPageSize
          !  CMP     r8, r12
          !  JG        SliceSoE2_opt_while_ptr_le_ThisPageSize_endloop
            !  MOV     rax, r8 ; If ptr & 1 = 0
            !  AND      rax, 1
            !  JNZ       @f
                !  MOV     rax, r8 ; *Page = *pPage + ptr << 3
                !  SAL      rax, 3
                !  INC      qword [r14+rax] ; *Page\q = 1 ; Page(ptr) = 1
! @@: ; EndIf
            !  ADD     r8, r9 ; ptr + thisprime
          !  JMP      SliceSoE2_opt_while_ptr_le_ThisPageSize_loop ; Wend
! SliceSoE2_opt_while_ptr_le_ThisPageSize_endloop:
          !  DEC     r10 ; iSoE - 1
        !  CMP     r10, 0 ; Until iSoE < 0
        !  JGE      SliceSoE2_opt_iSoE_nSoE_0_loop
! SliceSoE2_opt_iSoE_nSoE_0_endloop:

; i => r8   Gap => r9   nmin => r11   ThisPageSize => r12   @*Slice\Gaps[0] => r13   @Page() => r14   @*Slice\FirstPrimeGaps[0] => r15
        !  MOV     r8, 0
        !  MOV     r13, qword [p.p_pGaps]
        !  MOV     r15, qword [p.p_pFirstPrimeGaps]
! SliceSoE2_repeat_i_0_ThisPageSize_loop: ; Repeat
          !  CMP      qword [r14], 0 ; If *Page\q = 0 ; Page(i) = 0
          !  JNZ       SliceSoE2_opt_Page_i_eq_0_endtest
              !  INC      qword [p.v_nPrimes] ; nPrimes + 1
              !  CMP     qword [p.v_FirstPrime], 0 ; If FirstPrime = 0
              !  JNZ      @f
                  !  MOV     rax, r8 ; FirstPrime = i + nmin
                  !  ADD      rax, r11 ; qword [p.v_nmin]
                  !  MOV     qword [p.v_FirstPrime], rax
                  !  MOV     qword [p.v_LastPrime], rax ; LastPrime = FirstPrime
                  !  JMP      SliceSoE2_opt_FirstPrime_eq_0_endtest
! @@: ; Else
                  !  MOV     r9, r8 ; Gap.q = i + nmin - LastPrime
                  !  ADD      r9, r11 ; qword [p.v_nmin]
                  !  SUB      r9, qword [p.v_LastPrime]
                  !  SAL      r9, 3 ; *Gaps = *pGaps + Gap << 3
                  !  CMP     qword [r13+r9], 0 ; If *Gaps\q = 0 ; *Slice\Gaps[Gap] = 0
                  !  JNZ      @f
                      !  MOV     rbx, qword [p.v_LastPrime]  ; *FirstPrimeGaps = *pFirstPrimeGaps + Gap << 3 ; *FirstPrimeGaps\q = LastPrime ; *Slice\FirstPrimeGaps[Gap] = LastPrime
                      !  MOV     qword [r15+r9], rbx
! @@: ; EndIf
                  !  INC      qword [r13+r9] ; qword [p.p_Gaps] ; *Gaps\q + 1 ; *Slice\Gaps[Gap] + 1
                  !  MOV     rax, r8 ; LastPrime = i + nmin
                  !  ADD      rax, r11 ; qword [p.v_nmin]
                  !  MOV     qword [p.v_LastPrime], rax
! SliceSoE2_opt_FirstPrime_eq_0_endtest: ; EndIf
! SliceSoE2_opt_Page_i_eq_0_endtest: ; EndIf
          !  ADD     r8, 2 ; i + 2
          !  ADD     r14, 16 ; *Page + 16
        !  CMP      r8, r12 ; qword [p.v_ThisPageSize] ; Until i > ThisPageSize ; Next
        !  JLE       SliceSoE2_repeat_i_0_ThisPageSize_loop
    EndIf
    Gn + PageSize
  Until Gn > Gnmax
  *Slice\FirstPrime = FirstPrime
  *Slice\LastPrime = LastPrime
  *Slice\nPrimes + nPrimes
  *Slice\EndTime = ElapsedMilliseconds()
  *Slice\Status = 2
  nSlicesDone + 1
EndProcedure

;
; Main starts here
;

SetPriorityClass_(GetCurrentProcess_(), #BELOW_NORMAL_PRIORITY_CLASS)
If OpenConsole()
    Text.s = ""
    Env.s = ProgramFilename() + #CRLF$ +
            EnvHeader + #CRLF$ + #CRLF$ + 
            "nmin = " + FormatNumber(nmin, 0, ",", ".") + #CRLF$ +
            "nmax = " + FormatNumber(nmax, 0, ",", ".") + #CRLF$ +
            "maxStoredPrimes = " + FormatNumber(maxStoredPrimes, 0, ",", ".") + #CRLF$ +
            "ThreadedFlag = " + FormatNumber(ThreadedFlag, 0, ",", ".") + #CRLF$ +
            "Algorithm = " + sAlgorithm + #CRLF$ +
            "WriteDocument = " + FormatNumber(WriteDocument, 0, ",", ".") + #CRLF$ +
            "OpenDocument = " + FormatNumber(OpenDocument, 0, ",", ".") + #CRLF$ +
            "PageSizeLimit = " + FormatNumber(PageSizeLimit, 0, ",", ".") + #CRLF$ +
            "nSlicesMaxLimit = " + FormatNumber(nSlicesMaxLimit, 0, ",", ".") + #CRLF$ +
            "maxConcurrentThreads = " + FormatNumber(maxConcurrentThreads, 0, ",", ".")
    
    PrintN(Env)
    Text + Env + #CRLF$
    ConsoleTitle(FormatNumber(nmin, 0, ",", ".") + " - " + FormatNumber(nmax, 0, ",", "."))
    StartDate.q = Date()
    tz.q = ElapsedMilliseconds()
    nStoredPrimes.q = 0
    n.q = 0
    Repeat
      Select n
        Case 0, 1
          n + 1
        Case 2
          Primes(nStoredPrimes) = n
          nStoredPrimes + 1
          n + 1
        Case 3, 5, 7
          Primes(nStoredPrimes) = n
          nStoredPrimes + 1
          n + 2
        Default
          If n & 1 = 0
              n + 1
          EndIf
          If n > 5
              b.d = 1 / 3
              c.q = n * b
              If n - 3 * c = 0
                  IsPrime.q = #False
                  Goto IsPrime_test_point
              EndIf
              dmax.q = Sqr(n)
              d.q = 5
              While d <= dmax
                b.d = 1 / d
                c.q = n * b
                If n - d * c = 0
                    IsPrime.q = #False
                    Goto IsPrime_test_point
                EndIf
                d + 2
                b.d = 1 / d
                c.q = n * b
                If n - d * c = 0
                    IsPrime.q = #False
                    Goto IsPrime_test_point
                EndIf
                d + 4
              Wend
              IsPrime.q = #True
              Goto IsPrime_test_point
            Else
              Select n
                Case 3, 5
                  IsPrime.q = #True
                  Goto IsPrime_test_point
                Default
                  IsPrime.q = #False
                  Goto IsPrime_test_point
              EndSelect
          EndIf
IsPrime_test_point:          
          If IsPrime
            Primes(nStoredPrimes) = n
            nStoredPrimes + 1
          EndIf
          n + 2
      EndSelect
      If n & $FFFFF = $FFFFF
          PrintN(StrD(100 * n / maxStoredPrimes, 3) + " %" + #TAB$ + Str(nStoredPrimes))
      EndIf
      
    Until n > maxStoredPrimes
    nStoredPrimes - 1
    
    Time.d = (ElapsedMilliseconds() - tz) / 1000
    Line.s = "Stored primes array fetched with : " + FormatNumber(nStoredPrimes + 1, 0, ",", ".") + " primes in " + StrD(Time, 3) + " s"
    PrintN(Line)
    Text + Line + #CRLF$
    
    SliceSize.q = (nmax - nmin) / nSlices
    
    nSlices - 1
    Dim Slice2.SLICE(nSlices)
    Dim FirstPrimes.q(nSlices)
    Dim LastPrimes.q(nSlices)
    Dim Gaps.q(#maxGaps - 1)
    Dim GapsMod3.q(2, #maxGaps - 1)
    Dim FirstPrimeGaps.q(#maxGaps - 1)
    
    Line.s = "Slice2SoE - starting range" + #TAB$ + FormatNumber(nmin, 0, ",", ".") + " - " + FormatNumber(nmax, 0, ",", ".") + #TAB$ + "Using " + Str(nSlices + 1) + " threads" + #TAB$ + FormatDate("%yyyy/%mm/%dd %hh:%ii:%ss", Date())
    PrintN(Line)
    Text + Line + #CRLF$
    
    tz.q = ElapsedMilliseconds()
    
    PrintN("Effective slice size : " + FormatNumber(SliceSize, 0, ",", "."))

    nPrimes.q = 0
    nGaps.q = 0
    thisnmin.q = nmin
    thisnmax.q = nmin + SliceSize - 1
    
    kSlicesDone.d = 0
    iSlice.q = 0
    Repeat
      Slice2(iSlice)\nmin = thisnmin
      Slice2(iSlice)\nmax = thisnmax
      Slice2(iSlice)\Status = 0
      Slice2(iSlice)\ID = iSlice
      If ThreadedFlag
          AddElement(Threads())
          Select Algorithm
            Case #SliceSoE2
              Threads() = CreateThread(@SliceSoE2(), @Slice2(iSlice))
            Case #SliceSoE2_opt
              Threads() = CreateThread(@SliceSoE2_opt(), @Slice2(iSlice))
          EndSelect
        Else
          Select Algorithm
            Case #SliceSoE2
              SliceSoE2(Slice2(iSlice))
            Case #SliceSoE2_opt
              SliceSoE2_opt(Slice2(iSlice))
          EndSelect
      EndIf
      thisnmin = thisnmax + 1
      thisnmax = thisnmin + SliceSize - 1
      If iSlice = nSlices - 1
        thisnmax = nmax
      EndIf
      If ThreadedFlag
          Repeat
            nThreads.q = nSlicesStarted - nSlicesDone
            Delay(ThDelay1)
          Until nThreads < maxConcurrentThreads
      EndIf
      kSlicesDone.d = nSlicesDone / nSlices
      PastTime.d = (ElapsedMilliseconds() - tz) / 1000
      If ExpectedTime.d = #INFINITE Or PastTime = #INFINITE
          PrintN("Some cycles required to evaluate statistics")
        Else
          PrintN(StrD(kSlicesDone * 100, 3) + " %" + #TAB$ +
                 "Elaps. " + StrD(PastTime, 1) + " s    " +
                 "Est. " + sExPectedDate.s + " (" + StrD(ExpectedTime.d, 1) + " s " +
                 "/ still " + StrD(ExpectedTime - PastTime.d, 1) + " s)")
      EndIf
      ExpectedTime.d = PastTime * Pow(31, Log10(nmax - nmin) - Log10(nSlicesDone * SliceSize))
      sExPectedDate.s = FormatDate("%yyyy/%mm/%dd %hh:%ii:%ss", StartDate + ExpectedTime)
      iSlice + 1
    Until iSlice > nSlices
    
    If ThreadedFlag
        Repeat
          Ready.q = #True
          ForEach Threads()
            If IsThread(Threads())
              Ready = #False
            EndIf
          Next
          Delay(ThDelay2)
        Until Ready
    EndIf
    
    iSlice.q = 0
    Repeat
      FirstPrimes(iSlice) = Slice2(iSlice)\FirstPrime
      LastPrimes(iSlice) = Slice2(iSlice)\LastPrime
      iSlice + 1
    Until iSlice > nSlices
    
    iSlice.q = 0
    Repeat
      If iSlice > 0
        Gap.q = FirstPrimes(iSlice) - LastPrimes(iSlice - 1)
        Gaps(Gap) + 1
        If FirstPrimeGaps(Gap) > LastPrimes(iSlice - 1)
          FirstPrimeGaps(Gap) = LastPrimes(iSlice - 1)
        EndIf
      EndIf
      nPrimes + Slice2(iSlice)\nPrimes
      iGap.q = 0
      Repeat
        If Slice2(iSlice)\Gaps[iGap] <> 0
            nGaps + Slice2(iSlice)\Gaps[iGap]
            Gaps(iGap) + Slice2(iSlice)\Gaps[iGap]
            If FirstPrimeGaps(iGap) > Slice2(iSlice)\FirstPrimeGaps[iGap] Or FirstPrimeGaps(iGap) = 0
                FirstPrimeGaps(iGap) = Slice2(iSlice)\FirstPrimeGaps[iGap]
            EndIf
        EndIf
        iGap + 1
      Until iGap > #maxGaps - 1
      iSlice + 1
    Until iSlice > nSlices
    
    PrintN("")
    Time.d = (ElapsedMilliseconds() - tz) / 1000
    Line = "Slice2SoE - parsed range" + #TAB$ + FormatNumber(Slice2(0)\nmin, 0, ",", ".") + " - " + FormatNumber(Slice2(nSlices)\nmax, 0, ",", ".") + #TAB$ + "Done in " + StrD(Time, 3) + " s" + #TAB$ + FormatDate("%yyyy/%mm/%dd %hh:%ii:%ss", Date())
    PrintN(Line)
    Text + Line + #CRLF$
    Line = "Slice2SoE - found " + FormatNumber(nPrimes, 0, ",", ".") + " primes"
    PrintN(Line)
    Text + Line + #CRLF$
    PrintN("")
    Text + #CRLF$
    
    Line = "First prime in range : " + FormatNumber(Slice2(0)\FirstPrime, 0, ",", ".") + #CRLF$ + "Last prime in range : " + FormatNumber(Slice2(nSlices)\LastPrime, 0, ",", ".")
    PrintN(Line)
    Text + Line + #CRLF$
    
    nGaps.q = 0
    iGap.q = 0
    Repeat
      If Gaps(iGap) <> 0
          nGaps + Gaps(iGap)
        EndIf
        iGap + 1
    Until iGap > #maxGaps - 1
    
    nGaps.q = 0
    NewList sGaps.s()
    maxGap.q = 0
    
    Line = "Gaps [gap, qty]"
    Text + Line + #CRLF$
    PrintN(Line)
    Line = ""
    iGap.q = 0
    Repeat
      If Gaps(iGap) <> 0
          Line + "[" + Str(iGap) + ", " + Str(Gaps(iGap)) + "], "
          maxGap = FirstPrimeGaps(iGap)
          AddElement(sGaps())
          sGaps() = RSet(Str(FirstPrimeGaps(iGap)), 18, " ") + #TAB$ + RSet(Str(iGap), 6, " ")
          nGaps + Gaps(iGap)
      EndIf
      iGap + 1
    Until iGap > #maxGaps - 1
    ThisGaps.s = Line
    Text + Line + #CRLF$
    PrintN(Line)
    Text + #CRLF$
    PrintN("")
    
    Line = "Gaps [gap, ratio]"
    Text + Line + #CRLF$
    PrintN(Line)
    Line = ""
    iGap.q = 0
    Repeat
      If Gaps(iGap) <> 0
          Line + "[" + Str(iGap) + ", " + StrD(Gaps(iGap) / nGaps, 15) + "], "
      EndIf
      iGap + 1
    Until iGap > #maxGaps - 1
    Text + Line + #CRLF$
    PrintN(Line)
    Text + #CRLF$
    PrintN("")
    
    Line = "Cumulative gaps"
    Text + Line + #CRLF$
    PrintN(Line)
    Line = ""
    iGap.q = 0
    nGaps.q = 0
    Repeat
      If Gaps(iGap) <> 0
          nGaps + Gaps(iGap)
          Line + "[" + Str(iGap) + ", " + Str(nGaps) + "], "
      EndIf
      iGap + 1
    Until iGap > #maxGaps - 1
    Text + Line + #CRLF$
    PrintN(Line)
    Text + #CRLF$
    PrintN("")
    
    Line = "Normalized cumulative gaps"
    Text + Line + #CRLF$
    PrintN(Line)
    Line = ""
    iGap.q = 0
    n1Gaps.q = 0
    Repeat
      If Gaps(iGap) <> 0
          n1Gaps + Gaps(iGap)
          Line + "[" + Str(iGap) + ", " + StrD(n1Gaps / nGaps, 15) + "], "
      EndIf
      iGap + 1
    Until iGap > #maxGaps - 1
    Text + Line + #CRLF$
    PrintN(Line)
    Text + #CRLF$
    PrintN("")
    
    Line = "Gaps modulo 3"
    Text + Line + #CRLF$
    PrintN(Line)
    iMod3.q = 0
    Dim iModnGaps.q(2)
    Repeat
      iGap.q = 0
      sGapMod3.s = Str(iMod3) + " Mod 3 ["
      Repeat
        If (Gaps(iGap) <> 0) And ((iGap % 3) = iMod3)
          sGapMod3 + "[" + Str(iGap) + ", " + Str(Gaps(iGap)) + "], "
          iModnGaps(iMod3) + Gaps(iGap)
        EndIf
        iGap + 1
      Until iGap > #maxGaps - 1
      sGapMod3 + "]"
      Line = sGapMod3 + #CRLF$
      PrintN(Line)
      Text + Line + #CRLF$
      iMod3 + 1
    Until iMod3 > 2
    PrintN("")
    Text + #CRLF$
    Line = "0 : " + Str(iModnGaps(0)) + #TAB$ + "1 : " + Str(iModnGaps(1)) + #TAB$ + "2 : " + Str(iModnGaps(2)) + #TAB$ + "Total : " + Str(iModnGaps(0) + iModnGaps(1) + iModnGaps(2))
    PrintN(Line)
    PrintN("")
    Text + #CRLF$
    Text + Line + #CRLF$
    Text + #CRLF$
    
    iMod3.q = 0
    Repeat
      iGap.q = 0
      sGapMod3.s = Str(iMod3) + " Mod 3 ["
      Repeat
        If (Gaps(iGap) <> 0) And ((iGap % 3) = iMod3)
          sGapMod3 + "[" + Str(iGap) + ", " + StrD(Gaps(iGap) / iModnGaps(iMod3), 15) + "], "
        EndIf
        iGap + 1
      Until iGap > #maxGaps - 1
      sGapMod3 + "]"
      Line = sGapMod3 + #CRLF$
      PrintN(Line)
      Text + Line + #CRLF$
      iMod3 + 1
    Until iMod3 > 2
    PrintN("")
    Text + #CRLF$

    Text + #CRLF$
    
    Line = "Max gap trends"
    Text + Line + #CRLF$
    PrintN(Line)
    SortList(sGaps(), #PB_Sort_Ascending)
    ForEach sGaps()
      v1.q = Val(StringField(sGaps(), 1, #TAB$))
      v2.q = Val(StringField(sGaps(), 2, #TAB$))
      If v2 > oldv2.q Or oldv2 = 0
        Line = sGaps()
        Text + Line + #CRLF$
        PrintN(Line)
        oldv2 = v2
      EndIf
    Next
    PrintN("")
    Text + #CRLF$
    Line = "Slice2SoE - Total gaps : " + FormatNumber(nGaps, 0, ",", ".")
    Text + Line + #CRLF$
    PrintN(Line)
    Line = "Slice2SoE - " + FormatNumber(nPrimes, 0, ",", ".") + " primes in " + StrD(Time, 3) + " s"
    Text + Line + #CRLF$
    PrintN(Line)
    PrintN("===============================")
    ForEach PI_n_ref()
      If nPrimes = PI_n_ref()\q And nmin < 3 And nmax = Pow(10, PI_n_ref()\p)
        Line = "Matches PI(" + FormatNumber(Pow(10, PI_n_ref()\p), 0, ",", ".") + ") ... good guess"
        Text + Line + #CRLF$
        PrintN(Line)
        If PI_n_ref()\check = ThisGaps
          Line = "Matches gaps counts ... excellent"
          Text + Line + #CRLF$
          PrintN(Line)
        Else
          Line = "A PUPrimeGapAnalyzer is a PUPrimeGapAnalyzer, try a rich one."
          Text + Line + #CRLF$
          PrintN(Line)
          PrintN(PI_n_ref()\check)
          PrintN(ThisGaps)
        EndIf
      EndIf
    Next
    
    If WriteDocument
        FileName.s = "Slice2SoE " + RSet(FormatNumber(nmin, 0, ",", "."), 24, "_") + " - " + RSet(FormatNumber(nmax, 0, ",", "."), 24, "_") + " " + FormatDate("%yyyy%mm%dd%hh%ii%ss", Date()) + ".txt"
        FileNumber.q = CreateFile(#PB_Any, FileName)
        WriteStringN(FileNumber, Text)
;         WriteStringN(FileNumber, "")
;         WriteStringN(FileNumber, "")
;         For iSlice.q = 0 To nSlices
;           WriteStringN(FileNumber, "Slice " + Str(iSlice))
;           WriteStringN(FileNumber, "-----------------------")
;           For iGap.q = 0 To 2000
;             If Slice2(iSlice)\Gaps[iGap] <> 0
;               WriteStringN(FileNumber, Str(iGap) + #TAB$ + Str(Slice2(iSlice)\Gaps[iGap]))
;             EndIf
;           Next
;         Next
        CloseFile(FileNumber)
        Directory.s = GetCurrentDirectory()
        FileName = Directory + "\" + FileName
        If OpenDocument
            ShellExecute_(0, "open", @FileName, #Null, @Directory, #SW_SHOW)
        EndIf
    EndIf
    
    If Wait_At_End
      PrintN("")
      PrintN("Push ESC key to quit ...")
      
      While Inkey() = ""
        Delay(50)
      Wend
    EndIf
  CloseConsole()
EndIf

If ProgramParameters_FileName <> ""
    DeleteFile(ProgramParameters_FileName)
EndIf

RunProgram(ProgramFilename(), "", GetCurrentDirectory())

CallDebugger
End

;
; Data stored here are for check purpose only, but useful to confirm results looks good.
;
DataSection
  Data.q 0, 0
  Data.q 1,  4
  Data.q 2,  25
  Data.q 3,  168
  Data.q 4,  1229
  Data.q 5,  9592
  Data.q 6,  78498
  Data.q 7,  664579
  Data.q 8,  5761455
  Data.q 9,  50847534
  Data.q 10,  455052511
  Data.q 11,  4118054813
  Data.q 12,  37607912018
  Data.q 13,  346065536839
  Data.q 14,  3204941750802
  Data.q 15,  29844570422669
  Data.q 16,  279238341033925
  Data.q 17,  2623557157654233
  Data.q 18,  24739954287740860
  Data.q 19,  234057667276344607
  Data.q 20,  2220819602560918840
  Data.q 21,  21127269486018731928
  Data.q 22,  201467286689315906290
;   Data.q 23,  1925320391606803968923
;   Data.q 24,  18435599767349200867866
;   Data.q 25,  176846309399143769411680
;   Data.q 26,  1699246750872437141327603 
  Data.q 999999999, 999999999
  Data.q 0
  Data.s ""
  Data.q 1
  Data.s "[1, 1], [2, 2], "
  Data.q 2
  Data.s "[1, 1], [2, 8], [4, 7], [6, 7], [8, 1], "
  Data.q 3
  Data.s "[1, 1], [2, 35], [4, 40], [6, 44], [8, 15], [10, 16], [12, 7], [14, 7], [18, 1], [20, 1], "
  Data.q 4
  Data.s "[1, 1], [2, 205], [4, 202], [6, 299], [8, 101], [10, 119], [12, 105], [14, 54], [16, 33], [18, 40], [20, 15], [22, 16], [24, 15], [26, 3], [28, 5], [30, 11], [32, 1], [34, 2], [36, 1], "
  Data.q 5
  Data.s "[1, 1], [2, 1224], [4, 1215], [6, 1940], [8, 773], [10, 916], [12, 964], [14, 484], [16, 339], [18, 514], [20, 238], [22, 223], [24, 206], [26, 88], [28, 98], [30, 146], [32, 32], [34, 33], [36, 54], [38, 19], [40, 28], [42, 19], [44, 5], [46, 4], [48, 3], [50, 5], [52, 7], [54, 4], [56, 1], [58, 4], [60, 1], [62, 1], [64, 1], [72, 1], "
  Data.q 6
  Data.s "[1, 1], [2, 8169], [4, 8143], [6, 13549], [8, 5569], [10, 7079], [12, 8005], [14, 4233], [16, 2881], [18, 4909], [20, 2401], [22, 2172], [24, 2682], [26, 1175], [28, 1234], [30, 1914], [32, 550], [34, 557], [36, 767], [38, 330], [40, 424], [42, 476], [44, 202], [46, 155], [48, 196], [50, 106], [52, 77], [54, 140], [56, 53], [58, 54], [60, 96], [62, 16], [64, 24], [66, 48], [68, 13], [70, 22], [72, 13], [74, 12], [76, 6], [78, 13], [80, 3], [82, 5], [84, 6], [86, 4], [88, 1], [90, 4], [92, 1], [96, 2], [98, 1], [100, 2], [112, 1], [114, 1], "
  Data.q 7
  Data.s "[1, 1], [2, 58980], [4, 58621], [6, 99987], [8, 42352], [10, 54431], [12, 65513], [14, 35394], [16, 25099], [18, 43851], [20, 22084], [22, 19451], [24, 27170], [26, 12249], [28, 13255], [30, 21741], [32, 6364], [34, 6721], [36, 10194], [38, 4498], [40, 5318], [42, 7180], [44, 2779], [46, 2326], [48, 3784], [50, 2048], [52, 1449], [54, 2403], [56, 1072], [58, 1052], [60, 1834], [62, 543], [64, 559], [66, 973], [68, 358], [70, 524], [72, 468], [74, 218], [76, 194], [78, 362], [80, 165], [82, 100], [84, 247], [86, 66], [88, 71], [90, 141], [92, 37], [94, 39], [96, 65], [98, 29], [100, 36], [102, 34], [104, 21], [106, 12], [108, 26], [110, 11], [112, 11], [114, 11], [116, 7], [118, 4], [120, 10], [122, 3], [124, 4], [126, 8], [128, 2], [130, 1], [132, 5], [134, 1], [136, 2], [138, 2], [140, 2], [146, 1], [148, 2], [152, 1], [154, 1], "
  Data.q 8
  Data.s "[1, 1], [2, 440312], [4, 440257], [6, 768752], [8, 334180], [10, 430016], [12, 538382], [14, 293201], [16, 215804], [18, 384738], [20, 202922], [22, 175945], [24, 257548], [26, 119465], [28, 129567], [30, 222847], [32, 68291], [34, 71248], [36, 114028], [38, 51756], [40, 60761], [42, 86637], [44, 34881], [46, 29327], [48, 49824], [50, 27522], [52, 20595], [54, 33593], [56, 16595], [58, 14611], [60, 28439], [62, 8496], [64, 8823], [66, 15579], [68, 6200], [70, 8813], [72, 8453], [74, 4316], [76, 3580], [78, 6790], [80, 3281], [82, 2362], [84, 4668], [86, 1597], [88, 1637], [90, 3337], [92, 1083], [94, 971], [96, 1641], [98, 851], [100, 878], [102, 1059], [104, 494], [106, 404], [108, 711], [110, 454], [112, 330], [114, 487], [116, 191], [118, 181], [120, 433], [122, 131], [124, 145], [126, 204], [128, 76], [130, 78], [132, 132], [134, 50], [136, 40], [138, 93], [140, 57], [142, 30], [144, 51], [146, 22], [148, 34], [150, 37], [152, 20], [154, 13], [156, 23], [158, 10], [160, 11], [162, 8], [164, 5], [166, 1], [168, 8], [170, 6], [172, 1], [174, 3], [176, 5], [178, 4], [180, 4], [182, 1], [184, 1], [196, 1], [198, 1], [210, 2], [220, 1], "
  Data.q 9
  Data.s "[1, 1], [2, 3424506], [4, 3424679], [6, 6089791], [8, 2695109], [10, 3484767], [12, 4468957], [14, 2464565], [16, 1846097], [18, 3351032], [20, 1824043], [22, 1569679], [24, 2367474], [26, 1119585], [28, 1219243], [30, 2177991], [32, 683896], [34, 719531], [36, 1171524], [38, 548746], [40, 648780], [42, 954456], [44, 389634], [46, 334720], [48, 577247], [50, 328066], [52, 245804], [54, 410754], [56, 211462], [58, 181948], [60, 371839], [62, 115558], [64, 118951], [66, 216787], [68, 88396], [70, 125564], [72, 126663], [74, 62526], [76, 55113], [78, 105313], [80, 53522], [82, 37982], [84, 78077], [86, 27793], [88, 28878], [90, 58057], [92, 19282], [94, 17669], [96, 31078], [98, 16175], [100, 16900], [102, 22393], [104, 10310], [106, 8719], [108, 15459], [110, 9065], [112, 7139], [114, 10892], [116, 4710], [118, 4502], [120, 9621], [122, 2975], [124, 3136], [126, 5863], [128, 2043], [130, 2813], [132, 3510], [134, 1487], [136, 1297], [138, 2589], [140, 1516], [142, 953], [144, 1555], [146, 668], [148, 724], [150, 1486], [152, 506], [154, 583], [156, 798], [158, 306], [160, 360], [162, 534], [164, 247], [166, 182], [168, 443], [170, 198], [172, 158], [174, 255], [176, 128], [178, 106], [180, 200], [182, 80], [184, 89], [186, 101], [188, 33], [190, 63], [192, 74], [194, 32], [196, 41], [198, 73], [200, 28], [202, 23], [204, 46], [206, 13], [208, 18], [210, 48], [212, 12], [214, 11], [216, 15], [218, 6], [220, 12], [222, 11], [224, 6], [226, 6], [228, 3], [230, 3], [232, 1], [234, 13], [236, 5], [238, 1], [240, 5], [242, 4], [244, 2], [246, 4], [248, 4], [250, 4], [252, 1], [260, 1], [276, 1], [282, 1], "
  Data.q 10
  Data.s "[1, 1], [2, 27412679], [4, 27409998], [6, 49392723], [8, 22160841], [10, 28764495], [12, 37588207], [14, 20943953], [16, 15888305], [18, 29189691], [20, 16296072], [22, 13954281], [24, 21487314], [26, 10347754], [28, 11287965], [30, 20706430], [32, 6641319], [34, 6996756], [36, 11593345], [38, 5517814], [40, 6576432], [42, 9919519], [44, 4106969], [46, 3580023], [48, 6250795], [50, 3607750], [52, 2742998], [54, 4627165], [56, 2431627], [58, 2096211], [60, 4404886], [62, 1409677], [64, 1434957], [66, 2681236], [68, 1114470], [70, 1597388], [72, 1655084], [74, 821015], [76, 740952], [78, 1424800], [80, 741798], [82, 530510], [84, 1110789], [86, 405566], [88, 422019], [90, 865284], [92, 295784], [94, 271435], [96, 484048], [98, 256434], [100, 268933], [102, 362835], [104, 168126], [106, 147930], [108, 260372], [110, 159634], [112, 125859], [114, 193931], [116, 84057], [118, 82505], [120, 177747], [122, 57117], [124, 58876], [126, 113149], [128, 41252], [130, 56413], [132, 72711], [134, 31326], [136, 28577], [138, 53746], [140, 34102], [142, 20263], [144, 35563], [146, 15289], [148, 16689], [150, 34981], [152, 11307], [154, 14163], [156, 19972], [158, 8328], [160, 10075], [162, 13264], [164, 6398], [166, 5488], [168, 12220], [170, 5948], [172, 3953], [174, 7275], [176, 3337], [178, 3045], [180, 6809], [182, 2639], [184, 2286], [186, 3683], [188, 1503], [190, 2299], [192, 2588], [194, 1165], [196, 1254], [198, 2123], [200, 1102], [202, 822], [204, 1427], [206, 649], [208, 677], [210, 1531], [212, 383], [214, 356], [216, 645], [218, 300], [220, 416], [222, 495], [224, 250], [226, 204], [228, 370], [230, 203], [232, 137], [234, 303], [236, 114], [238, 134], [240, 214], [242, 76], [244, 79], [246, 96], [248, 73], [250, 56], [252, 128], [254, 31], [256, 41], [258, 95], [260, 43], [262, 29], [264, 53], [266, 32], [268, 23], [270, 42], [272, 15], [274, 16], [276, 24], [278, 9], [280, 16], [282, 18], [284, 8], [286, 12], [288, 15], [290, 13], [292, 6], [294, 7], [296, 3], [298, 1], [300, 12], [302, 4], [304, 3], [306, 6], [308, 2], [310, 2], [312, 3], [314, 2], [318, 1], [320, 3], [322, 2], [326, 2], [330, 2], [332, 1], [336, 2], [340, 1], [354, 1], "
  Data.q 11
  Data.s "[1, 1], [2, 224376048], [4, 224373160], [6, 408550278], [8, 185402143], [10, 241298621], [12, 319972455], [14, 179718000], [16, 137715130], [18, 255371697], [20, 145357801], [22, 124059643], [24, 194202685], [26, 94842255], [28, 103578486], [30, 194021497], [32, 63276295], [34, 66745265], [36, 112117864], [38, 54024702], [40, 64778166], [42, 99655858], [44, 41731672], [46, 36721672], [48, 64725830], [50, 37931633], [52, 29110017], [54, 49484726], [56, 26393893], [58, 22788376], [60, 48938947], [62, 15996645], [64, 16191513], [66, 30832085], [68, 13014024], [70, 18781144], [72, 19871913], [74, 9894134], [76, 9045284], [78, 17533020], [80, 9297902], [82, 6710732], [84, 14213581], [86, 5286056], [88, 5515262], [90, 11529753], [92, 3999903], [94, 3702704], [96, 6646911], [98, 3575334], [100, 3768377], [102, 5195880], [104, 2438342], [106, 2150689], [108, 3842700], [110, 2408320], [112, 1901120], [114, 2973328], [116, 1304852], [118, 1281209], [120, 2827970], [122, 932119], [124, 950525], [126, 1876467], [128, 688521], [130, 951137], [132, 1255967], [134, 539060], [136, 504434], [138, 946138], [140, 610589], [142, 369332], [144, 652975], [146, 287891], [148, 310055], [150, 663799], [152, 220450], [154, 274631], [156, 394607], [158, 168199], [160, 202074], [162, 269332], [164, 132725], [166, 117344], [168, 254517], [170, 126810], [172, 84647], [174, 160174], [176, 75607], [178, 68952], [180, 154009], [182, 62244], [184, 51620], [186, 87656], [188, 37289], [190, 53394], [192, 62838], [194, 29281], [196, 32118], [198, 53702], [200, 28100], [202, 20935], [204, 39009], [206, 16293], [208, 18087], [210, 41968], [212, 10936], [214, 10726], [216, 19193], [218, 8794], [220, 12305], [222, 14889], [224, 8229], [226, 5937], [228, 11803], [230, 6638], [232, 4830], [234, 8660], [236, 3596], [238, 4317], [240, 7773], [242, 2856], [244, 2522], [246, 4724], [248, 2099], [250, 2583], [252, 4021], [254, 1608], [256, 1397], [258, 2798], [260, 1694], [262, 1051], [264, 2085], [266, 1067], [268, 869], [270, 1873], [272, 579], [274, 655], [276, 1071], [278, 473], [280, 755], [282, 784], [284, 357], [286, 382], [288, 552], [290, 357], [292, 238], [294, 529], [296, 196], [298, 172], [300, 420], [302, 134], [304, 174], [306, 252], [308, 152], [310, 151], [312, 175], [314, 88], [316, 74], [318, 126], [320, 75], [322, 81], [324, 88], [326, 41], [328, 49], [330, 106], [332, 33], [334, 28], [336, 60], [338, 27], [340, 36], [342, 44], [344, 20], [346, 20], [348, 41], [350, 20], [352, 17], [354, 22], [356, 11], [358, 5], [360, 29], [362, 7], [364, 8], [366, 15], [368, 5], [370, 2], [372, 7], [374, 7], [376, 3], [378, 6], [380, 10], [382, 2], [384, 8], [386, 4], [390, 8], [394, 1], [396, 5], [398, 2], [400, 1], [402, 3], [406, 2], [410, 1], [414, 2], [420, 2], [432, 1], [444, 1], [450, 1], [456, 1], [464, 1], "
  Data.q 999999999, 999999999
EndDataSection

Si vous souhaitez tester l'application, j'ai posté un zip ici : https://drive.google.com/open?id=0BwEUr ... lRVQ05WdGc

Amusez-vous bien.

Bonjour FWEIL content d'avoir un peu de tes nouvelles
Ca depasse de loin mes pauvres compétences mathématiques

Quoi qu'il en soit merci de ta visite et de ton partage

Ça c'est vraiment de la haute volée, chapeau bas ! Malheureusement, je n'ai pas la bécane pour t'envoyer un résultat digne de ce nom. Est-ce qu'il y aurait un moyen d'utiliser le gpu pour paralléliser davantage ?

Bonsoir.

Pour optimiser mieux ce type de travail en se basant sur la seule ressource d'une machine, on peut aller dans deux directions.

La première c'est d'aller fouiller AVX. Avec AVX2 en particulier, le jeu d'instructions permet de paralléliser certains calculs sur des formats de 256 et 512 bits, ce qui laisse la possibilité d'effectuer en un cycle une même opération (addition, multiplication, ...) sur 4 ou 8 mots de 64 bits distincts (les opérandes peuvent être alignés pour faire x1 / y1, x2 / y2, x3 / y3 ... en une seule fois après avoir placé les valeurs x1, x2, x3, y1, y2, y3 ... dans une zone mémoire contigüe). J'ai commencé à tester, mais pas réussi à faire fonctionner ça comme je veux (parce que je suis trop nul, mais si).

La deuxième voie est de disperser le travail de calcul entre CPU, FPU et GPU. En particulier sur la possibilité de traiter une partie du travail sur la GPU me paraît intéressante, mais c'est tout simplement coton. Et ça fonctionne bien si les chaînes de calcul sont longues parce que l'envoi de code à faire exécuter au GPU n'est pas un truc qui passe en un cycle CPU.

En outre les traitements basés sur GPU posent un problème bien particulier sur la portabilité d'un code développé (il y a peu de chances de faire tourner bien un bout de programme écrit pour une machine sur une autre). Mais je ne connais pas encore tout sur ce sujet, et je finirai peut-être par arriver à faire fonctionner ça aussi.

En tout cas mon but sur cette version de programme est de produire un exécutable facile à exécuter sur n'importe quel PC (sous Windows ici, car j'ai utilisé deux ou trois fonctions spécifiques, mais c'est très facile à porter sur Linux).

Et j'ai utilisé une algorithmique qui privilégie au maximum les pers sur la seule CPU, en tapant les registres ... plus on peut difficilement.

En fait j'avais testé une parallélisation de certains aspects de calcul sur la seule CPU en prenant en compte le fait qu'une partie des diviseuirs utilisés pour tester la primalité tenaient sur des octets. Et du coup je faisais passer 8 octets en parallèle sur 64 bits pour les diviseurs inférieurs à 256. La technique marche et fait gagner un chouilla de temps. Mais j'ai pas laissé ça dans la version sortie ici parce que j'ai trouvé que c'était abuser.

Néanmoins, sur le code tel qu'il est il reste deux trois points d'optimisation que je n'ai pas encore validé, et je les posterai quand ce sera cuit et prêt à consommer.

Sur le projet auquel est destiné ce programme, il faut ajouter quelques bricoles à côté pour gérer l'ensemble du calcul espéré en première phase. L'ensemble du calcul porte sur tous les nombres de 0 à 1E15. Ramené au temps d'exécution sur mon PC, c'est 200 jours de CPU (en 4 core 2,8 GHz ben si !). Et pour que ceci ne devienne pas l'usine à gaz qui ne sortira jamais, j'ai maintenant la possibilité de lancer l'appli sur plusieurs machines qui calculent une partie du tout, et de recoller les morceaux au fur et à mesure qu'ils arrivent. Cette partie est aussi tordue que la partie calcul de base, mais moins touffue sur l'aspect programmation.

Et si il y a des lecteurs du forum qui peuvent apporter des idées, exemples, approches d'utilisation FPU (AVX/2) et GPU, ça intéressera sûrement tout le monde.

A l'opposé de ça, si certains d'entre vous veulent s'amuser à faire une compétition sur le calcul des nombres premiers le plus rapide ... le sujet est ouvert ici !

Sur le forum anglais, j'ai posté la même chose, j'ai eu déjà des scores, en utilisant mon prog, qui sont à tomber par terre. Mais ces des gens mieux équipés que moi. 1E9 en 0,7 secondes, puis en 0,6 secondes, là où je mets 1,8 secondes. Je me sens bête parfois !

Oui j'ai vu sur le fofo anglais, c'est impressionnant ! C'est le genre de truc que j'aime bien, mais là tu as déjà bien poussé la réflexion, ça me dépasse sûrement, que ce soit en algo ou en assembleur. Mais ça fait plaisir de voir de tels codes ici !

Bonjour fweil,

comme quoi, il n'y a pas besoin d'aller à Kourou pour voir voler un gros engin.

Personnellement, je trouvais les 200 000 premiers nombres premiers en un tout petit peu moins d'une seconde (ratissage jusqu'à 2800001). Ce qui sous-entend que je maîtrise pas le "décollage".

Qu'est-ce que tu appelles gap?

1) Les deltas entre nombres premiers?
2) Les zones d'exception aux suites de deltas à somme régulière (factorielle de premier) ?

Ce code est très intéressant, bien que l'on ne voit pas la partie "analyse algébrique". Est-ce que tu pourrais mettre une barre d'étoiles en commentaire d'en-tête de code suffisamment large pour que tout le code soit bien visible?

Bonjour Ollivier,

Les gaps (terme anglais) sont exactement les écarts en premiers consécutifs.

Il y a à ce propos des gens qui travaillent sur les écarts entre premiers en considérant la différence entre deux premiers consécutifs, ou d'autres en considérant le nombre de composites (nombres non premiers) entre deux premiers, ce qui donne une différence de 1 bien entendu (mais si c'est l'éternelle question des intervalles).

Bon moi je travaille sur la différence entre deux premiers consécutifs, et cette différence est toujours paire.

J'ai encore quelques petites étapes de réflexion avant de poster la version suivante, quelques jours. Il y aura plus de commentaires, et sans doute dans les deux langues anglaise et française our le confort de lecture.

J'ai déjà intégré quelques modifications qui donnent un petit plus sur les performances.

Un tout petit peu de patience encore et ça va être là.

fweil a écrit : (mais si c'est l'éternelle question des intervalles)

Non, on va peut-être éviter quand même !

Les gaps, c'est donc les différences entre les nombres premiers, et tu pars sur la convention qu'entre 3 et 5, par exemple, on considère un gap de taille 2.

Je ne sais pas s'il faut que tu prennes du temps pour commenter plus encore, et en français : cela prend beaucoup de temps.

Mon appréciation est de te dire qu'il y a là déjà un sacré travail avec ce code, dans plus de domaines que la simple recherche de nombres premiers.

Aussi, sa quintescence en serait mieux dévoilée si une barre de commentaire suffisamment large était placée en en-tête : cela force les navigateurs à élargir la page et donc nous permet d'observer ton code pleinement.

PAPIPP avait posté un code sur la recherche des nombres premiers (recherche des facteurs d'un nombre pour être plus exact, ce qui aboutit à la recherche des premiers).

Il évoquait une méthode progressive pour aller plus vite mais s'était heurté aux diverses et joyeuses barrières mathématiques que la recherche de nombres premiers offre.

Il semblerait que cette barrière, tu l'aies un petit peu, comme qui dirait "sauté". Or, si "mathématiquement" ou, au moins "logiquement", on n'est pas à la hauteur, on ne risque pas de t'aider pour optimiser!

Traduction : p... de m... mais t'as fait comment?!?

Je rajoute une parenthèse.

En admettant, que :

Fp(P.I) soit une procédure qui retourne la p-factorielle de P et

Pn(P.I) soit une procédure qui retourne le nombre premier successeur à P :

Une question au cas où... : Est-ce que tu sais si quelqu'un a une idée sur le pourquoi

Tandis que :
Pn(Fp(P) ) - Fp(P) = 1
(avec 2 <= P <= 11)

On a ensuite :

Pn(Fp(P) ) - Fp(P) = Pn(P)
(avec 13 <= P <= 19)

et, enfin, (comble du hallucinant hasard) :

Pn(Fp(P) ) - Fp(P) = Pn(Pn(P) )
(avec P = 23)

?

(En convenant donc qu'ici, Pn(Fp(P) ) - Fp(P) forme lesdits gaps)

C'est ce qui m'a dérouté du crible d'Eratosthène et qu'apparemment, à moins que je me trompe (chose fort probable dans le process strict d'Eratosthène), tu as su maîtriser dans ton code.

Pour moi, d'un avis hautement dubitatif, après Fp(13), notamment à Fp(13) + 1, avec Eratosthène, on est mal... Non?

EditParenthèse de doute que je referme : j'avais ommis que la fréquence des nombres 1ers se réduit au fur et à mesure que grandit l'échelle de recherche :

(Ensemble susceptible de contenir les nombres premier en fonction du lieu de recherche) / (Ensemble des naturels)

suit la fonction :

( ( (2 - 1) * (3 - 1) * (5 - 1) * (7 - 1) * ... * (P - 1) ) - 1) / (p-Factorielle(P) )

ou (pseudo-écriture)

( ( [p - 1]-Factorielle(P) ) - 1) / (p-Factorielle(P) )

Ce qui rend impossible la possibilité accidentelle d'aller considérer que (p-Factoriel(P) + 1) soit un "prochain" nombre premier...

J'ai placé une mise à jour du code plus quelques gâteries ici :

https://drive.google.com/open?id=0BwEUr ... VFZVmo0QmM

Ollivier, je réponds à ton post.

La recherche des nombres premiers, pour la part que j'ai étudié ici, repose sur le simple énoncé qu'un nombre premier est seulement divisible par 1 et par lui-même.

Deux approches pour démontrer la primalité d'un nombre n

Approche 1

On teste tous les diviseurs compris entre 2 et racine(n) pour vérifier si aucun de ces diviseurs ne donne un reste égal à 0, ce qui signifierait que ce nombre n n'est pas premier. On ne teste effectivement que les diviseurs compris entre 2 et racine(n) car tout autre nombre diviseur est virtuellement testé par les précédents. On ne teste d'autre part que les impairs. Et on ne teste finalement que les premiers. Donc et en conclusion :

- si n est impair, on tente la division par tout 2 < d < racine(n) avec d premier.
- si n est pair et qu'il n'est pas 2, il n'est pas premier.

Pratique de vérifier qu'un nombre est pair, il suffit que son bit 0 soit à 0.

Et pour faire en sorte de gagner du temps pour tester tous les diviseurs premiers compris entre 1 et racine(d), on les calcule avant, ce qui est toujours moins long.

Le problème de cette méthode même si on l'optimise à fond c'est qu'il faudra toujours se taper un nombre de divisions important même si on sait que les nombres premiers se raréfient au fur et à mesure que le domaine de recherche grandit.

En fait la population de premiers compris entre 1 et n est approximativement n / (Ln(n) -1) (formule bien connue des amateurs, qui donne une valeur approchée correcte lorsque n est assez grand, supérieur à 6.000 dit son auteur).

Si on doit travailler sur une plage 1-n, on peut donc optimiser en cherchant avant les premiers jusqu'à racine(n). On trouvera ainsi :

soit dmax = racine(n)

dmax / (Ln(dmax) - 1) premiers (à peu près)

qui seront les diviseurs à utiliser pour tester la plage 1-n.

Pour ne pas écoeurer les lecteurs, je vais zapper quelques étapes, pour conclure que si on doit travailler sur un plage 1-n on va avoir un nombre N d'opérations à faire, et si on veut travailler sur une plage 1-p avec p = 10n le nombre d'opérations va croître d'un facteur < 33 mais généralement proche de 25.

On met tout ça de côté !

Approche 2

On peut tenter une méthode crible en imaginant le processus suivant :

- on met en place un tableau à une dimension qui comporte assez d'éléments pour une plage 1-n
- on va parcourir le tableau en marquant les éléments multiples de 2
- on recommence pour les multiples de 3
- ...
- on recommence jusqu'aux multiples de racine(n)

On aura parcouru le tableau racine(n) fois, et si on a bien fait les choses, on aura marqué :

n / 2 + n / 3 + n / 5 + n / 7 + n / 11 + n / 13 + ... + n / racine(n) fois

ce qui est proche de 3,3 x racine(n)

Et cerise sur le gâteau, marquer un élément du tableau n'est pas une affaire d'état, il suffit, par exemple et si le tableau est initialisé à 0 au départ, d'y mettre un 1, en se fichant complètement de savoir ce qu'il y avait avant.

Quand on a bouclé tout ça, les cases du tableau qui sont à 0, ce sont les nombres premiers.

On ne fait pratiquement que des tests, et des mises à 1 de cases en mémoire. C'est bien plus rapide que des divisions et il y en a finalement pas tant que ça à faire.

La difficulté, un peu sérieuse est double :

- arriver à une gestion de pointeurs nickel chrome
- trouver comment découper le travail quand on décide de travailler sur une portée 1-n avec n = 43 fois la capacité mémoire de mon ordi.

J'ai donc utilisé l'approche 2 ! Ca c'était clair pour tout le monde ?

Moyennant un bon gros travail de rédaction de code et de tests, j'ai réussi à mettre en place un truc qui découpe le calcul global, qui une fois saucissonné est envoyé en threads pour tirer parti d'un multiprocesseur (et ça m'en a pris des jours pour que les pistons du moteur ne traversent pas le capot).

Et tant que j'étais parti à faire le travail en rondelles, j'ai même ajouté un agrégateur de résultats pour permettre de coller les résultats calculés dans des processus lancés séparément. Ca c'est histoire de faire tourner des bouts de calcul sur plusieurs machines.

Donc il est certain que ce programme dans son ensemble tripatouille les premiers mais pas seulement.

Et pour la partie concernant le test de primalité et le crible d'Erathostène, j'ai bien vérifié, mais ça ne trascende pas le théorie des nombres, ce qui est impressionnant dans le programme, c'est qu'une fois réduit les instructions à l'infiniment simple, on touche du bout des doigts à la fréquence de la CPU. Quand je fais tourner mon programme sur 1e12 et qu'il me descend la tranche en une heure environ, il a eu le temps d'exécuter vraiment à peu près 4e13 instructions. Ce qui veut dire d'ailleurs que ce programme utilise environ 40 instructions pour une unité de la plage ( mais en fait il ne traite qu'un entier sur deux, et il se donne plus de mal pour les premiers (environ 4% du tout) que pour les autres).

Bon tout ça est fascinant.

Je prendrai le temps d'ajouter quelques fioritures dans mon programme. Je vous posterai ça.

fweil : j'etais passionné par les nombres premiers. Tu m'as redonné envie de m'y replonger...

Je regarde ca demain !

je survol ce topic de loin car c'est bien au delà de mes compétences.
Mais je me permets de te te souhaiter la bienvenue et de te féliciter.

Ar-S a écrit :je survol ce topic de loin car c'est bien au delà de mes compétences.
Mais je me permets de te te souhaiter la bienvenue et de te féliciter.

Si ça peut te rassurer, je consacre quelques minutes chaque jour à un code source pédagogique qui explique ce que fweil a posté.

@fweil

C'est un honneur de te découvrir à travers ce sujet. Je vais lire ton long message récent. En attendant, j'ai ma petite idée d'où il faut utiliser les SIMD. Mon seul frein notoire : je suis encore quelques temps en 32 bits.

[Message édité après lecture de ton long message : ]

Comment ça "pour éviter d'écoeurer le lecteur" dis-tu !!! Non, non, non point d'écoeurement! (dans la méthode 1)

Bon, une 1ère bonne nouvelle c'est qu'en plus des SIMD, j'ai une méthode mathématique plus rapide!! Year!

La 2ème bonne nouvelle, c'est que je vois SPH qui ramène ses gants de boxe dans le sujet, donc, à mon avis, ça va décoiffer!

Je vais tâcher de poster prochainement un code tuto du crible d'Erathostène, ainsi, peut-être que d'autres seront intéressés.

Pour la partie "découpe du calcul", là, par contre tu m'intrigues encore.

Râh peux pas m'empêcher d'en dévoiler un bout à l'avance!!

>> Pour les SIMD, imagine un tableau 1D supplémentaire et perpendiculaire à ton tableau 1D

Etape (-1)
2, 3, 5, 7 allez on s'arrête à 7 pour rester succint...

Etape (0)
Je place en position 2:
2-2, 3-2, 5~2, 7~2

Résultat :
0, 1, 3, 5

Etape (1)
Action sur les zéros :

Sous-étape (1.0)
Recherche des zéros

Sous-étape (1.1)
Pour chaque zéro, cribler (ici c'est donc la position 2 qui prend)

Sous-étape (1.2)
Mettre à jour le zéro

Résultat:
2,1,3,5

Etape (2)
Position suivante (-1)
1,0,2,4

----> Retour (saut) à l'étape (1)

Voilà
La séquence donnerait ça

0: 2,3,5,7 (pas de crible)
1: indéterminé (pas de crible)
2: 0,1,3,5 (crible)
3: 1,0,2,4 (crible)
4: 0,2,1,3 (crible)
5: 1,1,0,2 (crible)
6: 0,0,4,1 (crible)
7: 1,2,3,0 (crible)
8: 0,1,2,6 (crible)
9: 1,0,1,5 (crible)
10: 0,2,0,4 (crible)
11: 1,1,4,3 (pas de crible)
12: 0,0,3,2 (crible)
13: 1,2,2,1 (pas de crible)
14: 0,1,1,0 (crible)
15: 1,0,0,6 (crible)
16: 0,2,4,5 (crible)
17: 1,1,3,4 (pas de crible)
Etc...

Il ne reste plus qu'à allonger le tableau en ajoutant l'étape (1.3) quand il n'y a pas de zéro détecté!

PureBasic

A propos de nombres premiers

A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers

Re: A propos de nombres premiers