Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Zhuowen Zhao 2017-08-18 12:00:24 -04:00
commit 690edaa8dc
9 changed files with 392 additions and 404 deletions

View File

@ -1 +1 @@
v2.0.1-833-ga28b4b3 v2.0.1-859-ge8723bc

View File

@ -41,7 +41,7 @@ parser.add_option('-f','--formula',
parser.add_option('-c','--condition', parser.add_option('-c','--condition',
dest = 'condition', metavar='string', dest = 'condition', metavar='string',
help = 'condition to filter rows') help = 'condition to alter existing column data')
parser.set_defaults(condition = None, parser.set_defaults(condition = None,
) )
@ -77,29 +77,28 @@ for name in filenames:
# --------------------------------------- evaluate condition --------------------------------------- # --------------------------------------- evaluate condition ---------------------------------------
if options.condition is not None: if options.condition is not None:
interpolator = []
condition = options.condition # copy per file, since might be altered inline condition = options.condition # copy per file, since might be altered inline
breaker = False breaker = False
for position,operand in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups for position,(all,marker,column) in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups
condition = condition.replace('#'+operand[0]+'#', idx = table.label_index(column)
{ '': '{%i}'%position, dim = table.label_dimension(column)
's#':'"{%i}"'%position}[operand[1]]) if idx < 0 and column not in specials:
if operand[2] in specials: # special label damask.util.croak('column "{}" not found.'.format(column))
interpolator += ['specials["%s"]'%operand[2]]
else:
try:
interpolator += ['%s(table.data[%i])'%({ '':'float',
's#':'str'}[operand[1]],
table.label_index(operand[2]))] # could be generalized to indexrange as array lookup
except:
damask.util.croak('column "{}" not found.'.format(operand[2]))
breaker = True breaker = True
else:
if column in specials:
replacement = 'specials["{}"]'.format(column)
elif dim == 1: # scalar input
replacement = '{}(table.data[{}])'.format({ '':'float',
's#':'str'}[marker],idx) # take float or string value of data column
elif dim > 1: # multidimensional input (vector, tensor, etc.)
replacement = 'np.array(table.data[{}:{}],dtype=float)'.format(idx,idx+dim) # use (flat) array representation
condition = condition.replace('#'+all+'#',replacement)
if breaker: continue # found mistake in condition evaluation --> next file if breaker: continue # found mistake in condition evaluation --> next file
evaluator_condition = "'" + condition + "'.format(" + ','.join(interpolator) + ")"
# ------------------------------------------ build formulas ---------------------------------------- # ------------------------------------------ build formulas ----------------------------------------
evaluator = {} evaluator = {}
@ -162,7 +161,7 @@ for name in filenames:
# -------------------------------------- evaluate formulas ----------------------------------------- # -------------------------------------- evaluate formulas -----------------------------------------
if options.condition is None or eval(eval(evaluator_condition)): # condition for veteran replacement fulfilled if options.condition is None or eval(condition): # condition for veteran replacement fulfilled
for veteran in veterans: # evaluate formulas that overwrite for veteran in veterans: # evaluate formulas that overwrite
table.data[table.label_index(veteran): table.data[table.label_index(veteran):
table.label_index(veteran)+table.label_dimension(veteran)] = \ table.label_index(veteran)+table.label_dimension(veteran)] = \

View File

@ -14,7 +14,7 @@ scriptID = ' '.join([scriptName,damask.version])
# -------------------------------------------------------------------- # --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Add data of selected column(s) from (first) row of second ASCIItable that shares the linking column value. Add data of selected column(s) from (first) row of linked ASCIItable that shares the linking column value.
""", version = scriptID) """, version = scriptID)
@ -25,7 +25,7 @@ parser.add_option('--link',
parser.add_option('-l','--label', parser.add_option('-l','--label',
dest = 'label', dest = 'label',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'column label(s) to be appended') help = 'column label(s) to add from linked ASCIItable')
parser.add_option('-a','--asciitable', parser.add_option('-a','--asciitable',
dest = 'asciitable', dest = 'asciitable',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
@ -69,7 +69,7 @@ for name in filenames:
try: table = damask.ASCIItable(name = name, try: table = damask.ASCIItable(name = name,
buffered = False) buffered = False)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,"{} <== {}".format(name,options.asciitable))
# ------------------------------------------ read header ------------------------------------------ # ------------------------------------------ read header ------------------------------------------

View File

@ -51,7 +51,7 @@ parser.add_option('-c','--condition',
dest = 'condition', metavar='string', dest = 'condition', metavar='string',
help = 'condition to filter rows') help = 'condition to filter rows')
parser.set_defaults(condition = '', parser.set_defaults(condition = None,
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -98,23 +98,29 @@ for name in filenames:
else: else:
order = range(len(labels)) # maintain original order of labels order = range(len(labels)) # maintain original order of labels
interpolator = [] # --------------------------------------- evaluate condition ---------------------------------------
condition = options.condition # copy per file, might be altered if options.condition is not None:
for position,operand in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups condition = options.condition # copy per file, since might be altered inline
condition = condition.replace('#'+operand[0]+'#', breaker = False
{ '': '{{{}}}' .format(position),
's#':'"{{{}}}"'.format(position)}[operand[1]])
if operand[2] in specials: # special label ?
interpolator += ['specials["{}"]'.format(operand[2])]
else:
try:
interpolator += ['{}(table.data[{}])'.format({ '':'float',
's#':'str'}[operand[1]],
table.label_index(operand[2]))]
except:
parser.error('column "{}" not found...\n'.format(operand[2]))
evaluator = "'" + condition + "'.format(" + ','.join(interpolator) + ")" for position,(all,marker,column) in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups
idx = table.label_index(column)
dim = table.label_dimension(column)
if idx < 0 and column not in specials:
damask.util.croak('column "{}" not found.'.format(column))
breaker = True
else:
if column in specials:
replacement = 'specials["{}"]'.format(column)
elif dim == 1: # scalar input
replacement = '{}(table.data[{}])'.format({ '':'float',
's#':'str'}[marker],idx) # take float or string value of data column
elif dim > 1: # multidimensional input (vector, tensor, etc.)
replacement = 'np.array(table.data[{}:{}],dtype=float)'.format(idx,idx+dim) # use (flat) array representation
condition = condition.replace('#'+all+'#',replacement)
if breaker: continue # found mistake in condition evaluation --> next file
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
@ -129,7 +135,7 @@ for name in filenames:
outputAlive = True outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
specials['_row_'] += 1 # count row specials['_row_'] += 1 # count row
if condition == '' or eval(eval(evaluator)): # valid row ? if options.condition is None or eval(condition): # valid row ?
table.data = [table.data[position] for position in positions] # retain filtered columns table.data = [table.data[position] for position in positions] # retain filtered columns
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line

View File

@ -52,7 +52,7 @@ parser.set_defaults(data = [],
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
if not options.vtk: parser.error('No VTK file specified.') if not options.vtk: parser.error('no VTK file specified.')
if not os.path.exists(options.vtk): parser.error('VTK file does not exist.') if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
if os.path.splitext(options.vtk)[1] == '.vtp': if os.path.splitext(options.vtk)[1] == '.vtp':
@ -66,16 +66,16 @@ elif os.path.splitext(options.vtk)[1] == '.vtk':
reader.Update() reader.Update()
Polydata = reader.GetPolyDataOutput() Polydata = reader.GetPolyDataOutput()
else: else:
parser.error('Unsupported VTK file type extension.') parser.error('unsupported VTK file type extension.')
Npoints = Polydata.GetNumberOfPoints() Npoints = Polydata.GetNumberOfPoints()
Ncells = Polydata.GetNumberOfCells() Ncells = Polydata.GetNumberOfCells()
Nvertices = Polydata.GetNumberOfVerts() Nvertices = Polydata.GetNumberOfVerts()
if Npoints != Ncells or Npoints != Nvertices: if Npoints != Ncells or Npoints != Nvertices:
parser.error('Number of points, cells, and vertices in VTK differ from each other.') parser.error('number of points, cells, and vertices in VTK differ from each other.')
damask.util.croak('{}: {} points, {} vertices, and {} cells...'.format(options.vtk,Npoints,Nvertices,Ncells)) damask.util.croak('{}: {} points/vertices/cells...'.format(options.vtk,Npoints))
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
@ -97,16 +97,19 @@ for name in filenames:
VTKarray = {} VTKarray = {}
active = defaultdict(list) active = defaultdict(list)
for datatype,dimension,label in [['data',99,options.data], for datatype,dimension,label in [['data',0,options.data],
['tensor',9,options.tensor], ['tensor',9,options.tensor],
['color' ,3,options.color], ['color' ,3,options.color],
]: ]:
for i,dim in enumerate(table.label_dimension(label)): for i,dim in enumerate(table.label_dimension(label)):
me = label[i] me = label[i]
if dim == -1: remarks.append('{} "{}" not found...'.format(datatype,me)) if dim == -1: remarks.append('{} "{}" not found...'.format(datatype,me))
elif dim > dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension)) elif dimension > 0 \
and dim != dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension))
else: else:
remarks.append('adding {} "{}"...'.format(datatype,me)) remarks.append('adding {}{} "{}"...'.format(datatype if dim > 1 else 'scalar',
'' if dimension > 0 or dim == 1 else '[{}]'.format(dim),
me))
active[datatype].append(me) active[datatype].append(me)
if remarks != []: damask.util.croak(remarks) if remarks != []: damask.util.croak(remarks)

