Nice to write some words here.
I am working, among other tricks, on a project consisting in editing some statistics about "primes gaps". The gaps between primes does matter to some scientists I know, and they asked me for some help in designing fast process to embed their database.
The project is limited to a 64 bits arithmetic, meaning that we only compute numbers up to 1e18, which is enough right now for the final results we want to get.
In this way, I wrote some code to make possible to calculate primes, count gaps between primes, and edit statistics about gaps quantities, normalized results (% vs quantities), and so on.
I also imagined a way to split a big work in slices, because covering all the scope from 0 to 1e18 is basically not accessible to humans working with PCs.
I actually limited the first step of results to 1e15 as a upper bound of computing.
I tested on a scope 1-1e12, for example, in the early code versions, and it was not possible to get end of process below a month. Now 1e12 is processed in an hour on my usual desktop.
It uses multithreading to take advantages with multicore CPUs.
To process very large scopes, like 1e15, it would take about 200 days to go. This is the main reason that made me writing something to split the work in shortest ranges that some other tools can merge at the end.
The program which is copied here contain all that stuff to ask you which range of numbers you want to test (default proposes a reasonable range 1-1E9).
They are some options to work with it, but this is another discussion not really helpful in Tricks 'n' Tips.
If you have any interest in more information about the project I am talking, or if you want more specific explanations except about the code which can be commented here, you can PM me. Also if you think you have some good time of PC to share with me, I can provide some part of computing in this project (just remind there is neither glory nor money in participating to this
).Well, not sure this preamble is clear and understandable, but here is the code ! I will maybe update it in the future at the same place for PureBasic coders.
Code: Select all
;
; 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
Enjoy.
