Kedarnath J Balakrishnan

Nur A Touba
email: touba@ece.utexas.edu

Reconfigurable Linear Decompressors Using Symbolic Gaussian Elimination

HAL is

Introduction

As circuits become increasingly complex with each generation, the amount of test data required to test each chip is also becoming humongous. Hence test data storage requirements on the tester and test data bandwidth requirements between the tester and chip are growing rapidly [START_REF] Khoche | I/O bandwidth bottleneck for test: Is it real?[END_REF]. Test data compression techniques provide a way to reduce these requirements thereby allowing less expensive testers to be used. Test time can also be simultaneously reduced with most test data compression techniques. The output response is relatively easy to compress since lossy compression techniques can be employed, e.g., using a multiple input signature register (MISR). However, compressing test stimuli is much more difficult because lossless compression techniques must be used. A lot of research has been done on lossless compression techniques for test vectors since reducing test data volume has become such an important problem.

Test vector compression schemes that use only linear operations to decompress the test vectors are called linear decompression schemes. Linear decompression techniques exploit the unspecified (don't care) bit positions in scan test cubes (i.e., deterministic scan test vectors where the unassigned bit positions are left as don't cares) to achieve large amounts of compression. A number of different techniques for designing linear decompressors have been proposed in * Kedarnath J. Balakrishnan is now with NEC Labs. in Princeton, NJ the literature. These include both techniques based on linear feedback shift register (LFSR) reseeding and combinational linear expansion circuits consisting of XOR gates. The original idea of using an LFSR as a linear decompressor and solving for test cubes using linear algebra was described in [START_REF] Konemann | LFSR-coded test patterns for scan designs[END_REF].

A decompressor based on XOR tree network was described in [START_REF] Bayraktaroglu | Test volume and application time reduction through scan chain concealment[END_REF]. This is an example of a continuous-flow decompressor. Continuous-flow linear decompressors are those that receive data from the tester in a continuous-flow manner i.e. every cycle. These operate very efficiently since they can be directly connected to the tester and they simply receive the data as fast as the tester can transfer it. This simulates scan chain concealment since for the tester it mimics the normal behavior of scan chains though the number of scan chains visible to the tester is much less than that in the design. Combinational continuous-flow linear decompressors are described in [START_REF] Hamzaoglu | Reducing test application time for full scan embedded cores[END_REF], [START_REF] Hsu | A case study on the implementation of the illinois scan architecture[END_REF], [START_REF] Bayraktaroglu | Test volume and application time reduction through scan chain concealment[END_REF][START_REF] Bayraktaroglu | Decompression hardware determination for test volume and time reduction through unified test pattern compaction and compression[END_REF], [START_REF] Mitra | Xmax: X-tolerant architectures for maximal test compression[END_REF], and [START_REF] Krishna | Adjustable width linear combinational scan vector decompression[END_REF]. Continuous-flow techniques that have sequential elements like LFSRs i.e. sequential continuous-flow linear decompressors are described in [START_REF] Jas | Virtual scan chains: A means for reducing scan length in cores[END_REF], [START_REF] Konemann | A smartBIST variant with guaranteed encoding[END_REF], [START_REF] Rajski | Embedded deterministic test for low cost manufacturing test[END_REF], [START_REF] Rao | Test application time and volume compression through seed overlapping[END_REF], and [START_REF] Krishna | 3-stage variable length continuous-flow scan vector decompression scheme[END_REF]. Most of the commercial tools for compressing test vectors are also based on linear decompressors. TestKompress from Mentor Graphics [START_REF] Rajski | Embedded deterministic test for low cost manufacturing test[END_REF], SmartBIST from IBM/Cadence [START_REF] Konemann | A smartBIST variant with guaranteed encoding[END_REF], and DBIST [START_REF] Chandramouli | How to implement deterministic logic built-in self-test (BIST). Complier: A Monthly Magazine for Technologists Worldwide[END_REF] In order to be able to compress a test set, the output space of the linear decompressor must contain all the test cubes in the test set. The output space of a linear decompressor is a linear subspace spanned by a Boolean matrix, A c×n . Each row in A corresponds to a scan cell and each column corresponds to a free-variable in the input sequence. To determine whether a particular test cube is contained in the out-

O 1 O 2 O 3 O 4 O 5 O 6 O 7 O 8 O 9 x 0 x 1 x 2 x 3 x 4
1 0 1 0 0 1 0 0 1 0 1 0 0 0 1 1 1 0 0 0 0 1 1 0 0 0 1 0 1 0 0 1 0 0 1 0 0 1 1 0 0 0 0 1 1 1 put space and the corresponding input sequence to generate it can be quickly done through Gaussian elimination.

Figure 2. XOR Network and Matrix for Figure

Figure 1 shows a combinational linear decompressor that receives b bits from the tester and expands it to n scan chains. The Boolean matrix for this decompressor can be constructed simply from the XOR network. Each row corresponding to a scan chain will have the inputs that are xored together to get the scan chain value as 1 and all the others as 0. An example of a XOR network with b = 5 and n = 9 and the corresponding matrix is shown in Figure 2. The columns of the matrix correspond to inputs with the first column representing x 0 and the last column representing x 4 .

Output O 1 is the XOR of inputs x 0 and x 2 and hence in the first row, the values corresponding to these two inputs are 1 while the rest are 0. Obtaining the Boolean matrix A by symbolic simulation of the linear decompressor is described in detail in [START_REF] Krishna | Reducing test data volume using LFSR reseeding with seed compression[END_REF].

When a LFSR is used as the decompressor, it has been proved that if the number of free variables used to generate a test cube is 20 more than the number of specified bits in test cube, then the probability of the test cube not being solvable is less than one in a million [START_REF] Konemann | LFSR-coded test patterns for scan designs[END_REF]. However, for a given test set, the number of free variables can be further reduced provided the corresponding equations for each test cube are solvable. Hence, traditionally linear decompressors are designed on the basis of the worst-case scenario. In LFSR reseeding the size of LFSR is usually proportional to s max -the maximum specified bits in any test cube of the test set.

Currently test compression methods are deployed in either of two formats. The first one is to decide on a compression scheme and the decompressor beforehand, during the circuit design stage itself and then generate the test patterns for the design based on that scheme. In this method, the decompressor circuit is also integrated into the design flow. Hence, fault coverage during test pattern generation may be affected due to the limitations of the decompressor. Also, last minute design changes may need test pattern regeneration and there may be some patterns that are not compressible using the given scheme. Techniques that can "reconfigure" the decompressor to take into account these changes will be very helpful. The other method is to design the decompressor after the circuit has been finalized and the test patterns available. Legacy circuits or intellectual property (IP) cores in a system-on-a-chip (SoC) environment are good examples of this. Each core in a SoC may require test compression scheme. Using a single decompressor that can be "reconfigured" for several cores instead of having several decompressors can result in a lot of area savings.

The idea of adding a reconfigurable part for test data reduction has been proposed earlier. [START_REF] Li | Deterministic BIST based on a reconfigurable interconnection network[END_REF] employs a reconfigurable interconnect network (RIN) within BIST to embed deterministic test patterns. The RIN has been proposed as an alternative to traditional phase shifters consisting of XOR gates. In [START_REF] Tang | On reducing test data volume and test application time for multiple scan chain designs[END_REF], the original Illinois scan architecture [START_REF] Hamzaoglu | Reducing test application time for full scan embedded cores[END_REF] is altered by using a reconfigurable switch to control the connection between external pins and the scan chains. In [START_REF] Krishna | Adjustable width linear combinational scan vector decompression[END_REF], control bits are used to determine the number of scan chains that are fed by the combinational network.

This paper describes a methodology for "reconfiguring" any linear decompressor. It can be implemented in conjunction with any decompressor that has a combinational linear network and significantly improve the compression obtained. The rest of the paper is organized as follows. Section 2 describes the idea of a reconfigurable linear decompressor. In Section 3, we present a symbolic Gaussian elimination technique that is used to design the reconfigurable linear decompressor and Section 4 shows how it can be used to improve the compression. Experimental results are presented in Section 5 and Section 6 has conclusions.

Reconfigurable Linear Decompressor

The idea of reconfiguring a linear decompressor is to modify the output space of the linear decompressor to ensure that a given test cube can be generated using the linear decompressor. The linear decompression network can be modified using the configuration bits. For each configuration, the output space of the linear decompressor will be different. Hence this increases the chances of a test cube being generated using the decompressor. The configuration bits can either be stored on the tester and transferred to the decompressor with the test cubes or stored on-chip using a ROM. The number of different configurations required depends on the test set. In the case when each test cube in a test set requires a different configuration, these need to be explicitly stored. In the case when all the test cubes can be generated using the same configuration, then the decompression network can be finalized for that configuration and no extra storage is required. This implies that the decompressor is simply redesigned.

There are several possible ways to modify a linear decompressor such that it can be reconfigured. A simple mod-ification would be to add multiplexers to each output of the network with the select input coming from the configuration bits. The other inputs could be some other output of the original decompressor. For example, consider a combinational decompressor using XOR gates shown in Figure 2. Instead of directly sending the outputs to the scan chains, there is another stage with a multiplexer before each scan chain as illustrated in Figure 3. These multiplexers will select which output is connected to which scan chain based on the configuration bits.

O 1 O 2 O 3 O 4 O 5 O 6 O 7 O 8 O 9 x 0 x 1 x 2 x 3 x 4
Configuration Bits

Figure 3. Making a Decompressor Reconfigurable