View File

@ -55,7 +55,7 @@ parser.set_defaults(data = [],
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
if not options.vtk: parser.error('No VTK file specified.') if not options.vtk: parser.error('no VTK file specified.')
if not os.path.exists(options.vtk): parser.error('VTK file does not exist.') if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
if os.path.splitext(options.vtk)[1] == '.vtr': if os.path.splitext(options.vtk)[1] == '.vtr':
@ -69,7 +69,7 @@ elif os.path.splitext(options.vtk)[1] == '.vtk':
reader.Update() reader.Update()
rGrid = reader.GetRectilinearGridOutput() rGrid = reader.GetRectilinearGridOutput()
else: else:
parser.error('Unsupported VTK file type extension.') parser.error('unsupported VTK file type extension.')
Npoints = rGrid.GetNumberOfPoints() Npoints = rGrid.GetNumberOfPoints()
Ncells = rGrid.GetNumberOfCells() Ncells = rGrid.GetNumberOfCells()
@ -96,16 +96,19 @@ for name in filenames:
VTKarray = {} VTKarray = {}
active = defaultdict(list) active = defaultdict(list)
for datatype,dimension,label in [['data',99,options.data], for datatype,dimension,label in [['data',0,options.data],
['tensor',9,options.tensor], ['tensor',9,options.tensor],
['color' ,3,options.color], ['color' ,3,options.color],
]: ]:
for i,dim in enumerate(table.label_dimension(label)): for i,dim in enumerate(table.label_dimension(label)):
me = label[i] me = label[i]
if dim == -1: remarks.append('{} "{}" not found...'.format(datatype,me)) if dim == -1: remarks.append('{} "{}" not found...'.format(datatype,me))
elif dim > dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension)) elif dimension > 0 \
and dim != dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension))
else: else:
remarks.append('adding {} "{}"...'.format(datatype,me)) remarks.append('adding {}{} "{}"...'.format(datatype if dim > 1 else 'scalar',
'' if dimension > 0 or dim == 1 else '[{}]'.format(dim),
me))
active[datatype].append(me) active[datatype].append(me)
if remarks != []: damask.util.croak(remarks) if remarks != []: damask.util.croak(remarks)
@ -141,7 +144,7 @@ for name in filenames:
if len(table.data) == Npoints: mode = 'point' if len(table.data) == Npoints: mode = 'point'
elif len(table.data) == Ncells: mode = 'cell' elif len(table.data) == Ncells: mode = 'cell'
else: else:
damask.util.croak('Data count is incompatible with grid...') damask.util.croak('data count is incompatible with grid...')
continue continue
damask.util.croak('{} mode...'.format(mode)) damask.util.croak('{} mode...'.format(mode))

View File