For a test cube to be solvable using the linear decompressor, a solution to the system of linear equations A x = b must exist where x is an assignment of values to the free-variables that are inputs to the decompressor when generating the test cube, and b is the value of each bit in the test cube. There is no need to solve the linear equations for the unspecified bits in the test cube, and hence only the linear equations (rows) corresponding to the specified bits in b need to be considered. Gaussian elimination [START_REF] Cullen | Linear Algebra with Applications[END_REF] can be used to perform rows operations that transform a set of columns into an identity matrix. The elements that make the identity matrix are called the pivots. An example of a system of linear equations for a test cube t 1 and the corresponding system after Gaussian elimination is shown in Fig. 4. The rows after Gaussian elimination can be classified as either pivoted rows or linearly dependent rows. The pivoted rows have pivots while the linearly dependent rows are all 0. For the example in Fig. 4, the first two rows and the last row are pivoted while the third row is linearly dependent. If all rows are pivoted, then a solution to the system of linear equations exists, and hence the test cube can be decompressed using the linear decompressor. If some of the rows are linearly dependent, then a solution only exists if all of the corresponding values in z (the vector b after Gaussian elimination) are equal to 0 for the linearly dependent rows. If there is a linearly dependent row whose corresponding value in z is equal to 1, then no solution exists. In Fig. 4, the third row is linearly dependent but the corresponding value in z is 1, and thus there is no solution.

For the example in Fig. 4, let the specified values in b correspond to scan cells c 1 through c 4 . If the decompressor is reconfigured such that the equations corresponding to scan cell c 2 and scan cell c 4 were exchanged, the system of linear equations would become solvable. This is shown in Fig. 5. Here, the last row is now linearly dependent and the corresponding value in z is 0. This is an example of how reconfiguration can be used to make a test cube solvable.

Symbolic Gaussian Elimination

The previous section described how reconfiguration of a linear decompressor can increase the chances of a test cube being solvable. In this section we describe a systematic procedure to do reconfiguration using Symbolic Gaussian Elimination. The key idea is to form the matrix A in terms of the configuration bits (i.e. each entry in the matrix is a function of the configuration bits) and then find an assignment of the configuration bits that makes the system of linear equations A x = b solvable.

In traditional Gaussian elimination, elementary row operations are used to reduce the coefficient matrix (the matrix A) to a set of pivoted and linearly dependent rows. This is done by going through each column of the matrix and choosing a non-zero element as the pivot. All other rows with a 1 in this column are xor-ed with this row as part of elementary row operation. The resultant matrix will either have pivoted rows or linearly dependent rows.

We extend this technique to co-efficient matrices that have functions as elements. The functions are Boolean functions on a given set of variables. This is illustrated in Fig. 6 where each element of the matrix F is a function of the variables {a, b, c}. In this matrix, each f i j is a function such that f i j : {0, 1} 3 → {0, 1}

In the trivial case when each function is a constant 1 or a constant 0, this matrix will degenerate to a boolean matrix with entries 0 or 1. The variables {a, b, c} correspond to the configuration bits of the linear decompressor. Solving for a system of linear equations such as F x = y, with F as the matrix shown in Fig. 6 and y is a vector implies that we are looking for a solution for at least one combination of a, b & c, or in other words, for at least one configuration. Given such a system of linear equations, the algorithm for symbolic Gaussian elimination is given below. First, each element in the y vector is converted to a function in terms of the variables in F. If the element in y is 1 then the function is identically equal to 1 and if the element is 0, the function is identically equal to 0. Then, the algorithm proceeds column wise choosing pivots in each column. Suppose f 11 is the first pivot. Row operations are performed next for each row and this pivot. In the row operations, every element in a row is XOR-ed with the result of the AND of its element in the pivot column and the corresponding element in the pivot row. This is illustrated in Fig. 7 that shows the matrix after row operations for pivot f 11 . The idea is to take into account both cases when f 11 = 0 and f 11 = 1. This operation is repeated for every pivot. The matrix after doing all the row operations for each pivot will look like Fig. 8 where f i j are again functions of {a, b, c}. For ease of explanation, assume that the pivot for each column i is given by element f ii . The next step would be to ensure that each pivot has atleast one minterm for which the function equates to one. If the pivot function does not have a single minterm that equates to one, then that row is equivalent to a linearly dependent row in the normal Gaussian eliminated matrix. Hence the corresponding element in the y matrix should be zero. This condition for a pivot f ii can be written in mathematical form as

f ii + ȳi = 1
Since this condition must be valid for each pivot, the overall condition can be written in a product of sum form as

( f 11 + ȳ1 )( f 22 + ȳ2 ) •••( f nn + ȳn ) = 1
If the above condition is satisfied, there exists a solution to the system of equations. The number of minterms of {a, b, c} for which the above condition is satisfied will indicate the number of different configurations possible.

Note that there may be more than one possible pivot for each column and hence the pivot is selected using a heuristic. For the first column, the element that has the maximum number of minterms is chosen as the pivot since it has the best chances of having a solution. For the next columns, the element that has the most number of minterms in common with the current pivots is chosen. This ensures that the algorithms proceeds in such a way that maximum number of solutions (configurations) are possible.

The algorithm described above can be implemented with very little overhead with respect to the basic Gaussian elimination method. Each function is stored in terms of its minterms. For example, if three configuration bits are used, then there are eight possible minterms. Each element in the matrix consists of eight values, one corresponding to each minterm. The row operations are performed on the corresponding minterms. The only additional step in this procedure is evaluating the final product of sums condition. This can be achieved by simply performing bitwise operations on the entries so that the corresponding minterms are evaluated together. The number of entries in the product of sums condition depends on the number of pivots and hence on the size of the matrix. The complexity of the evaluation step increases linearly with the size, since the bitwise operations need to be performed for each additional pivot. The number of minterms for each function depend on the number of configuration bits. The number of configuration bits is a design parameter that can be decided depending on the compression required and the maximum allowed running time of the algorithm.

Increasing Compression by Reconfiguration

Given a linear decompressor, the systematic procedure given in the previous section can be used to reconfigure it and increase the chances of compressing any test set. This reconfiguration can be done in two ways. The first one is to search for a single configuration by which all the cubes in the test set can be compressed. In this case, the configuration can be hardwired into the decompressor i.e. the decompressor is redesigned and no explicit configuration bits are required. The other method would be to have one configuration for each cube which is loaded into the decompressor every time a new test cube is loaded. The configuration bits for each test cube need to be stored explicitly.

Compression obtained using a linear decompression scheme can be improved using reconfiguration in several ways. One method would be to reduce the number of freevariables that the decompressor receives per test cube from the tester as much as possible while still keeping the test set compressible through reconfiguration. Any method can be used to design the initial decompressor. Then the number of free-variables that are input to the decompressor per test cube can be incrementally reduced and symbolic Gaussian elimination can be used to check whether it is possible to still solve for all the test cubes using reconfiguration. If so, then this process of incrementally reducing the number of free-variables and checking for a solution is repeated until a point is reached when no further reduction in the number of free-variables per test cube is possible while still being able to solve for all test cubes. The end result will be a linear decompressor that generates the exact same test set, but uses fewer tester channels thereby reducing tester storage and bandwidth requirements.

Another method is to keep the number of free-variables that the decompressor receives per test cube constant, but use reconfiguration to relax the constraints on ATPG (automatic test pattern generation) such that more specified bits per test cube can be generated. This will allow more static and dynamic compaction while still being able to solve for the test cubes. Some test compression methodologies (e.g., [START_REF] Bayraktaroglu | Test volume and application time reduction through scan chain concealment[END_REF][START_REF] Bayraktaroglu | Decompression hardware determination for test volume and time reduction through unified test pattern compaction and compression[END_REF], [START_REF] Rajski | Embedded deterministic test for low cost manufacturing test[END_REF]) involve fixing the decompressor design and then constraining the ATPG so that the resulting test cubes will be in the output space of the decompressor. The constraints on the ATPG reduce the amount of static and dynamic compaction that are performed and therefore can result in more test cubes and hence more test time. Reconfiguration of the linear decompressor can be used to allow more specified bits per test cube while still being able to solve for the test cube. This can be used to relax the constraints on the ATPG and thereby allow more static and dynamic compaction which will in turn reduce the total number of test cubes and hence result in a reduction of both test time and tester storage requirements.

Experimental Results

The effectiveness of the proposed method was evaluated by performing two sets of experiments. The first set of experiments consisting those described in the previous section were performed on randomly generated test cubes for large industrial-size scan architectures. Circuits were assumed to have either 512 or 1024 scan chains and the number of channels available from tester to be 32. The length of the scan chains varied from 24 to 128. The initial decompressor is an XOR network with 32 inputs and 512 or 1024 outputs. For reconfiguration, a four input multiplexer was added to the output of the decompressor network with the select bits coming from the configuration bits. All the experiments were performed with a total of eight configuration bits. The results are presented in Table 1. The first two columns show the number of scan chains and length of each scan chain. The columns under "without reconfiguration" show the maximum compression that can be obtained using the given decompressor. For each randomly generated test cube, the number of specified bits was incrementally increased until it could no longer be solved using the given decompressor. The maximum percentage of specified bits per test cube that could be solved by the given decompressor is shown in column 4. The corresponding encoding efficiency and compression ratios are given in columns 5 and 6. Encoding efficiency is defined as the ratio of number of specified bits in the test set to the number of compressed bits that need to be stored on the tester. Compression ratio is the ratio of the original bits in the test set to the compressed bits.

The columns under "Reducing channels" show the results for reduction in the number of tester channels using reconfiguration. The reduced number of tester channels along with the corresponding compression ratio are shown under "Red. Chan." and "Compr. Ratio". The reduced channels include one channel from the tester for the configuration bits to be loaded before each test cube. The column "% Impr." shows the improvement in compression ratio due to the reconfiguration. The average percentage improvement for all the different scan sizes is around 39.4 % which is a significant improvement for very little extra hardware overhead. The columns under "Increasing Spec. Bits" show the results for increasing the percentage of specified bits that can be handled by a given decompressor using reconfiguration. The tester channels are kept constant at 32 and the number of specified bits are increased as much as possible until it is no longer possible to solve for all the test cubes. The new percentage of specified bits and encoding efficiency are shown as well as the percentage improvement in the encoding efficiency. The experiments assume a single configuration for each test cube and the configuration bits are taken into account while calculating the compression results. In this experiment, the average percentage improvement for all the different scan sizes is around 36.5 %.

The other set of experiments were performed on 100 % stuck-at fault coverage test sets for the largest ISCAS '89 [START_REF] Brglez | Combinational profiles of sequential benchmark circuits[END_REF] benchmark circuits. Table 2 compares the compression results obtained using the proposed scheme with some of the combinational decompressor techniques proposed earlier. The number of tester channels required and the amount of test data that need to be stored on the tester for the XOR network [START_REF] Bayraktaroglu | Test volume and application time reduction through scan chain concealment[END_REF] and the adjustable width technique [START_REF] Krishna | Adjustable width linear combinational scan vector decompression[END_REF] are compared with those of the proposed scheme. Note that in [START_REF] Bayraktaroglu | Test volume and application time reduction through scan chain concealment[END_REF] the decompressor design is integrated into the test pattern generation, while in [START_REF] Krishna | Adjustable width linear combinational scan vector decompression[END_REF] compression is performed on already generated test patterns. The proposed scheme uses the same test patterns as [START_REF] Krishna | Adjustable width linear combinational scan vector decompression[END_REF]. As expected, the proposed scheme performs better than [START_REF] Bayraktaroglu | Test volume and application time reduction through scan chain concealment[END_REF] for all circuits, both in terms of the number of tester channels required and the amount of compressed data. In comparison with [START_REF] Krishna | Adjustable width linear combinational scan vector decompression[END_REF], the number of channels required are smaller or the same but the amount of compressed data is slightly higher. Note that the proposed scheme doesn't need to generate separate scan clocks for each group of scan chains as required in [START_REF] Krishna | Adjustable width linear combinational scan vector decompression[END_REF].

Conclusions

A symbolic Gaussian elimination method to solve a constrained boolean matrix equation is presented. An application of the above method to reconfigure any linear decompressor is proposed. The reconfiguration can be implemented with very little hardware overhead. Experimental results show that compression obtained using a linear decompressor can be significantly improved using reconfiguration.

Figure 1 .

 1 Figure 1. Combinational Linear Decompressor

Figure 4 . 1 1Figure 5 .

 415 Figure 4. System of Equations for Test Cube t 1