@ -37,12 +37,9 @@ def findClosestSeed(fargs):
return np.argmin(dist) # seed point closest to point return np.argmin(dist) # seed point closest to point
def laguerreTessellation(undeformed, coords, weights, grains, nonperiodic = False, cpus = 2): def laguerreTessellation(undeformed, coords, weights, grains, periodic = True, cpus = 2):
copies = \ copies = \
np.array([
[ 0, 0, 0 ],
]).astype(float) if nonperiodic else \
np.array([ np.array([
[ -1,-1,-1 ], [ -1,-1,-1 ],
[ 0,-1,-1 ], [ 0,-1,-1 ],
@ -71,7 +68,10 @@ def laguerreTessellation(undeformed, coords, weights, grains, nonperiodic = Fals
[ -1, 1, 1 ], [ -1, 1, 1 ],
[ 0, 1, 1 ], [ 0, 1, 1 ],
[ 1, 1, 1 ], [ 1, 1, 1 ],
]).astype(float)*info['size'] ]).astype(float)*info['size'] if periodic else \
np.array([
[ 0, 0, 0 ],
]).astype(float)
repeatweights = np.tile(weights,len(copies)).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) repeatweights = np.tile(weights,len(copies)).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
for i,vec in enumerate(copies): # periodic copies of seed points ... for i,vec in enumerate(copies): # periodic copies of seed points ...
@ -121,8 +121,8 @@ group.add_option('--cpus',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
help = 'number of parallel processes to use for Laguerre tessellation [%default]') help = 'number of parallel processes to use for Laguerre tessellation [%default]')
group.add_option('--nonperiodic', group.add_option('--nonperiodic',
dest = 'nonperiodic', dest = 'periodic',
action = 'store_true', action = 'store_false',
help = 'nonperiodic tessellation') help = 'nonperiodic tessellation')
parser.add_option_group(group) parser.add_option_group(group)
@ -144,6 +144,10 @@ group.add_option('-o',
dest = 'origin', dest = 'origin',
type = 'float', nargs = 3, metavar=' '.join(['float']*3), type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'origin of grid') help = 'origin of grid')
group.add_option('--nonnormalized',
dest = 'normalized',
action = 'store_false',
help = 'seed coordinates are not normalized to a unit cube')
parser.add_option_group(group) parser.add_option_group(group)
@ -206,7 +210,8 @@ parser.set_defaults(pos = 'pos',
phase = 1, phase = 1,
cpus = 2, cpus = 2,
laguerre = False, laguerre = False,
nonperiodic = False, periodic = True,
normalized = True,
config = True, config = True,
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -248,7 +253,7 @@ for name in filenames:
for i in range(3): for i in range(3):
if info['size'][i] <= 0.0: # any invalid size? if info['size'][i] <= 0.0: # any invalid size?
info['size'][i] = float(info['grid'][i])/max(info['grid']) # normalize to grid info['size'][i] = float(info['grid'][i])/max(info['grid']) # normalize to grid
remarks.append('rescaling size {} to {}...'.format({0:'x',1:'y',2:'z'}[i],info['size'][i])) remarks.append('rescaling size {} to {}...'.format(['x','y','z'][i],info['size'][i]))
if table.label_dimension(options.pos) != 3: if table.label_dimension(options.pos) != 3:
errors.append('seed positions "{}" have dimension {}.'.format(options.pos, errors.append('seed positions "{}" have dimension {}.'.format(options.pos,
@ -256,6 +261,7 @@ for name in filenames:
else: else:
labels += [options.pos] labels += [options.pos]
if not options.normalized: remarks.append('using real-space seed coordinates...')
if not hasEulers: remarks.append('missing seed orientations...') if not hasEulers: remarks.append('missing seed orientations...')
else: labels += [options.eulers] else: labels += [options.eulers]
if not hasGrains: remarks.append('missing seed microstructure indices...') if not hasGrains: remarks.append('missing seed microstructure indices...')
@ -272,7 +278,8 @@ for name in filenames:
# ------------------------------------------ read seeds --------------------------------------- # ------------------------------------------ read seeds ---------------------------------------
table.data_readArray(labels) table.data_readArray(labels)
coords = table.data[:,table.label_indexrange(options.pos)] * info['size'] coords = table.data[:,table.label_indexrange(options.pos)] * info['size'] if options.normalized \
else table.data[:,table.label_indexrange(options.pos)] - info['origin']
eulers = table.data[:,table.label_indexrange(options.eulers)] if hasEulers \ eulers = table.data[:,table.label_indexrange(options.eulers)] if hasEulers \
else np.zeros(3*len(coords)) else np.zeros(3*len(coords))
grains = table.data[:,table.label_indexrange(options.microstructure)].astype('i') if hasGrains \ grains = table.data[:,table.label_indexrange(options.microstructure)].astype('i') if hasGrains \
@ -291,7 +298,7 @@ for name in filenames:
damask.util.croak('tessellating...') damask.util.croak('tessellating...')
grid = np.vstack(meshgrid2(x, y, z)).reshape(3,-1).T grid = np.vstack(meshgrid2(x, y, z)).reshape(3,-1).T
indices = laguerreTessellation(grid, coords, weights, grains, options.nonperiodic, options.cpus) indices = laguerreTessellation(grid, coords, weights, grains, options.periodic, options.cpus)
# --- write header ------------------------------------------------------------------------ # --- write header ------------------------------------------------------------------------

View File

@ -31,7 +31,7 @@ parser.add_option('-b',
help = 'blacklist of grain IDs') help = 'blacklist of grain IDs')
parser.add_option('-p', parser.add_option('-p',
'--pos', '--seedposition', '--pos', '--seedposition',
dest = 'position', dest = 'pos',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'label of coordinates [%default]') help = 'label of coordinates [%default]')

View File

@ -6,7 +6,6 @@
!> @brief Mathematical library, including random number generation and tensor represenations !> @brief Mathematical library, including random number generation and tensor represenations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module math module math
use, intrinsic :: iso_c_binding
use prec, only: & use prec, only: &
pReal, & pReal, &
pInt pInt
@ -161,13 +160,10 @@ module math
math_rotate_forward3333, & math_rotate_forward3333, &
math_limit math_limit
private :: & private :: &
math_partition, &
halton, & halton, &
halton_memory, & halton_memory, &
halton_ndim_set, & halton_ndim_set, &
halton_seed_set, & halton_seed_set
i_to_halton, &
prime
contains contains
@ -289,58 +285,53 @@ recursive subroutine math_qsort(a, istart, iend)
integer(pInt) :: ipivot integer(pInt) :: ipivot
if (istart < iend) then if (istart < iend) then
ipivot = math_partition(a,istart, iend) ipivot = qsort_partition(a,istart, iend)
call math_qsort(a, istart, ipivot-1_pInt) call math_qsort(a, istart, ipivot-1_pInt)
call math_qsort(a, ipivot+1_pInt, iend) call math_qsort(a, ipivot+1_pInt, iend)
endif endif
end subroutine math_qsort
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
contains
!-------------------------------------------------------------------------------------------------
!> @brief Partitioning required for quicksort !> @brief Partitioning required for quicksort
!-------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
integer(pInt) function math_partition(a, istart, iend) integer(pInt) function qsort_partition(a, istart, iend)
implicit none implicit none
integer(pInt), dimension(:,:), intent(inout) :: a integer(pInt), dimension(:,:), intent(inout) :: a
integer(pInt), intent(in) :: istart,iend integer(pInt), intent(in) :: istart,iend
integer(pInt) :: d,i,j,k,x,tmp integer(pInt) :: i,j,k,tmp
d = int(size(a,1_pInt), pInt) ! number of linked data
! set the starting and ending points, and the pivot point
i = istart
j = iend
x = a(1,istart)
do do
! find the first element on the right side less than or equal to the pivot point ! find the first element on the right side less than or equal to the pivot point
do j = j, istart, -1_pInt do j = iend, istart, -1_pInt
if (a(1,j) <= x) exit if (a(1,j) <= a(1,istart)) exit
enddo enddo
! find the first element on the left side greater than the pivot point ! find the first element on the left side greater than the pivot point
do i = i, iend do i = istart, iend
if (a(1,i) > x) exit if (a(1,i) > a(1,istart)) exit
enddo enddo
if (i < j) then ! if the indexes do not cross, exchange values if (i < j) then ! if the indexes do not cross, exchange values
do k = 1_pInt,d do k = 1_pInt, int(size(a,1_pInt), pInt)
tmp = a(k,i) tmp = a(k,i)
a(k,i) = a(k,j) a(k,i) = a(k,j)
a(k,j) = tmp a(k,j) = tmp
enddo enddo
else ! if they do cross, exchange left value with pivot and return with the partition index else ! if they do cross, exchange left value with pivot and return with the partition index
do k = 1_pInt,d do k = 1_pInt, int(size(a,1_pInt), pInt)
tmp = a(k,istart) tmp = a(k,istart)
a(k,istart) = a(k,j) a(k,istart) = a(k,j)
a(k,j) = tmp a(k,j) = tmp
enddo enddo
math_partition = j qsort_partition = j
return return
endif endif
enddo enddo
end function math_partition end function qsort_partition
end subroutine math_qsort
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2189,6 +2180,52 @@ subroutine halton(ndim, r)
value_halton(1) = 1_pInt value_halton(1) = 1_pInt
call halton_memory ('INC', 'SEED', 1_pInt, value_halton) call halton_memory ('INC', 'SEED', 1_pInt, value_halton)
!--------------------------------------------------------------------------------------------------
contains
!-------------------------------------------------------------------------------------------------
!> @brief computes an element of a Halton sequence.
!> @details Only the absolute value of SEED is considered. SEED = 0 is allowed, and returns R = 0.
!> @details Halton Bases should be distinct prime numbers. This routine only checks that each base
!> @details is greater than 1.
!> @details Reference:
!> @details J.H. Halton: On the efficiency of certain quasi-random sequences of points in evaluating
!> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960.
!> @author John Burkardt
!-------------------------------------------------------------------------------------------------
subroutine i_to_halton (seed, base, ndim, r)
use IO, only: &
IO_error
implicit none
integer(pInt), intent(in) :: &
ndim, & !< dimension of the sequence
seed !< index of the desired element
integer(pInt), intent(in), dimension(ndim) :: base !< Halton bases
real(pReal), intent(out), dimension(ndim) :: r !< the SEED-th element of the Halton sequence for the given bases
real(pReal), dimension(ndim) :: base_inv
integer(pInt), dimension(ndim) :: &
digit, &
seed2
seed2 = abs(seed)
r = 0.0_pReal
if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=405_pInt)
base_inv(1:ndim) = 1.0_pReal / real (base(1:ndim), pReal)
do while ( any ( seed2(1:ndim) /= 0_pInt) )
digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim))
r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), pReal) * base_inv(1:ndim)
base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), pReal)
seed2(1:ndim) = seed2(1:ndim) / base(1:ndim)
enddo
end subroutine i_to_halton
end subroutine halton end subroutine halton
@ -2205,6 +2242,8 @@ end subroutine halton
!> @author John Burkardt !> @author John Burkardt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine halton_memory (action_halton, name_halton, ndim, value_halton) subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
use IO, only: &
IO_lc
implicit none implicit none
character(len = *), intent(in) :: & character(len = *), intent(in) :: &
@ -2214,8 +2253,8 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
integer(pInt), allocatable, save, dimension(:) :: base integer(pInt), allocatable, save, dimension(:) :: base
logical, save :: first_call = .true. logical, save :: first_call = .true.
integer(pInt), intent(in) :: ndim !< dimension of the quantity integer(pInt), intent(in) :: ndim !< dimension of the quantity
integer(pInt):: i
integer(pInt), save :: ndim_save = 0_pInt, seed = 1_pInt integer(pInt), save :: ndim_save = 0_pInt, seed = 1_pInt
integer(pInt) :: i
if (first_call) then if (first_call) then
ndim_save = 1_pInt ndim_save = 1_pInt
@ -2226,61 +2265,238 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Set ! Set
if(action_halton(1:1) == 'S' .or. action_halton(1:1) == 's') then actionHalton: if(IO_lc(action_halton(1:1)) == 's') then
if(name_halton(1:1) == 'B' .or. name_halton(1:1) == 'b') then
if(ndim_save /= ndim) then
deallocate(base)
ndim_save = ndim
allocate(base(ndim_save))
endif
base(1:ndim) = value_halton(1:ndim)
elseif(name_halton(1:1) == 'N' .or. name_halton(1:1) == 'n') then
nameSet: if(IO_lc(name_halton(1:1)) == 'b') then
if(ndim_save /= ndim) ndim_save = ndim
base = value_halton(1:ndim)
elseif(IO_lc(name_halton(1:1)) == 'n') then nameSet
if(ndim_save /= value_halton(1)) then if(ndim_save /= value_halton(1)) then
deallocate(base)
ndim_save = value_halton(1) ndim_save = value_halton(1)
allocate(base(ndim_save)) base = [(prime(i),i=1_pInt,ndim_save)]
do i = 1_pInt, ndim_save
base(i) = prime (i)
enddo
else else
ndim_save = value_halton(1) ndim_save = value_halton(1)
endif endif
elseif(name_halton(1:1) == 'S' .or. name_halton(1:1) == 's') then elseif(IO_lc(name_halton(1:1)) == 's') then nameSet
seed = value_halton(1) seed = value_halton(1)
endif endif nameSet
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Get ! Get
elseif(action_halton(1:1) == 'G' .or. action_halton(1:1) == 'g') then elseif(IO_lc(action_halton(1:1)) == 'g') then actionHalton
if(name_halton(1:1) == 'B' .or. name_halton(1:1) == 'b') then nameGet: if(IO_lc(name_halton(1:1)) == 'b') then
if(ndim /= ndim_save) then if(ndim /= ndim_save) then
deallocate(base)
ndim_save = ndim ndim_save = ndim
allocate(base(ndim_save)) base = [(prime(i),i=1_pInt,ndim_save)]
do i = 1_pInt, ndim_save
base(i) = prime(i)
enddo
endif endif
value_halton(1:ndim_save) = base(1:ndim_save) value_halton(1:ndim_save) = base(1:ndim_save)
elseif(name_halton(1:1) == 'N' .or. name_halton(1:1) == 'n') then elseif(IO_lc(name_halton(1:1)) == 'n') then nameGet
value_halton(1) = ndim_save value_halton(1) = ndim_save
elseif(name_halton(1:1) == 'S' .or. name_halton(1:1) == 's') then elseif(IO_lc(name_halton(1:1)) == 's') then nameGet
value_halton(1) = seed value_halton(1) = seed
endif endif nameGet
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Increment ! Increment
elseif(action_halton(1:1) == 'I' .or. action_halton(1:1) == 'i') then elseif(IO_lc(action_halton(1:1)) == 'i') then actionHalton
if(name_halton(1:1) == 'S' .or. name_halton(1:1) == 's') then if(IO_lc(name_halton(1:1)) == 's') seed = seed + value_halton(1)
seed = seed + value_halton(1) endif actionHalton
end if
!--------------------------------------------------------------------------------------------------
contains
!--------------------------------------------------------------------------------------------------
!> @brief returns any of the first 1500 prime numbers.
!> @details n = 0 is legal, returning PRIME = 1.
!> @details Reference:
!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions,
!> @details US Department of Commerce, 1964, pages 870-873.
!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae,
!> @details 30th Edition, CRC Press, 1996, pages 95-98.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
integer(pInt) function prime(n)
use IO, only: &
IO_error
implicit none
integer(pInt), intent(in) :: n !< index of the desired prime number
integer(pInt), dimension(0:1500), parameter :: &
npvec = int([&
1, &
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, &
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, &
73, 79, 83, 89, 97, 101, 103, 107, 109, 113, &
127, 131, 137, 139, 149, 151, 157, 163, 167, 173, &
179, 181, 191, 193, 197, 199, 211, 223, 227, 229, &
233, 239, 241, 251, 257, 263, 269, 271, 277, 281, &
283, 293, 307, 311, 313, 317, 331, 337, 347, 349, &
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, &
419, 421, 431, 433, 439, 443, 449, 457, 461, 463, &
467, 479, 487, 491, 499, 503, 509, 521, 523, 541, &
! 101:200
547, 557, 563, 569, 571, 577, 587, 593, 599, 601, &
607, 613, 617, 619, 631, 641, 643, 647, 653, 659, &
661, 673, 677, 683, 691, 701, 709, 719, 727, 733, &
739, 743, 751, 757, 761, 769, 773, 787, 797, 809, &
811, 821, 823, 827, 829, 839, 853, 857, 859, 863, &
877, 881, 883, 887, 907, 911, 919, 929, 937, 941, &
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, &
1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, &
1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, &
1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, &
! 201:300
1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, &
1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, &
1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, &
1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, &
1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, &
1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, &
1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, &
1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, &
1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, &
1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, &
! 301:400
1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, &
2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, &
2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, &
2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, &
2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, &
2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, &
2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, &
2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, &
2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, &
2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, &
! 401:500
2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, &
2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, &
2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, &
3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, &
3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, &
3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, &
3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, &
3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, &
3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, &
3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, &
! 501:600
3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, &
3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, &
3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, &
3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, &
3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, &
4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, &
4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, &
4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, &
4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, &
4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, &
! 601:700
4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, &
4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, &
4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, &
4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, &
4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, &
4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, &
4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, &
5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, &
5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, &
5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, &
! 701:800
5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, &
5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, &
5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, &
5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, &
5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, &
5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, &
5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, &
5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, &
5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, &
6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, &
! 801:900
6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, &
6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, &
6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, &
6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, &
6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, &
6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, &
6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, &
6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, &
6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, &
6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, &
! 901:1000
7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, &
7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, &
7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, &
7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, &
7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, &
7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, &
7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, &
7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, &
7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, &
7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, &
! 1001:1100
7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, &
8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, &
8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, &
8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, &
8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, &
8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, &
8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, &
8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, &
8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, &
8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831, &
! 1101:1200
8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, &
8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, &
9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, &
9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, &
9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, &
9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, &
9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, &
9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, &
9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, &
9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, &
! 1201:1300
9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, &
9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, &
9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973, 10007, &
10009, 10037, 10039, 10061, 10067, 10069, 10079, 10091, 10093, 10099, &
10103, 10111, 10133, 10139, 10141, 10151, 10159, 10163, 10169, 10177, &
10181, 10193, 10211, 10223, 10243, 10247, 10253, 10259, 10267, 10271, &
10273, 10289, 10301, 10303, 10313, 10321, 10331, 10333, 10337, 10343, &
10357, 10369, 10391, 10399, 10427, 10429, 10433, 10453, 10457, 10459, &
10463, 10477, 10487, 10499, 10501, 10513, 10529, 10531, 10559, 10567, &
10589, 10597, 10601, 10607, 10613, 10627, 10631, 10639, 10651, 10657, &
! 1301:1400
10663, 10667, 10687, 10691, 10709, 10711, 10723, 10729, 10733, 10739, &
10753, 10771, 10781, 10789, 10799, 10831, 10837, 10847, 10853, 10859, &
10861, 10867, 10883, 10889, 10891, 10903, 10909, 19037, 10939, 10949, &
10957, 10973, 10979, 10987, 10993, 11003, 11027, 11047, 11057, 11059, &
11069, 11071, 11083, 11087, 11093, 11113, 11117, 11119, 11131, 11149, &
11159, 11161, 11171, 11173, 11177, 11197, 11213, 11239, 11243, 11251, &
11257, 11261, 11273, 11279, 11287, 11299, 11311, 11317, 11321, 11329, &
11351, 11353, 11369, 11383, 11393, 11399, 11411, 11423, 11437, 11443, &
11447, 11467, 11471, 11483, 11489, 11491, 11497, 11503, 11519, 11527, &
11549, 11551, 11579, 11587, 11593, 11597, 11617, 11621, 11633, 11657, &
! 1401:1500
11677, 11681, 11689, 11699, 11701, 11717, 11719, 11731, 11743, 11777, &
11779, 11783, 11789, 11801, 11807, 11813, 11821, 11827, 11831, 11833, &
11839, 11863, 11867, 11887, 11897, 11903, 11909, 11923, 11927, 11933, &
11939, 11941, 11953, 11959, 11969, 11971, 11981, 11987, 12007, 12011, &
12037, 12041, 12043, 12049, 12071, 12073, 12097, 12101, 12107, 12109, &
12113, 12119, 12143, 12149, 12157, 12161, 12163, 12197, 12203, 12211, &
12227, 12239, 12241, 12251, 12253, 12263, 12269, 12277, 12281, 12289, &
12301, 12323, 12329, 12343, 12347, 12373, 12377, 12379, 12391, 12401, &
12409, 12413, 12421, 12433, 12437, 12451, 12457, 12473, 12479, 12487, &
12491, 12497, 12503, 12511, 12517, 12527, 12539, 12541, 12547, 12553],pInt)
if (n < size(npvec)) then
prime = npvec(n)
else
call IO_error(error_ID=406_pInt)
end if end if
end function prime
end subroutine halton_memory end subroutine halton_memory
@ -2325,252 +2541,6 @@ subroutine halton_seed_set(seed)
end subroutine halton_seed_set end subroutine halton_seed_set
!--------------------------------------------------------------------------------------------------
!> @brief computes an element of a Halton sequence.
!> @details Only the absolute value of SEED is considered. SEED = 0 is allowed, and returns R = 0.
!> @details Halton Bases should be distinct prime numbers. This routine only checks that each base
!> @details is greater than 1.
!> @details Reference:
!> @details J.H. Halton: On the efficiency of certain quasi-random sequences of points in evaluating
!> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
subroutine i_to_halton (seed, base, ndim, r)
use IO, only: &
IO_error
implicit none
integer(pInt), intent(in) :: ndim !< dimension of the sequence
integer(pInt), intent(in), dimension(ndim) :: base !< Halton bases
real(pReal), dimension(ndim) :: base_inv
integer(pInt), dimension(ndim) :: digit
real(pReal), dimension(ndim), intent(out) ::r !< the SEED-th element of the Halton sequence for the given bases
integer(pInt) , intent(in):: seed !< index of the desired element
integer(pInt), dimension(ndim) :: seed2
seed2(1:ndim) = abs(seed)
r(1:ndim) = 0.0_pReal
if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=405_pInt)
base_inv(1:ndim) = 1.0_pReal / real (base(1:ndim), pReal)
do while ( any ( seed2(1:ndim) /= 0_pInt) )
digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim))
r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), pReal) * base_inv(1:ndim)
base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), pReal)
seed2(1:ndim) = seed2(1:ndim) / base(1:ndim)
enddo
end subroutine i_to_halton
!--------------------------------------------------------------------------------------------------
!> @brief returns any of the first 1500 prime numbers.
!> @details n <= 0 returns 1500, the index of the largest prime (12553) available.
!> @details n = 0 is legal, returning PRIME = 1.
!> @details Reference:
!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions,
!> @details US Department of Commerce, 1964, pages 870-873.
!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae,
!> @details 30th Edition, CRC Press, 1996, pages 95-98.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
integer(pInt) function prime(n)
use IO, only: &
IO_error
implicit none
integer(pInt), intent(in) :: n !< index of the desired prime number
integer(pInt), parameter :: PRIME_MAX = 1500_pInt
integer(pInt), save :: icall = 0_pInt
integer(pInt), save, dimension(PRIME_MAX) :: npvec
if (icall == 0_pInt) then
icall = 1_pInt
npvec = [&
2_pInt, 3_pInt, 5_pInt, 7_pInt, 11_pInt, 13_pInt, 17_pInt, 19_pInt, 23_pInt, 29_pInt, &
31_pInt, 37_pInt, 41_pInt, 43_pInt, 47_pInt, 53_pInt, 59_pInt, 61_pInt, 67_pInt, 71_pInt, &
73_pInt, 79_pInt, 83_pInt, 89_pInt, 97_pInt, 101_pInt, 103_pInt, 107_pInt, 109_pInt, 113_pInt, &
127_pInt, 131_pInt, 137_pInt, 139_pInt, 149_pInt, 151_pInt, 157_pInt, 163_pInt, 167_pInt, 173_pInt, &
179_pInt, 181_pInt, 191_pInt, 193_pInt, 197_pInt, 199_pInt, 211_pInt, 223_pInt, 227_pInt, 229_pInt, &
233_pInt, 239_pInt, 241_pInt, 251_pInt, 257_pInt, 263_pInt, 269_pInt, 271_pInt, 277_pInt, 281_pInt, &
283_pInt, 293_pInt, 307_pInt, 311_pInt, 313_pInt, 317_pInt, 331_pInt, 337_pInt, 347_pInt, 349_pInt, &
353_pInt, 359_pInt, 367_pInt, 373_pInt, 379_pInt, 383_pInt, 389_pInt, 397_pInt, 401_pInt, 409_pInt, &
419_pInt, 421_pInt, 431_pInt, 433_pInt, 439_pInt, 443_pInt, 449_pInt, 457_pInt, 461_pInt, 463_pInt, &
467_pInt, 479_pInt, 487_pInt, 491_pInt, 499_pInt, 503_pInt, 509_pInt, 521_pInt, 523_pInt, 541_pInt, &
! 101:200
547_pInt, 557_pInt, 563_pInt, 569_pInt, 571_pInt, 577_pInt, 587_pInt, 593_pInt, 599_pInt, 601_pInt, &
607_pInt, 613_pInt, 617_pInt, 619_pInt, 631_pInt, 641_pInt, 643_pInt, 647_pInt, 653_pInt, 659_pInt, &
661_pInt, 673_pInt, 677_pInt, 683_pInt, 691_pInt, 701_pInt, 709_pInt, 719_pInt, 727_pInt, 733_pInt, &
739_pInt, 743_pInt, 751_pInt, 757_pInt, 761_pInt, 769_pInt, 773_pInt, 787_pInt, 797_pInt, 809_pInt, &
811_pInt, 821_pInt, 823_pInt, 827_pInt, 829_pInt, 839_pInt, 853_pInt, 857_pInt, 859_pInt, 863_pInt, &
877_pInt, 881_pInt, 883_pInt, 887_pInt, 907_pInt, 911_pInt, 919_pInt, 929_pInt, 937_pInt, 941_pInt, &
947_pInt, 953_pInt, 967_pInt, 971_pInt, 977_pInt, 983_pInt, 991_pInt, 997_pInt, 1009_pInt, 1013_pInt, &
1019_pInt, 1021_pInt, 1031_pInt, 1033_pInt, 1039_pInt, 1049_pInt, 1051_pInt, 1061_pInt, 1063_pInt, 1069_pInt, &
1087_pInt, 1091_pInt, 1093_pInt, 1097_pInt, 1103_pInt, 1109_pInt, 1117_pInt, 1123_pInt, 1129_pInt, 1151_pInt, &
1153_pInt, 1163_pInt, 1171_pInt, 1181_pInt, 1187_pInt, 1193_pInt, 1201_pInt, 1213_pInt, 1217_pInt, 1223_pInt, &
! 201:300
1229_pInt, 1231_pInt, 1237_pInt, 1249_pInt, 1259_pInt, 1277_pInt, 1279_pInt, 1283_pInt, 1289_pInt, 1291_pInt, &
1297_pInt, 1301_pInt, 1303_pInt, 1307_pInt, 1319_pInt, 1321_pInt, 1327_pInt, 1361_pInt, 1367_pInt, 1373_pInt, &
1381_pInt, 1399_pInt, 1409_pInt, 1423_pInt, 1427_pInt, 1429_pInt, 1433_pInt, 1439_pInt, 1447_pInt, 1451_pInt, &
1453_pInt, 1459_pInt, 1471_pInt, 1481_pInt, 1483_pInt, 1487_pInt, 1489_pInt, 1493_pInt, 1499_pInt, 1511_pInt, &
1523_pInt, 1531_pInt, 1543_pInt, 1549_pInt, 1553_pInt, 1559_pInt, 1567_pInt, 1571_pInt, 1579_pInt, 1583_pInt, &
1597_pInt, 1601_pInt, 1607_pInt, 1609_pInt, 1613_pInt, 1619_pInt, 1621_pInt, 1627_pInt, 1637_pInt, 1657_pInt, &
1663_pInt, 1667_pInt, 1669_pInt, 1693_pInt, 1697_pInt, 1699_pInt, 1709_pInt, 1721_pInt, 1723_pInt, 1733_pInt, &
1741_pInt, 1747_pInt, 1753_pInt, 1759_pInt, 1777_pInt, 1783_pInt, 1787_pInt, 1789_pInt, 1801_pInt, 1811_pInt, &
1823_pInt, 1831_pInt, 1847_pInt, 1861_pInt, 1867_pInt, 1871_pInt, 1873_pInt, 1877_pInt, 1879_pInt, 1889_pInt, &
1901_pInt, 1907_pInt, 1913_pInt, 1931_pInt, 1933_pInt, 1949_pInt, 1951_pInt, 1973_pInt, 1979_pInt, 1987_pInt, &
! 301:400
1993_pInt, 1997_pInt, 1999_pInt, 2003_pInt, 2011_pInt, 2017_pInt, 2027_pInt, 2029_pInt, 2039_pInt, 2053_pInt, &
2063_pInt, 2069_pInt, 2081_pInt, 2083_pInt, 2087_pInt, 2089_pInt, 2099_pInt, 2111_pInt, 2113_pInt, 2129_pInt, &
2131_pInt, 2137_pInt, 2141_pInt, 2143_pInt, 2153_pInt, 2161_pInt, 2179_pInt, 2203_pInt, 2207_pInt, 2213_pInt, &
2221_pInt, 2237_pInt, 2239_pInt, 2243_pInt, 2251_pInt, 2267_pInt, 2269_pInt, 2273_pInt, 2281_pInt, 2287_pInt, &
2293_pInt, 2297_pInt, 2309_pInt, 2311_pInt, 2333_pInt, 2339_pInt, 2341_pInt, 2347_pInt, 2351_pInt, 2357_pInt, &
2371_pInt, 2377_pInt, 2381_pInt, 2383_pInt, 2389_pInt, 2393_pInt, 2399_pInt, 2411_pInt, 2417_pInt, 2423_pInt, &
2437_pInt, 2441_pInt, 2447_pInt, 2459_pInt, 2467_pInt, 2473_pInt, 2477_pInt, 2503_pInt, 2521_pInt, 2531_pInt, &
2539_pInt, 2543_pInt, 2549_pInt, 2551_pInt, 2557_pInt, 2579_pInt, 2591_pInt, 2593_pInt, 2609_pInt, 2617_pInt, &
2621_pInt, 2633_pInt, 2647_pInt, 2657_pInt, 2659_pInt, 2663_pInt, 2671_pInt, 2677_pInt, 2683_pInt, 2687_pInt, &
2689_pInt, 2693_pInt, 2699_pInt, 2707_pInt, 2711_pInt, 2713_pInt, 2719_pInt, 2729_pInt, 2731_pInt, 2741_pInt, &
! 401:500
2749_pInt, 2753_pInt, 2767_pInt, 2777_pInt, 2789_pInt, 2791_pInt, 2797_pInt, 2801_pInt, 2803_pInt, 2819_pInt, &
2833_pInt, 2837_pInt, 2843_pInt, 2851_pInt, 2857_pInt, 2861_pInt, 2879_pInt, 2887_pInt, 2897_pInt, 2903_pInt, &
2909_pInt, 2917_pInt, 2927_pInt, 2939_pInt, 2953_pInt, 2957_pInt, 2963_pInt, 2969_pInt, 2971_pInt, 2999_pInt, &
3001_pInt, 3011_pInt, 3019_pInt, 3023_pInt, 3037_pInt, 3041_pInt, 3049_pInt, 3061_pInt, 3067_pInt, 3079_pInt, &
3083_pInt, 3089_pInt, 3109_pInt, 3119_pInt, 3121_pInt, 3137_pInt, 3163_pInt, 3167_pInt, 3169_pInt, 3181_pInt, &
3187_pInt, 3191_pInt, 3203_pInt, 3209_pInt, 3217_pInt, 3221_pInt, 3229_pInt, 3251_pInt, 3253_pInt, 3257_pInt, &
3259_pInt, 3271_pInt, 3299_pInt, 3301_pInt, 3307_pInt, 3313_pInt, 3319_pInt, 3323_pInt, 3329_pInt, 3331_pInt, &
3343_pInt, 3347_pInt, 3359_pInt, 3361_pInt, 3371_pInt, 3373_pInt, 3389_pInt, 3391_pInt, 3407_pInt, 3413_pInt, &
3433_pInt, 3449_pInt, 3457_pInt, 3461_pInt, 3463_pInt, 3467_pInt, 3469_pInt, 3491_pInt, 3499_pInt, 3511_pInt, &
3517_pInt, 3527_pInt, 3529_pInt, 3533_pInt, 3539_pInt, 3541_pInt, 3547_pInt, 3557_pInt, 3559_pInt, 3571_pInt, &
! 501:600
3581_pInt, 3583_pInt, 3593_pInt, 3607_pInt, 3613_pInt, 3617_pInt, 3623_pInt, 3631_pInt, 3637_pInt, 3643_pInt, &
3659_pInt, 3671_pInt, 3673_pInt, 3677_pInt, 3691_pInt, 3697_pInt, 3701_pInt, 3709_pInt, 3719_pInt, 3727_pInt, &
3733_pInt, 3739_pInt, 3761_pInt, 3767_pInt, 3769_pInt, 3779_pInt, 3793_pInt, 3797_pInt, 3803_pInt, 3821_pInt, &
3823_pInt, 3833_pInt, 3847_pInt, 3851_pInt, 3853_pInt, 3863_pInt, 3877_pInt, 3881_pInt, 3889_pInt, 3907_pInt, &
3911_pInt, 3917_pInt, 3919_pInt, 3923_pInt, 3929_pInt, 3931_pInt, 3943_pInt, 3947_pInt, 3967_pInt, 3989_pInt, &
4001_pInt, 4003_pInt, 4007_pInt, 4013_pInt, 4019_pInt, 4021_pInt, 4027_pInt, 4049_pInt, 4051_pInt, 4057_pInt, &
4073_pInt, 4079_pInt, 4091_pInt, 4093_pInt, 4099_pInt, 4111_pInt, 4127_pInt, 4129_pInt, 4133_pInt, 4139_pInt, &
4153_pInt, 4157_pInt, 4159_pInt, 4177_pInt, 4201_pInt, 4211_pInt, 4217_pInt, 4219_pInt, 4229_pInt, 4231_pInt, &
4241_pInt, 4243_pInt, 4253_pInt, 4259_pInt, 4261_pInt, 4271_pInt, 4273_pInt, 4283_pInt, 4289_pInt, 4297_pInt, &
4327_pInt, 4337_pInt, 4339_pInt, 4349_pInt, 4357_pInt, 4363_pInt, 4373_pInt, 4391_pInt, 4397_pInt, 4409_pInt, &
! 601:700
4421_pInt, 4423_pInt, 4441_pInt, 4447_pInt, 4451_pInt, 4457_pInt, 4463_pInt, 4481_pInt, 4483_pInt, 4493_pInt, &
4507_pInt, 4513_pInt, 4517_pInt, 4519_pInt, 4523_pInt, 4547_pInt, 4549_pInt, 4561_pInt, 4567_pInt, 4583_pInt, &
4591_pInt, 4597_pInt, 4603_pInt, 4621_pInt, 4637_pInt, 4639_pInt, 4643_pInt, 4649_pInt, 4651_pInt, 4657_pInt, &
4663_pInt, 4673_pInt, 4679_pInt, 4691_pInt, 4703_pInt, 4721_pInt, 4723_pInt, 4729_pInt, 4733_pInt, 4751_pInt, &
4759_pInt, 4783_pInt, 4787_pInt, 4789_pInt, 4793_pInt, 4799_pInt, 4801_pInt, 4813_pInt, 4817_pInt, 4831_pInt, &
4861_pInt, 4871_pInt, 4877_pInt, 4889_pInt, 4903_pInt, 4909_pInt, 4919_pInt, 4931_pInt, 4933_pInt, 4937_pInt, &
4943_pInt, 4951_pInt, 4957_pInt, 4967_pInt, 4969_pInt, 4973_pInt, 4987_pInt, 4993_pInt, 4999_pInt, 5003_pInt, &
5009_pInt, 5011_pInt, 5021_pInt, 5023_pInt, 5039_pInt, 5051_pInt, 5059_pInt, 5077_pInt, 5081_pInt, 5087_pInt, &
5099_pInt, 5101_pInt, 5107_pInt, 5113_pInt, 5119_pInt, 5147_pInt, 5153_pInt, 5167_pInt, 5171_pInt, 5179_pInt, &
5189_pInt, 5197_pInt, 5209_pInt, 5227_pInt, 5231_pInt, 5233_pInt, 5237_pInt, 5261_pInt, 5273_pInt, 5279_pInt, &
! 701:800
5281_pInt, 5297_pInt, 5303_pInt, 5309_pInt, 5323_pInt, 5333_pInt, 5347_pInt, 5351_pInt, 5381_pInt, 5387_pInt, &
5393_pInt, 5399_pInt, 5407_pInt, 5413_pInt, 5417_pInt, 5419_pInt, 5431_pInt, 5437_pInt, 5441_pInt, 5443_pInt, &
5449_pInt, 5471_pInt, 5477_pInt, 5479_pInt, 5483_pInt, 5501_pInt, 5503_pInt, 5507_pInt, 5519_pInt, 5521_pInt, &
5527_pInt, 5531_pInt, 5557_pInt, 5563_pInt, 5569_pInt, 5573_pInt, 5581_pInt, 5591_pInt, 5623_pInt, 5639_pInt, &
5641_pInt, 5647_pInt, 5651_pInt, 5653_pInt, 5657_pInt, 5659_pInt, 5669_pInt, 5683_pInt, 5689_pInt, 5693_pInt, &
5701_pInt, 5711_pInt, 5717_pInt, 5737_pInt, 5741_pInt, 5743_pInt, 5749_pInt, 5779_pInt, 5783_pInt, 5791_pInt, &
5801_pInt, 5807_pInt, 5813_pInt, 5821_pInt, 5827_pInt, 5839_pInt, 5843_pInt, 5849_pInt, 5851_pInt, 5857_pInt, &
5861_pInt, 5867_pInt, 5869_pInt, 5879_pInt, 5881_pInt, 5897_pInt, 5903_pInt, 5923_pInt, 5927_pInt, 5939_pInt, &
5953_pInt, 5981_pInt, 5987_pInt, 6007_pInt, 6011_pInt, 6029_pInt, 6037_pInt, 6043_pInt, 6047_pInt, 6053_pInt, &
6067_pInt, 6073_pInt, 6079_pInt, 6089_pInt, 6091_pInt, 6101_pInt, 6113_pInt, 6121_pInt, 6131_pInt, 6133_pInt, &
! 801:900
6143_pInt, 6151_pInt, 6163_pInt, 6173_pInt, 6197_pInt, 6199_pInt, 6203_pInt, 6211_pInt, 6217_pInt, 6221_pInt, &
6229_pInt, 6247_pInt, 6257_pInt, 6263_pInt, 6269_pInt, 6271_pInt, 6277_pInt, 6287_pInt, 6299_pInt, 6301_pInt, &
6311_pInt, 6317_pInt, 6323_pInt, 6329_pInt, 6337_pInt, 6343_pInt, 6353_pInt, 6359_pInt, 6361_pInt, 6367_pInt, &
6373_pInt, 6379_pInt, 6389_pInt, 6397_pInt, 6421_pInt, 6427_pInt, 6449_pInt, 6451_pInt, 6469_pInt, 6473_pInt, &
6481_pInt, 6491_pInt, 6521_pInt, 6529_pInt, 6547_pInt, 6551_pInt, 6553_pInt, 6563_pInt, 6569_pInt, 6571_pInt, &
6577_pInt, 6581_pInt, 6599_pInt, 6607_pInt, 6619_pInt, 6637_pInt, 6653_pInt, 6659_pInt, 6661_pInt, 6673_pInt, &
6679_pInt, 6689_pInt, 6691_pInt, 6701_pInt, 6703_pInt, 6709_pInt, 6719_pInt, 6733_pInt, 6737_pInt, 6761_pInt, &
6763_pInt, 6779_pInt, 6781_pInt, 6791_pInt, 6793_pInt, 6803_pInt, 6823_pInt, 6827_pInt, 6829_pInt, 6833_pInt, &
6841_pInt, 6857_pInt, 6863_pInt, 6869_pInt, 6871_pInt, 6883_pInt, 6899_pInt, 6907_pInt, 6911_pInt, 6917_pInt, &
6947_pInt, 6949_pInt, 6959_pInt, 6961_pInt, 6967_pInt, 6971_pInt, 6977_pInt, 6983_pInt, 6991_pInt, 6997_pInt, &
! 901:1000
7001_pInt, 7013_pInt, 7019_pInt, 7027_pInt, 7039_pInt, 7043_pInt, 7057_pInt, 7069_pInt, 7079_pInt, 7103_pInt, &
7109_pInt, 7121_pInt, 7127_pInt, 7129_pInt, 7151_pInt, 7159_pInt, 7177_pInt, 7187_pInt, 7193_pInt, 7207_pInt, &
7211_pInt, 7213_pInt, 7219_pInt, 7229_pInt, 7237_pInt, 7243_pInt, 7247_pInt, 7253_pInt, 7283_pInt, 7297_pInt, &
7307_pInt, 7309_pInt, 7321_pInt, 7331_pInt, 7333_pInt, 7349_pInt, 7351_pInt, 7369_pInt, 7393_pInt, 7411_pInt, &
7417_pInt, 7433_pInt, 7451_pInt, 7457_pInt, 7459_pInt, 7477_pInt, 7481_pInt, 7487_pInt, 7489_pInt, 7499_pInt, &
7507_pInt, 7517_pInt, 7523_pInt, 7529_pInt, 7537_pInt, 7541_pInt, 7547_pInt, 7549_pInt, 7559_pInt, 7561_pInt, &
7573_pInt, 7577_pInt, 7583_pInt, 7589_pInt, 7591_pInt, 7603_pInt, 7607_pInt, 7621_pInt, 7639_pInt, 7643_pInt, &
7649_pInt, 7669_pInt, 7673_pInt, 7681_pInt, 7687_pInt, 7691_pInt, 7699_pInt, 7703_pInt, 7717_pInt, 7723_pInt, &
7727_pInt, 7741_pInt, 7753_pInt, 7757_pInt, 7759_pInt, 7789_pInt, 7793_pInt, 7817_pInt, 7823_pInt, 7829_pInt, &
7841_pInt, 7853_pInt, 7867_pInt, 7873_pInt, 7877_pInt, 7879_pInt, 7883_pInt, 7901_pInt, 7907_pInt, 7919_pInt, &
! 1001:1100
7927_pInt, 7933_pInt, 7937_pInt, 7949_pInt, 7951_pInt, 7963_pInt, 7993_pInt, 8009_pInt, 8011_pInt, 8017_pInt, &
8039_pInt, 8053_pInt, 8059_pInt, 8069_pInt, 8081_pInt, 8087_pInt, 8089_pInt, 8093_pInt, 8101_pInt, 8111_pInt, &
8117_pInt, 8123_pInt, 8147_pInt, 8161_pInt, 8167_pInt, 8171_pInt, 8179_pInt, 8191_pInt, 8209_pInt, 8219_pInt, &
8221_pInt, 8231_pInt, 8233_pInt, 8237_pInt, 8243_pInt, 8263_pInt, 8269_pInt, 8273_pInt, 8287_pInt, 8291_pInt, &
8293_pInt, 8297_pInt, 8311_pInt, 8317_pInt, 8329_pInt, 8353_pInt, 8363_pInt, 8369_pInt, 8377_pInt, 8387_pInt, &
8389_pInt, 8419_pInt, 8423_pInt, 8429_pInt, 8431_pInt, 8443_pInt, 8447_pInt, 8461_pInt, 8467_pInt, 8501_pInt, &
8513_pInt, 8521_pInt, 8527_pInt, 8537_pInt, 8539_pInt, 8543_pInt, 8563_pInt, 8573_pInt, 8581_pInt, 8597_pInt, &
8599_pInt, 8609_pInt, 8623_pInt, 8627_pInt, 8629_pInt, 8641_pInt, 8647_pInt, 8663_pInt, 8669_pInt, 8677_pInt, &
8681_pInt, 8689_pInt, 8693_pInt, 8699_pInt, 8707_pInt, 8713_pInt, 8719_pInt, 8731_pInt, 8737_pInt, 8741_pInt, &
8747_pInt, 8753_pInt, 8761_pInt, 8779_pInt, 8783_pInt, 8803_pInt, 8807_pInt, 8819_pInt, 8821_pInt, 8831_pInt, &
! 1101:1200
8837_pInt, 8839_pInt, 8849_pInt, 8861_pInt, 8863_pInt, 8867_pInt, 8887_pInt, 8893_pInt, 8923_pInt, 8929_pInt, &
8933_pInt, 8941_pInt, 8951_pInt, 8963_pInt, 8969_pInt, 8971_pInt, 8999_pInt, 9001_pInt, 9007_pInt, 9011_pInt, &
9013_pInt, 9029_pInt, 9041_pInt, 9043_pInt, 9049_pInt, 9059_pInt, 9067_pInt, 9091_pInt, 9103_pInt, 9109_pInt, &
9127_pInt, 9133_pInt, 9137_pInt, 9151_pInt, 9157_pInt, 9161_pInt, 9173_pInt, 9181_pInt, 9187_pInt, 9199_pInt, &
9203_pInt, 9209_pInt, 9221_pInt, 9227_pInt, 9239_pInt, 9241_pInt, 9257_pInt, 9277_pInt, 9281_pInt, 9283_pInt, &
9293_pInt, 9311_pInt, 9319_pInt, 9323_pInt, 9337_pInt, 9341_pInt, 9343_pInt, 9349_pInt, 9371_pInt, 9377_pInt, &
9391_pInt, 9397_pInt, 9403_pInt, 9413_pInt, 9419_pInt, 9421_pInt, 9431_pInt, 9433_pInt, 9437_pInt, 9439_pInt, &
9461_pInt, 9463_pInt, 9467_pInt, 9473_pInt, 9479_pInt, 9491_pInt, 9497_pInt, 9511_pInt, 9521_pInt, 9533_pInt, &
9539_pInt, 9547_pInt, 9551_pInt, 9587_pInt, 9601_pInt, 9613_pInt, 9619_pInt, 9623_pInt, 9629_pInt, 9631_pInt, &
9643_pInt, 9649_pInt, 9661_pInt, 9677_pInt, 9679_pInt, 9689_pInt, 9697_pInt, 9719_pInt, 9721_pInt, 9733_pInt, &
! 1201:1300
9739_pInt, 9743_pInt, 9749_pInt, 9767_pInt, 9769_pInt, 9781_pInt, 9787_pInt, 9791_pInt, 9803_pInt, 9811_pInt, &
9817_pInt, 9829_pInt, 9833_pInt, 9839_pInt, 9851_pInt, 9857_pInt, 9859_pInt, 9871_pInt, 9883_pInt, 9887_pInt, &
9901_pInt, 9907_pInt, 9923_pInt, 9929_pInt, 9931_pInt, 9941_pInt, 9949_pInt, 9967_pInt, 9973_pInt,10007_pInt, &
10009_pInt,10037_pInt,10039_pInt,10061_pInt,10067_pInt,10069_pInt,10079_pInt,10091_pInt,10093_pInt,10099_pInt, &
10103_pInt,10111_pInt,10133_pInt,10139_pInt,10141_pInt,10151_pInt,10159_pInt,10163_pInt,10169_pInt,10177_pInt, &
10181_pInt,10193_pInt,10211_pInt,10223_pInt,10243_pInt,10247_pInt,10253_pInt,10259_pInt,10267_pInt,10271_pInt, &
10273_pInt,10289_pInt,10301_pInt,10303_pInt,10313_pInt,10321_pInt,10331_pInt,10333_pInt,10337_pInt,10343_pInt, &
10357_pInt,10369_pInt,10391_pInt,10399_pInt,10427_pInt,10429_pInt,10433_pInt,10453_pInt,10457_pInt,10459_pInt, &
10463_pInt,10477_pInt,10487_pInt,10499_pInt,10501_pInt,10513_pInt,10529_pInt,10531_pInt,10559_pInt,10567_pInt, &
10589_pInt,10597_pInt,10601_pInt,10607_pInt,10613_pInt,10627_pInt,10631_pInt,10639_pInt,10651_pInt,10657_pInt, &
! 1301:1400
10663_pInt,10667_pInt,10687_pInt,10691_pInt,10709_pInt,10711_pInt,10723_pInt,10729_pInt,10733_pInt,10739_pInt, &
10753_pInt,10771_pInt,10781_pInt,10789_pInt,10799_pInt,10831_pInt,10837_pInt,10847_pInt,10853_pInt,10859_pInt, &
10861_pInt,10867_pInt,10883_pInt,10889_pInt,10891_pInt,10903_pInt,10909_pInt,19037_pInt,10939_pInt,10949_pInt, &
10957_pInt,10973_pInt,10979_pInt,10987_pInt,10993_pInt,11003_pInt,11027_pInt,11047_pInt,11057_pInt,11059_pInt, &
11069_pInt,11071_pInt,11083_pInt,11087_pInt,11093_pInt,11113_pInt,11117_pInt,11119_pInt,11131_pInt,11149_pInt, &
11159_pInt,11161_pInt,11171_pInt,11173_pInt,11177_pInt,11197_pInt,11213_pInt,11239_pInt,11243_pInt,11251_pInt, &
11257_pInt,11261_pInt,11273_pInt,11279_pInt,11287_pInt,11299_pInt,11311_pInt,11317_pInt,11321_pInt,11329_pInt, &
11351_pInt,11353_pInt,11369_pInt,11383_pInt,11393_pInt,11399_pInt,11411_pInt,11423_pInt,11437_pInt,11443_pInt, &
11447_pInt,11467_pInt,11471_pInt,11483_pInt,11489_pInt,11491_pInt,11497_pInt,11503_pInt,11519_pInt,11527_pInt, &
11549_pInt,11551_pInt,11579_pInt,11587_pInt,11593_pInt,11597_pInt,11617_pInt,11621_pInt,11633_pInt,11657_pInt, &
! 1401:1500
11677_pInt,11681_pInt,11689_pInt,11699_pInt,11701_pInt,11717_pInt,11719_pInt,11731_pInt,11743_pInt,11777_pInt, &
11779_pInt,11783_pInt,11789_pInt,11801_pInt,11807_pInt,11813_pInt,11821_pInt,11827_pInt,11831_pInt,11833_pInt, &
11839_pInt,11863_pInt,11867_pInt,11887_pInt,11897_pInt,11903_pInt,11909_pInt,11923_pInt,11927_pInt,11933_pInt, &
11939_pInt,11941_pInt,11953_pInt,11959_pInt,11969_pInt,11971_pInt,11981_pInt,11987_pInt,12007_pInt,12011_pInt, &
12037_pInt,12041_pInt,12043_pInt,12049_pInt,12071_pInt,12073_pInt,12097_pInt,12101_pInt,12107_pInt,12109_pInt, &
12113_pInt,12119_pInt,12143_pInt,12149_pInt,12157_pInt,12161_pInt,12163_pInt,12197_pInt,12203_pInt,12211_pInt, &
12227_pInt,12239_pInt,12241_pInt,12251_pInt,12253_pInt,12263_pInt,12269_pInt,12277_pInt,12281_pInt,12289_pInt, &
12301_pInt,12323_pInt,12329_pInt,12343_pInt,12347_pInt,12373_pInt,12377_pInt,12379_pInt,12391_pInt,12401_pInt, &
12409_pInt,12413_pInt,12421_pInt,12433_pInt,12437_pInt,12451_pInt,12457_pInt,12473_pInt,12479_pInt,12487_pInt, &
12491_pInt,12497_pInt,12503_pInt,12511_pInt,12517_pInt,12527_pInt,12539_pInt,12541_pInt,12547_pInt,12553_pInt]
endif
if(n < 0_pInt) then
prime = PRIME_MAX
else if (n == 0_pInt) then
prime = 1_pInt
else if (n <= PRIME_MAX) then
prime = npvec(n)
else
prime = -1_pInt
call IO_error(error_ID=406_pInt)
end if
end function prime
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief factorial !> @brief factorial
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------