f 11 (Figure 6 .

 116 Figure 6. A Matrix with Functions as Elements

Figure 7. Matrix after Row Operations for Pivot f 11

	f 11	f 12		. . .	f 1m
	f 21 f 21 f 11	f 22 f 21 f 12	. . .	f 2m f 21 f 1m
	f 31 f 31 f 11	f 32 f 31 f 12	. . .	f 3m f 31 f 1m
	. . .	. . .		. . .	. . .
	f n1 f n1 f 11	f n2 f n1 f 12	. . .	f nm f n1 f 1m
	f ' 11	f ' 12	f ' 13	. . .	f ' 1m
	f ' 21	f ' 22	f ' 23	. . .	f ' 2m
	f ' 31	f ' 32	f ' 33	. . .	f ' 3m
	. . .	. . .	. . .	. . .	. . .
	f ' n1	f ' n2	f ' n3	. . .	f ' nm
	Figure 8.

Matrix after All Row Operations

Table 1 . Results for Randomly Generated Test Cubes

 1

	Num.	Scan		Without Reconfiguration		Reducing Channels	Increasing Spec. Bits
	Scan	Chain	Tester Percent. Encod. Compr.	Red.	Compr.	%	Incr. % Encod.	%
	Chains	Length	Chan.	Spec.	Effic.	Ratio	Chan.	Ratio	Impr.	Spec.	Effic.	Impr.
		24	32	3.61 %	0.578	16.0	23	22.3	39.4 % 4.69 %	0.727	25.8%
	512	32 64	32 32	3.42 % 2.83 %	0.547 0.453	16.0 16.0	23 22	22.3 23.3	39.4 % 4.30 % 45.5 % 4.00 %	0.667 0.621	21.9% 37.1%
		128	32	2.67 %	0.427	16.0	24	21.3	33.3 % 2.91 %	0.661	54.8 %
		24	32	1.66 %	0.531	32.0	24	42.7	33.3 % 2.34 %	0.727	36.9%
	1024	32 64	32 32	1.66 % 1.49 %	0.531 0.476	32.0 32.0	24 22	42.7 46.5	33.3 % 2.25 % 45.5% 2.18 %	0.697 0.677	31.3% 42.2%
		128	32	1.50 %	0.480	32.0	22	46.5	45.5%	2.20 %	0.680	41.7%

Table 2 . Results for ISCAS'89 Benchmarks

 2

		[1]		[13]	Prop. Scheme
	Circuit	Tester Chan.	Test Data	Tester Chan.	Test Data	Tester Chan.	Test Data
	s13207	24	25344	19	14307	15	14098
	s15850	32	22784	19	15067	18	18080
	s38417	32	89856	19	49001	19	54020
	s38584	24	38976	19	28994	19	31436

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE


Blondeau Da

Silva Stéphane

Benford or not Benford: a systematic but not always well-founded use of an elegant law in experimental fields

In this paper, we will see that the proportion of d as leading digit, d ∈ 1, 9 , in data (obtained thanks to the hereunder developed model) is more likely to follow a law whose probability distribution is determined by a specific upper bound, rather than Benford's Law. These probability distributions fluctuate around Benford's value as can often be observed in the literature in many naturally occurring collections of data (where the physical, biological or economical quantities considered are upper bounded). Knowing beforehand the value of the upper bound can be a way to find a better adjusted law than Benford's one.

Introduction

Benford's Law, also called Newcomb-Benford's Law, is noteworthy to say the least: according to it, the first digit d, d ∈ 1, 9 , of numbers in many naturally occurring collections of data does not follow a discrete uniform distribution, as might be thought, but a logarithmic distribution. Discovered by the astronomer Newcomb in 1881 ( [START_REF] Newcomb | Note on the frequency of use of the different digits in natural numbers[END_REF]), this law was definitively brought to light by the physicist Benford in 1938 ( [START_REF] Benford | The law of anomalous numbers[END_REF]). He proposed the following probability distribution: the probability for d to be the first digit of a number seems to be equal to log(d + 1) -log(d), i.e. log(1 + 1 d ). Benford tested it over data set from 20 different domains (surface areas of rivers, sizes of american populations, physical constants, molecular weights, entries from a mathematical handbook, numbers contained in an issue of Reader's Digest, the street addresses of the first persons listed in American Men of Science, death rates, etc.). Most of the empirical data, as physical data (Knuth in [START_REF] Knuth | The Art of Computer Programming 2[END_REF] or Burke and Kincanon in [START_REF] Burke | Benford's law and physical constants: the distribution of initial digits[END_REF]), economic and demographic data (Nigrini and Wood in [START_REF] Nigrini | Assessing the integrity of tabulated demographic data[END_REF]) or genome data (Friar et al. in [START_REF] Friar | Genome sizes and the Benford distribution[END_REF]), follow approximately Benford's Law. To such an extent that this law is used to detect possible frauds in lists of socio-economic data ( [START_REF] Varian | Benford's law (letters to the editor)[END_REF]) or in scientific publications ( [START_REF] Diekmann | Not the first digit! using benford's law to detect fraudulent scientific data[END_REF]).

First restricted to the experimental field, it is now established that this law holds for various mathematical sequences (see for example [START_REF] Bunimovich | One-dimensional dynamical systems and benford's law[END_REF]). In the situation, where the distribution of first digits is scale, unit or base invariant, this distribution is always given by Benford's Law ( [START_REF] Pinkham | On the distribution of first significant digits[END_REF] and [START_REF] Hill | Base-invariance implies benford's law[END_REF]). Selecting different samples in different populations, under certain constraints, leads also to construct a sequence that follows the Benford's Law ( [START_REF] Hill | A statistical derivation of the significant-digit law[END_REF]). Furthermore independant variables multiplication conducts to this law ( [START_REF] Boyle | An application of fourier series to the most significant digit problem[END_REF]). One might add that some sequences satisfy Benford's Law exactly (for example see [START_REF] Sarkar | An observation on the significant digits of binomial coefficients and factorials[END_REF], [START_REF] Washington | Benford's law for fibonacci and lucas numbers[END_REF] or [START_REF] Jolissaint | Loi de benford, relations de récurrence et suites équidistribuées[END_REF]).

We can note that there also exist distributions known to disobey Benford's Law ( [START_REF] Raimi | The first digit problem[END_REF] and [START_REF] Beer | Terminal digit preference: beware of benford's law[END_REF]). And even concerning empirical data sets, this law appears to be a good approximation of the reality, but no more than an approximation ( [START_REF] Gauvrit | Pourquoi la loi de benford n'est pas mystérieuse[END_REF]).

In the model we build in the article, the naturally occurring data will be considered as the realizations of independant random variables following the hereinafter constraints: (a) the data is strictly positive and is upper-bounded by an integer n, constraint which is often valid in data sets, the physical, biological and economical quantities being limited ; (b) each random variable is considered to follow a discrete uniform distribution whereby the first i strictly positive integers are equally likely to occur (i being uniformly randomly selected in 1, n ). This model relies on the fact that the random variables are not always the same. The article is divided into two parts. In the first one, we will accurately study the case where the leading digit is 1. In the second one, we will generalize our results to the eight last cases.

Through this article we will demonstrate that the predominance of 1 as first digit (followed by those of 2 and so on) is all but surprising, and that the observed fluctuations around the values of probability determined by Benford's Law are also predictible. The point is that, since 1938, Benford's Law probabilities became standard values that should exactly be followed by most of naturally occurring collections of data. However the reality is that the proportion of each d as leading digit, d ∈ 1, 9 , structurally fluctuates. There is not a single Benford's Law but numerous distinct laws that we will hereafter examine.

1 The chosen probability space

Notations

In order to determine the proportion of numbers whose leading digit is d ∈ 1, 9 , we will first build our probability space and further explain the model we choose.

Let i be a strictly positive integer. Let U {i} denote the discrete uniform distribution whereby the first i strictly positive integers are equally likely to be observed.

Let n be a strictly positive integer. Let us consider the random experiment E n of tossing two independent dice. The first one is a fair n-sided die showing n different numbers from 1 to n. The number i rolled on it defines the number of faces on the second die. It thus shows i different numbers from 1 to i.

Let us now define the probability space Ω n as follows: Ω n = {(i, j) : i ∈ 1, n and j ∈ 1, i }. Our probability measure is denoted by P.

Let us denote by L n the random variable from Ω n to 1, 9 that maps each element ω of Ω n to the leading digit of the second component of ω.

Why such a model?

Let us imagine a perfect consumer shopping in a perfectly structured supermarket: (a) in the i th (i being a strictly positive integer) section of this supermarket, the products prices range between 1 and i cents of the considered currency; (b) the prices in a section are uniformly distributed; (c) each section contains the same quantity of products and (d) the consumer randomly chooses his products in the whole store.

Under these constraining hypotheses, these perfect entities enable us to use our model. Note that conditions (c) and (d) gathered avoid us to conduct a double drawing every time: first the section then the product. In that respect, the sales receipt will verify the following results in terms of proportion of d as leading digit, d ∈ 1, 9 .

Among the different domains studied by Benford ([2]), some could be well adapted to our model: sizes of populations (sections here gathering all the populations having the same usable areas, the geographic constraints preventing the surface area to be broader; populations being not neccessary settled on the entire area, their sizes fluctuate) or street adresses for example (sections here gathering the adresses of a selected street; the lenght of the considered streets might be uniformly distributed to fit model criteria).

Hence the defined model is relevant when the studied data can be considered as realizations of a homogeneous and expanded range of random variables approximately following discrete uniform distributions.

Proportion of d

Through the below proposition, we will express the probability P(L n = d), for each n ∈ N * , i.e. the probability that the leading digit of our second throw in our random experiment is d.

Proposition 2.1. Let k denote the positive integer such that k = min{i ∈ N : d × 10 i > n}. If n < (d + 1) × 10 k-1 , the value of P(L n = d) is: 1 n k-2 l=0 ( (d+1)×10 l -1 b=d×10 l b - (9d-1)×10 l -8 9 b + d×10 l+1 -1 a=(d+1)×10 l 10 l+1 -1 9 a )+ n b=d×10 k-1 b - (9d-1)×10 k-1 -8 9 b .
Otherwise the value of P(L n = d) is:

1 n k-2 l=0 ( (d+1)×10 l -1 b=d×10 l b - (9d-1)×10 l -8 9 b + d×10 l+1 -1 a=(d+1)×10 l 10 l+1 -1 9 a ) + (d+1)×10 k-1 -1 b=d×10 k-1 b - (9d-1)×10 k-1 -1 9 b + n a=(d+1)×10 k-1 10 k -1 9
a .

Proof. Let us denote by D n the random variable from Ω n to 1, n that maps each element ω of Ω n to the first component of ω. It returns the number obtained on the first throw of the unbiased n-sided die. For each i ∈ 1, n , we have:

P(D n = i) = 1 n . (1)
According to the Law of total probability, we state:

P(L n = d) = n i=1 P(L n = d|D n = i) P(D n = i) . (2)
Thereupon two cases appear in determining the value of P(L n = d|D n = i), for i ∈ 1, n . Let us study the first case where the leading digit of i is d. Let k i be the positive integer such that k i = min{k ∈ N : d × 10 k > i} in both cases. Among the first d × 10 ki-1 -1 non-zero integers (all lower than i), the number of integers whose leading digit is

d is (if k i ≥ 2): ki-2 t=0 10 t = 1 × 1 -10 ki-1 1 -10 = 10 ki-1 -1 9 .
This equality still holds true for k i = 1. From d × 10 ki-1 to i, there exist i -d × 10 ki-1 + 1 additional integers whose leading digit is d. It may be inferred that:

P(L n = d|D n = i) = 1 i ( 10 ki-1 -1 9 + i -d × 10 ki-1 + 1) = i -(9d-1)10 k i -1 -8 9 i , (3)
the leading digit of i being here d.

In the second case, we consider the integers i whose leading digits are different from d. Among the first (d + 1) × 10 ki-1 -1 non-zero integers (i is greater than or equal to (d + 1) × 10 ki-1 ), the number of integers whose leading digit is d is: From 2 × 10 ki-1 to i, there exists no additional integers whose leading digit is d. It can be concluded that:

P(L n = d|D n = i) = 10 k i -1 9 i , (4)
the leading digit of i being here different from d.

Using equalities (1), ( 2), ( 3) and (4), we get our result.

For example, we get: b -8

Study of two subsequences

It is natural that we take a specific look at the values of n positioned just before a long sequence of numbers whose leading digit is d ; or conversely, at those positioned just before a long sequence of numbers whose leading digit is all but d.

To this end we will consider the sequence (P(L n = d)) n∈N * . In the interests of simplifying notation, we will denote by (P (d,n) ) n∈N * this sequence. Let us study two of its subsequences.

The first subsequence

The first one is the subsequence (P (d,φ d (n)) ) n∈N * where φ d is the function from N * to N that maps n to d × 10 n -1. We get the below result: Proof. Let n be a positive integer such that n ≥ 2. According to Proposition 2.1, we have:

P (d,φ d (n)) = P (d,d×10 n -1) = 1 d × 10 n -1 n-1 l=0 ( (d+1)×10 l -1 b=d×10 l b -(9d-1)×10 l -8 9 b + d×10 l+1 -1 a=(d+1)×10 l 10 l+1 -1 9 a ) .
Let us first find an appropriate lower bound of P (d,φ d (n)) :

P (d,φ d (n)) ≥ 1 d × 10 n n-1 l=1 ( (d+1)×10 l -1 b=d×10 l 1 - (9d -1)10 l 9 (d+1)×10 l -1 b=d×10 l 1 b + 10 l+1 -1 9 d×10 l+1 -1 a=(d+1)×10 l 1 a ) ≥ 1 d × 10 n n-1 l=1 10 l - 9d -1 9 10 l ln( (d + 1) × 10 l -1 d × 10 l -1 ) + 10 l+1 -1 9 ln( d × 10 l+1 (d + 1) × 10 l ) ,
knowing that for all integers (p, q), such that 1 < p < q:

ln( q + 1 p ) ≤ q k=p 1 k ≤ ln( q p -1
) .

(

) 5
Therefore we have:

P (d,φ d (n)) ≥ 1 d × 10 n n-1 l=1 10 l - 9d -1 9 n-1 l=1 10 l ln( d + 1 d ) + ln( 10 l -1 d+1 10 l -1 d ) + ln( 10d d+1 ) 9 n-1 l=1 (10 l+1 -1) ≥ 1 d × 10 n 10(10 n-1 -1) 9 - (9d -1) ln( d+1 d ) 9 × 10(10 n-1 -1) 9 - 9d -1 9 n-1 l=1 10 l ln(1 + 1 d(d+1) 10 l -1 d ) + ln( 10d d+1 ) 9 × ( 10 2 (10 n-1 -1) 9 -(n -1)) ≥ 9 -(9d -1) ln( d+1 d ) + 10 ln( 10d d+1 ) 81d - 9 + 10 ln( 10d d+1 ) + 9 ln( 10d d+1 ) 10 n 81d × 10 n-1 - 9d -1 9 n-1 l=1 10 l ln(1 + 1 d+1 d×10 l -1 ) d × 10 n
we know that for all x ∈] -1; +∞[ we have: ln(1 + x) ≤ x, thus:

- 9d -1 9 n-1 l=1 10 l ln(1 + 1 d+1 d×10 l -1 ) d × 10 n ≥ - 1 10 n n-1 l=1 1 d+1 × 10 l d × 10 l -1 ≥ - 1 10 n n-1 l=1 1 ≥ - n 10 n .
Consequently, we obtain this lower bound:

P (d,φ d (n)) ≥ 9 -(9d -1) ln( d+1 d ) + 10 ln( 10d d+1 ) 81d - 90 + 100 ln( 10d d+1 ) + 9 ln( 10d d+1 )n + 81dn 81d × 10 n . (6)
Let us now find an appropriate upper bound of P (d,φ d (n)) :

P (d,φ d (n)) ≤ 1 d × 10 n -1 n-1 l=0 ( (d+1)×10 l -1 b=d×10 l 1 - (9d -1)10 l -8 9 (d+1)×10 l -1 b=d×10 l 1 b + 10 l+1 9 d×10 l+1 -1 a=(d+1)×10 l 1 a ) ≤ 1 d × 10 n -1 10 n -1 9 - 1 9 n-1 l=1 (9d -1)10 l -8 ln( (d + 1) × 10 l d × 10 l ) + n-1 l=0 10 l+1 9 ln( d × 10 l+1 -1 (d + 1) × 10 l -1
) , thanks to inequalities (5). Thus we get:

P (d,φ d (n)) ≤ 1 d × 10 n -1 10 n 9 - ln( d+1 d ) 9 n-1 l=1 (9d -1)10 l -8 + n-1 l=0 10 l+1 9 (ln( 10d d + 1 ) + ln( 10 l -1 10d 10 l -1 d+1
))

≤ 1 d × 10 n -1 10 n 9 - ln( d+1 d ) 9 (9d -1)10(10 n-1 -1) 9 -8n + 10 ln( 10d d+1 ) 9 × 10 n -1 9 + 10 9 n-1 l=0 10 l ln(1 + 1 d+1 -1 10d 10 l -1 d+1 ) ≤ 1 d(10 n -1 d ) 10 n -1 d + 1 d 9 - ln( d+1 d ) 9 × ( (9d -1)(10 n -1 d + 1 d -10) 9 -8n) + 10 ln( 10d d+1 ) 9 × 10 n -1 d + 1 d -1 9 + 10 9 n-1 l=0 10 l ln(1 + 1 d+1 -1 10d 10 l -1 d+1 ) ≤ 9 -(9d -1) ln( d+1 d ) + 10 ln( 10d d+1 ) 81d + 9 d + 90d ln( d+1 d ) + 72n ln( d+1 d ) + 10 ln( 10d d+1 ) d 81(d × 10 n -1) + 10 9(d × 10 n -1) n-1 l=0 10 l (1 -d+1 10d ) (d + 1)10 l -1 ≤ 9 -(9d -1) ln( d+1 d ) + 10 ln( 10d d+1 ) 81d + 9 d + 90d ln( d+1 d ) + 72n ln( d+1 d ) + 10 ln( 10d d+1 ) d 81(d × 10 n -1) + 10 9(d × 10 n -1) n-1 l=0 1 .
The last step is easy to demonstrate even for l = 0. Consequently, we obtain this upper bound: The bound just above and the one brought to light in inequality ( 6), added to the following limits: Let us denote by α d the limit of (P (d,φ d (n)) ) n∈N * :

P (d,φ d (n)) ≤ 9 -(9d -
α d = 9 + 10 ln 10 + 9(d + 1) ln(1 -1 d+1 ) 81d
.

Here is the first values of

P (1,φ(n)) (α 1 ≈ 0.241):
Here is a few values of

P (d,φ d (n)) , for d ∈ 2, 9 : n φ(n) P (i,φ(n)) P (i,φ(n)) -

The second subsequence

The second subsequence we will consider is (P (d,ψ d (n)) ) n∈N * where ψ d is the function from N * to N that maps n to (d + 1) × 10 n -1. We get the following result:

Proposition 3. Proof. Let n be a positive integer such that n ≥ 2. According to Proposition 2.1, we have:

P (d,ψ d (n)) = P (d,(d+1)×10 n -1) = 1 (d + 1) × 10 n -1 n l=0 (d+1)×10 l -1 b=d×10 l b - (9d-1)10 l -8 9 b + n-1 l=0 d×10 l+1 -1 a=(d+1)×10 l 10 l+1 -1 9 a .
Let us first find an appropriate lower bound of P (d,ψ d (n)) in a way very similar to that used in the proof of Proposition 3.1:

P (d,ψ d (n)) ≥ 1 (d + 1) × 10 n 10(10 n -1) 9 - (9d -1) ln( d+1 d ) 9 × 10(10 n -1) 9 - 9d -1 9 n l=1 10 l ln(1 + 1 d(d+1) 10 l -1 d ) + ln( 10d d+1 ) 9 × 10 2 (10 n-1 -1) 9 -(n -1)
Let us now find an appropriate upper bound of P (d,ψ d (n)) using the proof of Proposition 3.1:

P (d,ψ d (n)) ≤ 1 (d + 1) × 10 n -1 n l=0 10 l - ln( d+1 d ) 9 n l=0 (9d -1)10 l -8 + n-1 l=0 10 l+1 9 ln( 10d d + 1 ) + ln(1 + 1 d+1 -1 10d 10 l -1 d+1 ) ≤ 1 (d + 1)(10 n -1 d+1 ) 10 n+1 -1 9 - ln( d+1 d ) 9 (9d -1)(10 n+1 -1) 9 -8(n + 1)) + 10 ln( 10d d+1 ) 9 × 10 n -1 9 + 10 9 n-1 l=0 10 l (1 -d+1 10d ) (d + 1)10 l -1 ≤ 1 (d + 1)(10 n -1 d+1 ) 10(10 n -1 d+1 ) -1 + 10 d+1 9 + 10 ln( 10d d+1 )(10 n -1 d+1 + 1 d+1 -1) 81 - ln( d+1 d ) 9 (9d -1)(10(10 n -1 d+1 ) -1 + 10 d+1 ) 9
-8(n + 1) + 10 9

n-1 l=0

1 .

Thereby:

P (d,ψ d (n)) ≤ 90 -10(9d -1) ln( d+1 d ) + 10 ln( 10d d+1 ) 81(d + 1) + 90 d+1 + 72(n + 1) ln( d+1 d ) + 90n 81 (d + 1) × 10 n -1 . (8)
Bounds brought to light in inequalities ( 7) and ( 8) and the fact that: Let us denote by β d the limit of (P (d,ψ d (n)) ) n∈N * :

β d = 10 9 + ln 10 + 9d ln(1 -1 d+1 ) 81(d + 1)
.

Here is the first values of

P (1,ψ(n)) (β 1 ≈ 0.313):
Here is a few values of

P (d,ψ d (n)) , for d ∈ 2, 9 :
4 The graph of (P

(d,n) ) n∈N *
Let us first plot the graph of the sequence (P (1,n) ) n∈N * for values of n from 1 to 1200 (Figure 1). Then we plot a second graph of P (1,n) versus log(n), for n ∈ 1, 32000 (Figure 2). On this graph, the four dots represented by red circles are associated with the first values of (P Through Figure 2, it is clear that the proportion of 1 as leading digit structurally fluctuate and does not follow Benford's Law.

Let us additionally plot graphs of sequences (P (d,n) ) n∈N * for values of n from 1 to 400 (Figure 3). Then we plot graphs of P (d,n) versus log(n), for n ∈ 1, 32000 (Figure 4).

Through Figure 4, it is once more clear that the proportion of each d as leading digit, d ∈ 1, 9 , structurally fluctuate and does not follow Benford's Law.

For each d ∈ 1, 9 the values seem to fluctuate between two values, under the following constraint: Proposition 4.1. For all n ∈ N * such that n ≥ 10 and for all (p, q) ∈ 1, 9 2 such that p < q, we have:

P (p,n) > P (q,n) .
The relative position of graphs of P (d,n) , for d ∈ 1, 9 , can be observed on Figures 3 and4.

Proof. For all n ∈ N * such that n ≥ 10 and for all (p, q) ∈ 1, 9 2 such that p < q, we denote by k p and k q the positive integers such that k p = min{i ∈ N : p×10 i > n} and k q = min{i ∈ N : q × 10 i > n}. We note that k p ≥ k q . For (r, l) ∈ N 2 , let A l and B r,l be real numbers such that A l = 10 l+1 -1 9 and B r,l = (9r-1)×10 l -8

9

. We consider that A -1 = 0. Third cases can be distinguished: In the first case (q + 1) × 10 kp-2 ≤ p × 10 kp-1 ≤ n < (p + 1) × 10 kp-1 ≤ q × 10 kp-1 . Thus we have k p = k q + 1. Thanks to Proposition 2.1 we obtain:

P (p,n) -P (q,n) = 1 n kp -2 l=0 ( (p+1)×10 l -1 b=p×10 l b -B p,l b + p×10 l+1 -1 a=(p+1)×10 l A l a ) + n b=p×10 kp -1 b -B p,kp -1 b - 1 n kp -3 l=0 ( (q+1)×10 l -1 b=q×10 l b -B q,l b + q×10 l+1 -1 a=(q+1)×10 l A l a ) + (q+1)×10 kp-2 -1 b=q×10 kp -2 b -B q,kp-2 b + n a=(q+1)×10 kp -2 A kp -2 a = 1 n kp -3 l=0 (p+1)×10 l -1 b=p×10 l b -B p,l b + p×10 l+1 -1 a=(p+1)×10 l A l a - (q+1)×10 l -1 b=q×10 l b -B q,l b - q×10 l+1 -1 a=(q+1)×10 l A l a + (p+1)×10 kp -2 -1 b=p×10 kp -2 b -B p,kp -2 b + p×10 kp-1 -1 a=(p+1)×10 kp -2 A kp -2 a + n b=p×10 kp-1 b -B p,kp b - (q+1)×10 kp -2 -1 b=q×10 kp -2 b -B q,kp b - n a=(q+1)×10 kp -2 A kp -2 a = 1 n kp -2 l=0 (p+1)×10 l -1 b=p×10 l b -B p,l -A l-1 b + q×10 l -1 a=(p+1)×10 l A l -A l-1 a + (q+1)×10 l -1 b=q×10 l A l -(b -B q,l ) b + p×10 l+1 -1 a=(q+1)×10 l A l -A l a + n b=p×10 kp -1 b -B p,kp -1 -A kp -2 b .
Furthermore, for all (r, l) ∈ N 2 , we have A l > A l-1 and:

∀s ∈ r × 10 l , (r + 1) × 10 l+1 -1 , s -B r,l -A l-1 = s - (9r -1) × 10 l -8 9 - 10 l -1 9 = s -(r × 10 l -1) > 0 A l -(s -B r,l ) = 10 l+1 -1 9 -s + (9r -1) × 10 l -8 9 = (r + 1) × 10 l -1 -s ≥ 0 .
Consequently, in this case, P (p,n) -P (q,n) > 0.

In the second case (p + 1) × 10 kp-1 ≤ n < q × 10 kp-1 . Thus we have k p = k q + 1. Thanks to Proposition 2.1 we obtain:

P (p,n) -P (q,n) = 1 n kp -2 l=0 ( (p+1)×10 l -1 b=p×10 l b -B p,l b + p×10 l+1 -1 a=(p+1)×10 l A l a ) + (p+1)×10 kp -1 -1 b=p×10 kp -1 b -B p,kp-1 b + n a=(p+1)×10 kp -1 A kp -1 a - 1 n kp-3 l=0 ( (q+1)×10 l -1 b=q×10 l b -B q,l b + q×10 l+1 -1 a=(q+1)×10 l A l a ) + (q+1)×10 kp-2 -1 b=q×10 kp -2 b -B q,kp-2 b + n a=(q+1)×10 kp -2 A kp -2 a = 1 n kp -2 l=0 (p+1)×10 l -1 b=p×10 l b -B p,l -A l-1 b + q×10 l -1 a=(p+1)×10 l A l -A l-1 a + (q+1)×10 l -1 b=q×10 l A l -(b -B q,l ) b + (p+1)×10 kp-1 -1 b=p×10 kp -1 b -B p,kp -1 -A kp -2 b + n a=(p+1)×10 kp -1 A kp -1 -A kp -2 a .
Consequently, in this case, P (p,n) -P (q,n) > 0.

In the third case q × 10 kp-1 ≤ n < p × 10 kp . Thus we have k p = k q . Thanks to Proposition 2.1 we obtain:

P (p,n) -P (q,n) = 1 n kp -2 l=0 (p+1)×10 l -1 b=p×10 l b -B p,l -A l-1 b + q×10 l -1 a=(p+1)×10 l A l -A l-1 a + (q+1)×10 l -1 b=q×10 l A l -(b -B q,l ) b + (p+1)×10 kp-1 -1 b=p×10 kp -1 b -B p,kp -1 -A kp -2 b + q×10 kp -1 a=(p+1)×10 kp -1 A kp -1 -A kp -2 a + min((q+1)×10 kp -1 -1,n) b=q×10 kp -1 A l -(b -B q,l ) b .
Consequently, in this latter case, P (p,n) -P (q,n) > 0.

Remark 4.2. For n ∈ N * , we have, if d > n, P (d,n) = 0. Hence for all n ∈ N * and for all (p, q) ∈ 1, 9 2 such that p < q, we have:

P (p,n) ≥ P (q,n) .
Let us from now on denote by k d the positive integer such that k d = min{i ∈ N : d × 10 i > n}. Through Figures 2 and4, "delayed effects" appear to exist, in particular after a long sequence of numbers whose leading digit is d. Let us examine these features in more detail. We study for this purpose the increasing or decreasing nature of sequences (P (d,n) ) n∈N * , for d ∈ 1, 9 :

Proposition 4.3. i ∈ {0, 1}. P (d,n+1) -P (d,n) = 1 n + 1 ( c i n + 1 -P (d,n) ),
where:

ci = n + 1 -(9d-1)10 k d -8 9 if n + 1 ∈ d × 10 k d , (d + 1) × 10 k d -1 10 k d +1 -1 9 if n + 1 ∈ (d + 1) × 10 k d , d × 10 k d +1 -1
Proof. It is based on the formulas of Proposition 2.1. Indeed:

P (d,n+1) = 1 n + 1 (nP (d,n) + c i n + 1 ) = P (d,n) + 1 n + 1 ( c i n + 1 -P (d,n) ).
Through this proposition the obvious condition regarding the increasing or decreasing nature of the sequence is underscored: whether P (d,n) value, for d ∈ 1, 9 , is less or greater than the appropriate ci n+1 value, for n ∈ N * and i ∈ {0, 1}. Finding the approximate values of n for which the increasing or decreasing nature of the sequence (P (d,n) ) n∈N * appears is henceforth the aim of that section. We first provide a proposition similar to the previous one:

P (d,n) = 1 n P (d,d×10 k d -1) × (d × 10 k d -1) + n b=d×10 k d b -(9d-1)10 k d -8 9
b .

If n ∈ (d + 1) × 10 k d , d × 10 k d +1 -1 , then we have:

P (d,n) = 1 n P (d,(d+1)×10 k d -1) × ((d + 1) × 10 k d -1) + n a=(d+1)×10 k d 10 k d +1 -1 9
a .

Proof. Results are directly derived from Proposition 2.1.

We now consider the sequence ( P (d,n) ) n∈N * defined as follows. If n ∈ d × 10 k d , (d + 1) × 10 k d -1 :

P (d,n) = d(α d -1) 10 k d n + 1 + 9d -1 9 10 k d n ln( 10 k d n ) + ln d .
If n ∈ (d + 1) × 10 k d , d × 10 k d +1 -1 :

P (d,n) = (d + 1)β d 10 k d n - 10 9
10 k d n ln( 10 k d n ) + ln(d + 1) .
We denote by γ d a real number such that

γ d ∈]1; 1+ 1 d [ and I (d,γ d ) the set such that I (d,γ d ) = +∞ i=1 dγ d ×10 i , (d+1)×10 i -1 ∪ (d+1)γ d ×10 i , d×10 i+1 -1 .
The below proposition can thereupon be stated: Proposition 4.5.

P (d,n) ∼ n→+∞ n∈I (d,γ d ) P (d,n) .
Proof. Let us study both cases. In the first one, n ∈ γ d d × 10 k d , (d + 1) × 10 k d -1 . Let I 1 be the interval such that, I 1 = +∞ i=1 dγ d × 10 i , (d + 1) × 10 i -1 . We have:

1 d + 1 ≤ 10 k d n ≤ 1 γ d d . (9)
Before we go any further, let us prove the following lemma:

Lemma 4.6. For n > d × 10 k d :

ln( n 10 k d ) -ln d + ln(1 + 1 n ) ≤ n b=d×10 k d 1 b ≤ ln( n 10 k d ) -ln d + ln(1 + 1 d × 10 k d -1
) .

Proof. This result is directly related to inequalities 5. Indeed for n > d × 10 k d > 1, we have:

ln( n + 1 d × 10 k d ) ≤ n b=d×10 k d 1 b ≤ ln( n d × 10 k d -1
) .

Thanks to Proposition 4.4 we get:

P (d,n) = P (d,d×10 k d -1) d × 10 k d -1 n + n -d × 10 k d + 1 n - (9d -1)10 k d -8 9n n b=d×10 k d 1 b = d(P (d,d×10 k d -1) -1) 10 k d n + 1 - (9d -1)10 k d 9n n b=d×10 k d 1 b - P (1,d×10 k d -1) n + 1 n + 8 9n n b=d×10 k d 1 b .
Thanks to inequalities ( 9) and knowing that We have:

n b=d×10 k d 1 b = ln( n 10 k d ) -ln d + o n→+∞ n∈I1 (1)
1 10d ≤ 10 k d n ≤ 1 (d + 1)γ d . ( 10
)
Before we go any further, let us prove this additional lemma:

Lemma 4.7. For n > (d + 1) × 10 k d : ln( n 10 k d ) -ln(d + 1) + ln(1 + 1 n ) ≤ n a=(d+1)×10 k d 1 a and, n a=(d+1)×10 k d 1 a ≤ ln( n 10 k d ) -ln(d + 1) + ln(1 + 1 (d + 1) × 10 k d -1
) .

Proof. This result is directly related to inequalities 5. Indeed for n > (d + 1) × 10 k d > 1, we have:

ln( n + 1 (d + 1) × 10 k d ) ≤ n a=(d+1)×10 k d 1 a ≤ ln( n (d + 1) × 10 k d -1
) .

Thanks to Proposition 4.4 we get:

P (d,n) = P (d,(d+1)×10 k d -1) (d + 1) × 10 k d -1 n + 1 9 × 10 k d +1 -1 n n a=(d+1)×10 k d 1 a = (d + 1)P (d,(d+1)×10 k d -1) 10 k d n + 10 9 × 10 k d n n a=(d+1)×10 k d 1 a - P (1,(d+1)×10 k d -1) n - 1 9n n a=(d+1)×10 k d 1 a .
Thanks the inequalities [START_REF] Hill | A statistical derivation of the significant-digit law[END_REF] and knowing that

n a=(d+1)×10 k d 1 a = ln( n 10 k d ) - ln(d + 1) + o n→+∞ n∈I2 (
P (d,n) ∼ n→+∞ n∈I2 (d + 1)β d 10 k d n - 10 9
10 k d n ln( 10 k d n ) + ln(d + 1) .
To find the approximate values of n for which the increasing or decreasing nature of (P (d,n) ) n∈N * appears, we need to study two distinct functions:

Lemma 4.8. The minimum m d of the function f d from [ 1 d+1 , 1 d ] that maps x onto d(α d -1)x + 1 + 9d-1 9 x(ln(x) + ln d) is reached when: x = 10 -10 9(9d-1) 1 -1 d+1 -d+1 9d-1 d .

Its value is m

d = 1 - (9d-1)10 - 10
9(9d-1) 1-1 d+1 -d+1 9d-1 9d . Proof. ∀x ∈ [ 1 d+1 , 1 d ], f d (x) = d(α d -1) + 9d-1 9 (1 + ln(x) + ln d).
By solving the equation d(α d -1) + 9d-1 9 (1 + ln(x) + ln d) = 0, we have:

x = 10 -10 9(9d-1) 1 -1 d+1 -d+1 9d-1 d . Finally f d ( 10 - 10
9(9d-1) 1-1 d+1 -d+1 9d-1 d ) = 1 - (9d-1)10 - 10
9(9d-1) 1-1 d+1 -d+1 9d-1 9d
.

Remark 4.9. We note that: 1-

f d ( 1 d ) = α d -1 + 1 = α d and f d ( 1 d + 1 ) = (α d -1) × d d + 1 + 1 + 9d -1 9 × 1 d +
1 d+1 d 9 (d+1)
.

Remark 4.11. We note that:

g d ( 1 d + 1 ) = β d and g d ( 1 10d ) = (d + 1)β d × 1 10d - 10 9 × 1 10d ln( 1 10d
) + ln(d + 1)

= 9 + ln 10 + 9d ln(1 -1 d+1 ) 81d - 1 9d ln( 1 10d
) + ln(d + 1) = 9 + 10 ln 10 + 9(d + 1) ln(1

-1 d+1 ) 81d = α d .
Using the approximation of Proposition 4.5, the definition of (

P (d,n) ) n∈N * (for n ∈ (d + 1) × 10 ki , d × 10 k d +1 -1 )
M (1,i) M (1,i) P (1,M (1,i
= 10 1 9 1-1 2+1 2 2+1 = 4×10 1 9
27 , i.e. M (2,4) ≈ 52263.

Both values of maxima and ranks for which these maxima are reached are correctly approximated by our results as illustrated in Table 8.

When considering the random experiment defined in the beginning of the article, we have thus determined the values of the proportions of selected numbers whose leading digit is d and its bounds: these values seem to fluctuate between m d and M d .

Central values

From previous Figures, we notice that there exist fluctuations in the graph of (P (d,n) ) n∈N * . We can calculate over each "pseudo-cycle", i.e. for all n ∈ d × 10 i , d × 10 i+1 -1 where i ∈ N, the mean value C (d,i) of P (d,n) . For example, we obtain: We will now consider the sequence (C (d,n) ) n∈N and will demonstrate that it converges. Before we go any further, let us prove the following lemma: Lemma 5.2. For all (p, q) ∈ N 2 , such that 4 < p < q, we have: ln q(p -1) ln( q p-1 ) 2 ≤ q n=p n ln n ≤ ln (q + 1)p ln( q+1 p ) 2 .

Proof. The function from [3; +∞[ to R that maps x onto x ln x is increasing on [3; +∞[. Thus:

q p-1 x ln x dx ≤ q n=p n ln n ≤ q+1 p x ln x dx (ln x) 2 2 q p-1 ≤ q n=p n ln n ≤ (ln x) 2 2 q+1 p .
The result follows.

The below proposition can thereupon be stated:

Proposition 5.3. Proof. Let E be a real number such that 0 < E < 1. For all n ∈ N * , we have:

C (d,n) ∼ +∞ (18d(α d -1) -(9d -
C (d,n) = 1 9d × 10 n d×10 n+1 -1 i=d×10 n P (d,i) = 1 9d × 10 n (d+E)×10 n -1 i=d×10 n P (d,i) + (d+1)×10 n -1 i= (d+E)×10 n P (d,i) + ((d+1)+E)×10 n -1 i=(d+1)×10 n P (d,i) + d×10 n+1 -1 i= ((d+1)+E)×10 n P (d,i) .
For all n ∈ N * , let us consider C (d,n) = 1 9d×10 n d×10 n+1 -1 i=d×10 n P (d,i) . We know that:

∀i ∈ N * , |P (d,i) -P (d,i) | ≤ |P (d,i) | + | P (d,i) | ≤ 2 .
There also exists an integer t such that for all i ≥ t and i ∈ I E , |P (d,i) -P (d,i) | ≤ E (see Proposition 4.5).

Thus, for all n ∈ N * such that 10 n ≥ t, we have:

|C (d,n) -C (d,n) | ≤ 1 9d × 10 n d×10 n+1 -1 i=d×10 n |P (d,i) -P (d,i) | ≤ 1 9d × 10 n (d+E)×10 n -1 i=d×10 n |P (d,i) -P (d,i) | + (d+1)×10 n -1 i= (d+E)×10 n |P (d,i) -P (d,i) | + ((d+1)+E)×10 n -1 i=(d+1)×10 n |P (d,i) -P (d,i) | + d×10 n+1 -1 i= ((d+1)+E)×10 n |P (d,i) -P (d,i) | ≤ 1 9d × 10 n 2( (d + E) × 10 n -d × 10 n ) + ((d + 1) × 10 n -(d + E) × 10 n )E + 2( ((d + 1) + E) × 10 n -(d + 1) × 10 n ) + (d × 10 n+1 -((d + 1) + E) × 10 n )E ≤ 1 9d × 10 n 2E × 10 n + 10 n E + 2E × 10 n + 8 × 10 n E ≤ (5 + (9d -1))E 9 . Consequently C (d,n) ∼ +∞ C (d,n
) . We will henceforth study ( C (d,n) ) n∈N . We have for all n ∈ N * :

C (d,n) = 1 9d (d+1)×10 n -1 i=d×10 n P (d,i) 10 n + d×10 n+1 -1 i=(d+1)×10 n P (d,i) 10 n . (11)
We consider the first term of this sum:

1 9d (d+1)×10 n -1 i=d×10 n P (d,i) 10 n = 1 9d (d+1)×10 n -1 i=d×10 n d(α d -1) i + 1 10 n + 9d -1 9d (ln( 10 n i ) + ln d) = 9d(α d -1) + (9d -1)n ln(10) + (9d -1) ln d 81d (d+1)×10 n -1 i=d×10 n 1 i - 9d -1 81d (d+1)×10 n -1 i=d×10 n ln i i + 1 9d .
The proofs of Lemma 4.6 and Proposition 3.1 allow us to state:

ln( (d + 1) × 10 n d × 10 n ) ≤ (d+1)×10 n -1 i=d×10 n 1 i ≤ ln( (d + 1) × 10 n -1 d × 10 n -1 ) ln( d + 1 d ) ≤ (d+1)×10 n -1 i=d×10 n 1 i ≤ ln( d + 1 d ) + ln(1 + 1 d(d+1) 10 n -1 d ) , i.e. (d+1)×10 n -1 i=d×10 n 1 i = ln( d+1 d ) + O +∞ ( 1 10
10 n -1 d )) 2 ≤ (d+1)×10 n -1 i=d×10 n ln i i .
We have: We obtain:

ln(d(d + 1)) + 2n ln(10) + O +∞ ( 1 10 n ) ln( d+1 d ) + O +∞ ( 1 10 n ) 2 ≤ (d+1)×10 n -1 i=d×10 n ln i i (ln(d(d + 1)) + 2n ln(10)) ln( d+1 d ) + o +∞ (1) 2 ≤ (d+1)×10 n -1 i=d×10 n ln i i .
Thanks to Lemma 5.2 we also have:

(d+1)×10 n -1 i=d×10 n ln i i ≤ ln d(d + 1) × 10 2n ln( (d+1)×10 n d×10 n ) 2 (d+1)×10 n -1 i=d×10 n ln i i ≤ (ln(d(d + 1)) + 2n ln(10)) ln( d+1 d ) 2 .
Finally the first term of equality [START_REF] Jolissaint | Loi de benford, relations de récurrence et suites équidistribuées[END_REF] .

We consider henceforth the second term of the equality (11):

1 9d d×10 n+1 -1 i=(d+1)×10 n P (d,i) 10 n = 1 9d d×10 n+1 -1 i=(d+1)×10 n (d + 1)β d 1 i - 10
(d + 1) × 10 n ) ≤ d×10 n+1 -1 i=(d+1)×10 n 1 i ≤ ln( d × 10 n+1 -1 (d + 1) × 10 n -1 ) ln( d + 1 ) ≤ d×10 n+1 -1 i=(d+1)×10 n 1 i ≤ ln( 10d d + 1 ) + ln(1 + 1 d+1 -1 10i 10 n -1 d+1
) ,

i.e.

d×10 n+1 -1 i=(d+1)×10 n 1 i = ln( 10d d+1 ) + O +∞ ( 1 10
n ). Then thanks to Lemma 5.2:

ln (d × 10 n+1 -1)((d + 1) × 10 n -1) ln( d×10 n+1 -1 (d+1)×10 n -1 ) 2 ≤ d×10 n+1 -1 i=(d+1)×10 n ln i i ln(10d(d + 1) × 10 2n -(11d + 1) × 10 n + 1)(ln( 10d d+1 ) + ln(1 + 1 d+1 -1 10d 10 n -1 d+1
))

2

≤ d×10 n+1 -1 i=(d+1)×10 n ln i i .
We have:

22

Thanks to Lemma 5.2 we also have:

d×10 n+1 -1 i=(d+1)×10 n ln i i ≤ ln(d(d + 1) × 10 2n+1 ) ln( 10d d+1 ) 2 d×10 n+1 -1 i=(d+1)×10 n ln i i ≤ (ln(10d(d + 1)) + 2n ln(10)) ln( 10d d+1 ) 2 .
Finally the second term of equality (11) verifies: Hence:

C (1,n) ∼ +∞ (18d(α d -1) -(9d -1) ln( d+1 d )) ln( d+1 d ) + 18 162d + (9(d + 1)β d + 5 ln( 10d d+1 )) ln( 10d d+1 ) 81d ∼ +∞ (18d(α d -1) -(9d -1) ln( d+1 d )) ln( d+1 d ) + 18 + 2(9(d + 1)β d + 5 ln( 10d d+1 )) ln( 10d d+1 ) 162d .
The result follows.

Let us denote by C d the limit of (C (d,n) ) n∈N . Note that our first choice of "pseudo-cycle" give more weight at the values of proportion situated at the end of intervals d×10 i , d×10 i+1 -1 where i ∈ N. We can also have defined the sequence ( C(d,n) ) n∈N as follows:

C(d,n) = 1 9(d + 1) × 10 n (d+1)×10 n+1 -1 i=(d+1)×10 n P (1,i) .
If so, we would state that the limit Cd of the sequence ( C(d,n) ) n∈N is: (

) 12
The first term of this sum is similar to the second term of the associated sum in the proof of Proposition 5.3. The second one verifies: Finally the second term of equality ( 12) is: The result follows.

Once more, means values over both "pseudo-cycles" are very close to the theoric value highlighted by Benford: log(1 + 1 d ) ( [START_REF] Benford | The law of anomalous numbers[END_REF]). Indeed, according to Hill ([10]), it is absolutely normal. In a way, it can be considered as an equivalent to the central limit theorem ([8]).

Conclusion

To conclude, through our model, we have seen that the proportion of d as leading digit, d ∈ 1, 9 , in certain naturally occurring collections of data is more likely to follow a law whose probability distribution is (d, P (d,n) ) d∈ 1,9 , where n is the smaller integer upper bound of the physical, biological or economical quantities considered, rather than Benford's Law. These probability distributions fluctuate around Benford's value as can be seen in the literature (see [START_REF] Knuth | The Art of Computer Programming 2[END_REF], [START_REF] Burke | Benford's law and physical constants: the distribution of initial digits[END_REF], [START_REF] Nigrini | Assessing the integrity of tabulated demographic data[END_REF] or [START_REF] Friar | Genome sizes and the Benford distribution[END_REF] for example) in accordance with our model. Knowing beforehand the value of the upper bound n can be a way to find a better adjusted law than Benford's one.

The results of the article would have been the same in terms of fluctuations of the proportion of d ∈ 1, 9 as leading digit, of limits of subsequences, or of results on central values, if our discrete uniform distributions uniformly randomly selected were lower bounded by a positive integer different from 0: first terms in proportion formulas become rapidly negligible. Through our model we understand that the predominance of 1 as first digit (followed by those of 2 and so on) is all but surprising in experimental data: it is only due to the fact that, in the lexicographical order, 1 appears before 2, 2 appears before 3, etc.

However the limits of our model rest on the assumption that the random variables used to obtain our data are not the same and follow discrete uniform distributions that are uniformly randomly selected. In certain naturally occurring collections of data it cannot conceivably be justified. Studying the cases where the random variables follow other distributions (and not necessarily randomly selected) sketch some avenues for future research on the subject.

Examples 2 . 2 .

 22 If n = 20, we have k = 2. The value of P(L 20 = 1) is then (second case of Proposition 2.

Proposition 3 . 1 .

 31 The subsequence (P (d,φ d (n)) ) n∈N * converges to: 9 -(9d -1) ln( d+1 d ) + 10 ln(

  × 10 n -1) .

2 .

 2 The subsequence (P (d,ψ d (n)) ) n∈N * converges to: 10 9 -(9d -1) ln( d+1 d ) + ln( 10d d+1 ) 81(d + 1) = 10 9 + ln 10 + 9d ln(1 -1 d+1 ) 81(d + 1).

≥ 90 - 1 (

 901 10(9d -1) ln( d+1 d ) + 10 ln( 10d d+1 ) 81(d + 1) -90 + 100 ln( 10d d+1 ) + 9n ln( 10d d+1 ) 81(d + 1) × 10 n d + 1) × 10 n . Then: P (d,ψ d (n)) ≥ 90 -10(9d -1) ln( d+1 d ) + 10 ln( 10d d+1 ) 81(d + 1) -90 + 100 ln( 10d d+1 ) + 9n ln( 10d d+1 ) + 81dn (d + 1) × 10 n .

  ln( 10d d+1 ) + 9n ln( 10d d+1 ) + 81dn (d + 1) × 10 n ) = 0 and lim n→+∞ ( 90 d+1 + 72(n + 1) ln( d+1 d ) + 90n 81 (d + 1) × 10 n -1 ) = 0 lead to the expected result.

  associated with the subsequence (P (1,ψ(n)) ) n∈N * . Their distances to the samecoloured horizontal dotted asymptote tend towards 0. The linear equations of these lines are y = α 1 and y = β 1 respectively (see Propositions 3.1 and 3.2).

Figure 1 :

 1 Figure 1: Graph of (P (1,n) ) n∈N * .

Figure 2 :

 2 Figure 2: Graph of P (1,n) versus log(n), values of P (1,φ(n)) being in red and those of P (1,ψ(n)) being in deep blue. Limits α 1 and β 1 of these two subsequences are represented by horizontal asymptotes.

Figure 3 :

 3 Figure 3: Graphs of (P (d,n) ) n∈N * , for d ∈ 1, 9 .

Figure 4 :

 4 Figure 4: For d ∈ 1, 9 , graphs of P (d,n) versus log(n). Note that points have not been all represented.

Proposition 4 . 4 .

 44 If n ∈ d × 10 k d , (d + 1) × 10 k d -1 , then we have:

  (thanks to Lemma 4.6, n d×10 k d being greater than or equal to γ d > 1) and that P (d,d×10 k d -1) ∼ +∞ α d (see Proposition 3.1), we finally have: In the second case n ∈ (d + 1)γ d × 10 k d , d × 10 k d +1 -1 . Let I 2 be the interval such that, I 2 = +∞ i=1 (d + 1)γ d × 10 i , d × 10 i+1 -1 .

3 Lemma 4 . 10 .

 3410 The maximum M d of the function g d from [ 1 10d , 1 d+1 ] that maps x onto (d + 1)β d x -10 9 x ln(x) + ln(d + 1) is reached when x = g d (x) = (d + 1)β d -10 9 (1 + ln(x) + ln(d + 1)). By solving the equation (d+1)β d -10 9 (1+ln(x)+ln(d+1)) = 0, we obtain x =

  and the properties of g d (Lemma 4.10), we can approximate the ranks of local maxima. Let i be a strictly positive integer. Let us denote by M (d,i) the rank of the local maximum of P (d,n) in (d + 1) × 10 i , (d + 1) × 10 i+1 -1 and by M (d,i) the estimate of this rank. Let us gather in the following table the first values of those ranks and values of local maxima, for d = 1.

  i

P

  (9,i) ≈ 0.043 .

1 )

 1 n ). Then thanks to Lemma 5.2: ln ((d + 1) × 10 n -1)(d × 10 n -1) ln( × 10 2n -(2d + 1) × 10 n + 1)(ln( d+1 d ) + ln(1 + 1 d(d+1)

  ln(d(d + 1) × 10 2n -(2d + 1) × 10 n + 1) = ln(d(d + 1) × 10 2n ) + ln(1 -(2d + 1) -1 10 n d(d + 1) × 10 n ) = ln(d(d + 1)) + 2n ln(10) + O +∞ ( 2d + 1 d(d + 1) × 10 n ) .

1 )

 1 × 10 n+1 ) .

(

  90d(α d -1) -5(9d -1) ln( d+1 d )) ln( d+1 d ) + 90 + (9(d + 1)β d + 5 ln( 10d d+1 )) ln( 10d d+1 ) 81(d + 1) .Proof.C(d,n) ∼ +∞ C (d,n)where, for all n ∈ N * :

  d -1) + 10(9d -1)((n + 1) ln(10) + ln d)

  α d -1) + 10(9d -1)((n + 1) ln 10 + ln d)) ln( d+1 d ) 81(d + 1) -5(9d -1)(ln(d(d + 1)) + 2(n + 1) ln(10)) ln( α d -1) -5(9d -1) ln( d+1 d )) ln( d+1 d α d -1) -5(9d -1) ln( d+1 d )) ln( d+1 d ) + 90 81(d + 1) ∼ +∞ (90d(α d -1) -5(9d -1) ln( d+1 d )) ln( d+1 d ) + 90 + (9(d + 1)β d + 5 ln( 10d d+1 )) ln( 10d d+1 ) 81(d + 1).

Table 1 :

 1 First five values P (1,φ(n)) and P (1,φ(n)) -α 1 . We round off these values to three significant digits. d P (d,φ d (1)) P (d,φ d (2)) P (d,φ d (3)) P (d,φ d (4))

	α 1

Table 2 :

 2

Values of P (d,φ d (n)) and α d , for n ∈ 1, 4 . These values are rounded to the nearest thousandth.

Table 3 :

 3 First five values P (1,ψ(n)) and P (1,ψ(n)) -β 1 . We round off these values to three significant digits.

	d P (d,ψ d (1)) P (d,ψ d (2)) P (d,ψ d (3)) P (d,ψ d (4))	β d
	2	0.176	0.166	0.165	0.165	0.165
	3	0.110	0.109	0.109	0.109	0.109
	4	0.078	0.080	0.081	0.0.081	0.081
	5	0.060	0.063	0.064	0.064	0.064
	6	0.049	0.052	0.052	0.053	0.053
	7	0.041	0.044	0.045	0.045	0.045
	8	0.035	0.038	0.039	0.039	0.039
	9	0.031	0.034	0.034	0.034	0.034

Table 4 :

 4

Values of P (d,ψ d (n)) and β d , for n ∈ 1, 4 . These values are rounded to the nearest thousandth.

Table 5 :

 5 Using the approximation of Proposition 4.5, the definition of ( P (d,n) ) n∈N * (for n ∈ d × 10 k d , (d + 1) × 10 k d -1 ) and the properties of f d (Lemma 4.8), we can approximate the ranks of local minima. Let i be a strictly positive integer. Let us denote by m (d,i) the rank of the local minimum of P (d,n) in d × 10 i , (d + 1) × 10 i -1 and by m (d,i) the estimate of this rank. Let us gather in the following table the first values of those ranks and values of local minima, for d = 1. First four values of above defined ranks and associated values of local minima of P (1,n) . We round off values of P (1,n) to three significant digits and values of m (1,i) to unity. Note that m 1 ≈ 0.232. Indeed, for n ∈ 10, 19 , our approximation of the rank m (1,1) for which the minimum is reached verifies: 10 1 m (1,1) = 2 Both values of minima (approximately 0.232 according to Lemma 4.8) and ranks for which these minima are reached are correctly approximated by our results as illustrated in Table 5, for d = 1. Let us similarly gather in the below table the first values of those ranks and values of local minima, for d ∈ 2, 9 .

	1	ln(	d d + 1	)

5 36 , i.e. m (1,1) ≈ 12.

Table 6 :

 6 Values of above defined ranks and associated values of local minima of P (d,n) , for n ∈ d × 10 4 , (d + 1) × 10 4 -1 . We round off values of P (d,n) and m d to three significant digits and values of m (d,4) to unity.

	Indeed, for n ∈ 20000, 29999 , our approximation of the rank m (2,4) for
	which the minimum is reached verifies: 10 4 m (2,4) =

Table 7 :

 7 First four values of above defined ranks and associated values of local maxima of P (1,n) . We round off values of P (1,n) to three significant digits and values of M (1,i) to unity. Note that M 1 ≈ 0.359. Indeed, for n ∈ 20, 99 , our approximation of the rank M (1,1) for which the maximum is reached verifies: 10 1 n M = 5

	) )

Table 8 :

 8 Values of above defined ranks and associated values of local maxima of P (d,n) , for n ∈ (d + 1) × 10 4 , d × 10 5 -1 . We round off values of P (d,n) and M d to three significant digits and values of M (d,4) to unity.Indeed, for n ∈ 30000, 199999 , our approximation of the rank M (2,4) for which the maximum is reached verifies:10 4

	M (2,4)

Table 9

 9

	below gathers

Table 9 :

 9 Values of C d , Cd and Benford's Law probabilities, these values being rounded to the nearest thousandth.

2

, i.e. m (2,4) ≈ 21643. Both values of minima and ranks for which these minima are reached are correctly approximated by our results as illustrated in Table 6.

The second function we need to study is defined below:

Appendix: Python script Using Propositions 2.1, we can determine the terms of (P (d,n) ) n∈N * , for d ∈ 1, 9 . To this end, we have created a script with the Python programming language (Python Software Foundation, Python Language Reference, version 3.4. available at http://www.python.org, see [START_REF] Van Rossum | Python tutorial[END_REF]). The implemented function expvalProp has two parameters: the rank n of the wanted term of the sequence and the value ld of the considered leading digit. Here is the used algorithm: