Criteria for the numerical constant recognition

02/28/2020
by   Andrzej Odrzywolek, et al.
Jagiellonian University
0

The need for recognition/approximation of functions in terms of elementary functions/operations emerges in many areas of experimental mathematics, numerical analysis, computer algebra systems, model building, machine learning, approximation and data compression. One of the most underestimated methods is the symbolic regression. In the article, reductionist approach is applied, reducing full problem to constant functions, i.e, pure numbers (decimal, floating-point). However, existing solutions are plagued by lack of solid criteria distinguishing between random formula, matching approximately or literally decimal expansion and probable ”exact” (the best) expression match in the sense of Occam's razor. In particular, convincing STOP criteria for search were never developed. In the article, such a criteria, working in statistical sense, are provided. Recognition process can be viewed as (1) enumeration of all formulas in order of increasing Kolmogorov complexity K (2) random process with appropriate statistical distribution (3) compression of a decimal string. All three approaches are remarkably consistent, and provide essentially the same limit for practical depth of search. Tested unique formulas count must not exceed 1/sigma, where sigma is relative numerical error of the target constant. Beyond that, further search is pointless, because, in the view of approach (1), number of equivalent expressions within error bounds grows exponentially; in view of (2), probability of random match approaches 1; in view of (3) compression ratio much smaller than 1.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 4

page 15

page 16

01/17/2020

FRaZ: A Generic High-Fidelity Fixed-Ratio Lossy Compression Framework for Scientific Floating-point Data

With ever-increasing volumes of scientific floating-point data being pro...
01/27/2021

Probabilistic Error Analysis For Sequential Summation of Real Floating Point Numbers

We derive two probabilistic bounds for the relative forward error in the...
09/23/2016

A Computer Algebra Package for Polynomial Sequence Recognition

The software package developed in the MS thesis research implements func...
11/16/2019

Counting solutions to random CNF formulas

We give the first efficient algorithm to approximately count the number ...
03/30/2021

How to hunt wild constants

There are now several comprehensive web applications, stand-alone comput...
06/18/2020

AI Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity

We present an improved method for symbolic regression that seeks to fit ...
05/11/2020

Computationally Inequivalent Summations and Their Parenthetic Forms

Floating-point addition on a finite-precision machine is not associative...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Experimental mathematics become very important nowadays [1]. Using modern computers we are able to perform sequential brute-force searches of various mathematical objects and models at stunning scale. Multi-core processors [2]

, GPUs, vectorization, FPGA and future quantum computers allow us to reach depths incomprehensible by humans. In essence, this is some form of artificial intelligence

[3], because obtained results are indistinguishable from those provided by reasoning or more sophisticated algorithms [4]. Particularly interesting problem is a numerical constant recognition [5]. We are given decimal expansion, and ask if there is equivalent formula? There is strong demand from users of mathematical software to provide such a feature (identify/Maple [6], nsimplify/SymPy [8], RIES [9], FindFormula [7] ). In typical situation we encounter some decimal number, resulting from numerical simulation or experiment, e.g: 1.82263, and ask ourselves if it is equivalent to some symbolic expression, like e.g: . Unfortunately, such a problem without additional constraints is ill-posed, and provided ”answers” nonsense. Often ridiculously complicated [10]. This is not a surprise, as Cardinality of real/complex transcendental numbers is uncountably infinite (continuum ), while number of formulas and symbols is countable (). Therefore probability for randomly chosen real number to be equivalent to some formula is zero. We must therefore restrict to numbers with finite decimal expansion. This is also practical approach, due to widespread of floating-point hardware, with double precision [11] being de facto standard. However, in this case we encounter infinity from the other side. Assuming floating-point constant is a truncated real number111If it is not, then it is a rational number, e.g 1.82263 = 182263/100000., we can produce arbitrary number of formulas which numerically differ only at more distant decimal places: etc. Fortunately, one still can ask which one of these formulas has lowest Kolmogorov complexity [12]. But to do this, we must specify ”programing language” first. This lead to the following formulation of the numerical constant recognition problem. Imagine a person with hand-held scientific calculator, who secretly push a sequence of buttons, and print out the decimal/numerical result. We ask: could this process be reversed? In other words, given numerical result, are we able to recover sequence of some particular calculator buttons? Since every meaningful combination of buttons is equivalent to some explicit, closed-form [13], mathematical formula, above thought experiment becomes practical formulation of the constant recognition problem. In real world implementation, human and calculator are replaced by a computer program222Alternatively, database of the pre-computed results [10]., increasing speed billion-fold. Mathematically, results which can be explicitly computed using hand-held scientific calculator are members of so-called (EL) numbers [14]. Above formulation has many advantages over set-theoretic or decimal string matching approach. It is well-posed, tractable, and can be either answered:
a) positively,
b) in terms of probability, or
c) falsified.

This depend on both maximum length of button sequence (code length, Kolmogorov complexity), as well as precision of the numerical result. One could anticipate, that for high-precision numbers (arbitrary precision in particular) and short code length identification must be unambiguous. For intermediate case (double precision, unspecified code length) we expect to provide some probability measure for identification. For low-accuracy numbers (e.g. of experimental or Monte Carlo origin) failure seems inevitable. Article is devoted to quantify above considerations. Certain statistical criteria are presented for practical number identification. Hopefully, this will convince researchers, that constant recognition problem can be solved practically. Instead of software used mainly for recreational mathematics and random guessing, it could become reliable tool for numerical analysis and computer algebra systems.

Article is organized as follows. In Sect. 2 we discuss simplified computer language (in replacement for physical scientific calculator) used to define Kolmogorov complexity of the formulas. In Sect. 3 we discuss some practical issues regarding numerical implementation of the sequential formula generators. Then, we combine experience/knowledge of Sect. 2 and 3 to find properties and typical behavior of the constructed EL numbers. Sections 4-6 propose three criteria for number recognition: (1) instant error drop-off to machine epsilon instead statistically expected -folding; exponential growth of number of formulas indicates failure of search (2) maximum likelihood formula in the view of statistical process (3) compression ratio of the decimal constant in terms of RPN calculator code.

2 Elementary complex numbers

Our first task is to precisely define set of formulas/numbers we want to identify using decimal expansion. We restrict ourselves to explicit formulas [14] i.e. those computable using hand-held scientific calculator. We assume root-finding procedure is absent. Therefore, as notable exception, most of algebraic numbers are not in above class, because polynomials of the order 5 and larger are not explicitly solvable in general. To be specific, we consider any complex number created from:

  1. integers:

  2. rationals:

  3. mathematical constants:

  4. addition/subtraction, multiplication/division, exponentiation/logarithm (arbitrary base)

  5. elementary functions of one variable:

    • trigonometric:

    • inverse trigonometric (cyclometric):

    • hyperbolic:

    • usually unnamed functions:

  6. function composition

The choice of complex instead of real field is justified by mathematical simplicity. In practical numerical implementation it might be sub-optimal, due to reduced numerical performance and still missing implementation of some elementary functions (e.g: clog2) in standard libraries [15]. Last item in the above list, function composition, is important for completeness. E.g: is well-defined elementary expression, although rarely used in practice. No special functions are used.

Particularly troublesome are sets of symbols related to integers and rationals, as they both are potentially infinite. Without proper handling enumeration quickly get stuck on rational approximations with large integer numerator/denominator. Therefore, we must restrict to some ”small” subset of them. Standard computer languages (C, C++, Fortran) use ”small” integers still too large from our point of view, e.g. 8-bit signed and unsigned ones. Complexity of numbers 0, , 188, …is identical using such an approach. The other extreme follows reductionist definition [17], staring from e.g: , and constructing all other integers as follows:

In calculators, ten digits are present: 0,1,2,3,4,5,6,7,8,9 and integers are entered in sequence using standard positional numeral system. RIES [9] follow similar way, except for 0 and 1. Use of positional base-10 numerals is simple reflection of human anatomy/history, and difficult to include in our virtual calculator. Therefore we restrict to really small magnitude integers, in particular. All other integers must be constructed from above four. Noteworthy, these integers are also among essential mathematical constants included in famous [18] Euler formula:

(1)

Two is also smallest possible integer base for logarithm . Binary system is now used in virtually all modern computer hardware.

The same reasoning applies to rationals. One might enumerate them using Cantor diagonal method, Stern-Brocot tree, or generate unambiguously by repeated composition333See also Appendix B in Supplemental Material. of a function:

starting from zero:

Usually the best is to leave language without explicit rationals, and let them appear by division of integers. Alternatively, addition of reciprocal and successor functions allow for construction of continued fractions, see Appendix B in Supplemental Material.

Set of (named) mathematical constant certainly must include , , and due to (1), as long as we consider complex numbers. Besides this, some mathematicians consider other constants [19] as ”important”, like first square roots (), golden ratio [9], , and so on. Noteworthy, all of them are itself in class, i.e., we can generate them from integers and elementary functions. They do not define any new numbers. Their inclusion only alter language definition and therefore Kolmogorov complexity. Other situation is with inclusion of mathematical constants of unknown status/memberships (unknown to be rational or transcendental), like Euler gamma or Glaisher constant . They may, or may not, extend class. This is itself an application for constant recognition.

Commutative operations like addition and multiplication are in general -ary. To simulate calculator behavior we must treat them as repeated binary operations:

Elementary functions are not independent, e.g:

Scientific calculators sometimes distinguish less important, secondary functions which require pressing two buttons. But mathematically they are completely equivalent. You can compute using , or vice versa. This is especially true in complex domain, where all elementary functions are reducible to and .

Reduction is possible for binary operations as well, e.g:

Replacing multiplication by logarithms and addition is an achievement of medieval mathematics [20] used without changes for 300 years. It was replaced with slide rule used in XX century. They both went extinct with modern computers. However, we point out, that using logarithms and addition you can not only do bottom-up translation to multiplication/exponentiation. You might go opposite way as well, from high-rank (grade) hyper-operation down do additions. Expressing addition/multiplication by exponentiation/logarithms only (up-bottom) is therefore possible, but tricky, see Appendix C in Supplemental Material.

Taking above considerations into account, important question arises: how many constants, functions and binary operations are required for our virtual calculator to be still fully operational? How far reduction process can go, without impairing our computational abilities? This is also known as ,,broken calculator problem” [9]. Another related question is, if such a reduced calculator is optimal for the task of constant recognition. Obviously, maximally reduced button set makes theoretical and statistical analysis convenient, therefore it is very useful as mathematical model. In practice, as we will show later, Kolmogorov complexity of formulas in maximally reduced language, which are considered as simple by humans, might become surprisingly large. On the contrary, formulas simple in reduced operation set, like nested power-towers, look inhuman, and are out of scope traditional mathematical aesthetics, despite small Kolmogorov complexity.

The simplest possible language we found has a length three (). It is still able to perform all operations of the scientific calculator. It includes:

either (2a)
binary exponentiation (2b)
arbitrary base (two-argument) logarithm (2c)

Another very simple language of length four is:

any constant (3a)
natural exponential function (3b)
natural logarithm (3c)
subtraction (3d)

Detailed proofs by exhaustion are presented in Appendices C,D in Supplemental Material. Above two examples clearly justify name for class of explicit elementary numbers discussed in this section. We were unable to find any shorter languages. At least one (noncommutative?) binary operation seem required to start abstract syntax (AST) tree growth (see however Appendix B), and at least one constant to terminate leafs. One of the essential constants or probably also is required, either explicitly, like in (2a), or hidden in function/operation definition (3b,3c). But we cannot prove that further reduction to 2 buttons is impossible. We also cannot prove, that above language of length three (2

) is unique. Therefore, one cannot estimate Kolmogorov complexity of EL number unambigously. Existence of both smallest and unique language, generating all

numbers, would provide natural enumeration of mathematical formulas, similar to Peano arithmetic for natural numbers. Anyway, (2) or (3) are the sets of irreducible operations, for now. Calculator/language can be extended, but under no circumstances any of the buttons defined by (2) or (3) can be removed. Failure to abide by above requirement might result in catastrophic failure of the number recognition software even in the simplest test cases. This is plague of existing implementations, letting unaware end-users to choose base building blocks on their own. This leaves impression of random failures without any obvious reason. Users can use their own set of constants, functions and binary operations, but they must be merged with irreducible set (2) or (3).

It is illustrative to compare (2) and (3) with Mathematica core language [7], composed of:

addition (Plus) (4a)
multiplication (Times) (4b)
exponentiation (Power) (4c)
natural logarithm (Log) (4d)
(4e)

augmented with arbitrarily large integers, including their complex combinations (Complex) and rationals (Rational). Visual comparison of (2), (3) and (4) is presented in Fig. 1.

Fig. 1: Illustration of base constants, functions and binary operations used to define number class. Irreducible set (2) is labeled ”up-botton”, because it is based on highest rank complex (hyper-)operation currently implemented (exponentiation), while (3) (bottom-up) starts with low-rank arithmetic (subtraction). Mathematica core language is included for comparison. All three sets allow generation of any explicit elementary transcendental constant in the EL class. Question mark symbolize unknown ”holy-grail” complex function, which could, if exists, by repeated composition, generate all numbers in unique order.

Accounting for the above considerations, for further analysis we have selected four base sets (buttons) to attack the inverse calculator problem:

  1. E, POW, LOG [ buttons, eq. (2) ]

    (5a)
  2. X, EXP, LN, - [ buttons, eq. (3)]

    (5b)
  3. Pi, E, I, -1, 2, 1/2, Log, Plus, Times, Power [ buttons, eq. (4) ]

    (5c)
  4. full scale scientific calc [ buttons, see below]

    (5d)

Calculator 1 (up-bottom), based on (2) is the simplest found. It use one constant (, PUSH on stack operation), no functions, and only two non-commutative binary operations: arbitrary base logarithm and exponentiation . Noteworthy, AST is a binary tree in case of (5a).

Calculator 2 (bottom-up), based on (3) is the second shortest, and have remarkable property: it can use any numerical constant/symbol to generate all other numbers. For example, , , and so on. We can use any of , , , , , …This property can be exploited in two applications: (1) use our calculator to operate on vectors, (2) use with large Kolmogorov complexity to ”shift” formula generator. The former can be used e.g. to take advantage of modern CPU AVX extensions, the latter to use as some form of random seed for enumeration procedure. This can be done extending any of calculators (5) with as well, but only for (3) it is visible explicite in definition. In presentation of results we used .

Calculator 3 is designed to mimic Mathematica behavior. Constants beyond and were chosen as follows. To enable rapid integer generation via addition, multiplication and exponentiation we use -1 and 2. Incidentally, these numbers are also are initial values for Lucas numbers. Rational numbers are easily generated via continued fractions, thanks to constant, allowing for construction of reciprocal . Rational constant together with form trigonometric functions and square roots . Number of instructions equal to ten is selected intentionally. It will allow for instant estimate of the compression factor in Sect. 6. Both target numerical constant and RPN calculator code can be expressed as string of base-10 digits 0123456789, see next Sect. 3.

Fourth calculator is the largest one. Maximum number of 36 buttons is limited by a current implementation of the fast itoa function used to convert string variables into base-36 numbers, i.e. alphanumeric lowercase digits and letters. It is the most close to what people usually expect from scientific hand-held device. Full list of buttons:

  • constants:

  • functions:

  • binary operations: .

”Full” calculator defined above use 13 constants, 18 functions of one variable, and 5 binary operations.

3 Efficient formula enumeration

Once calculator buttons, equivalent to some truncated computer language specification, are fixed, we face the next task: enumeration of all possible formulas. We exclude recursive implementation. While short and elegant, it quickly consumes available memory, and is hard to parallelize [9]. For sequential generation, we notice that formulas are equivalent to set of all abstract syntax trees (AST) or valid reverse polish [21] notation (RPN, [22]) calculator codes. The former description has an advantage in case where efficient tree enumeration algorithm exists [16], e.g: binary trees [23]. Unfortunately, this is applicable only to the simplest calculator (5a). Other mix binary and unary trees. Therefore, to handle variety of possible calculators, including future extension to genuine -ary special functions (e.g: hypergeometric, Meijer G, Painleve transcendents) our method of choice is enumeration of RPN codes. This has two major advantages. First, enumeration is trivially provided by itoa function with base- numbers (including leading zeros), where is number of buttons. Second, while majority of enumerated codes is invalid, checking RPN syntax is very fast, almost negligible compared to itoa itself, let alone to computation of complex exponential and logarithmic functions. Procedure has two loops: outer for code length , inner enumerating codes of length . Inner loop is trivial to parallelize using OpenMP directives without effort, scaling linearly with number of physical cores. Moreover, hyper-threading is utilized as well, although scaling is only at quarter of core count. Therefore, prospects for high-utilization of modern high-end multi-core CPU’s [2] are looking good. Digits are associated with RPN calculator buttons. For simplest case (2) they are ternary base digits 012, which were assigned to three RPN calculator buttons as: 0 E (), 1 LOG (), 2 POW (). Detailed example of the algorithm for (2) is presented in Appendix A of Supplemental Material.

Combinatorial growth of the enumeration is characterized, in addition to Kolmogorov complexity (RPN code length) , by a total number of possible codes tested so far . It grows with as:

Total number of syntactically correct codes and total number of unique numbers are additonal useful characterizations of the search depth. Obviously, . For perfectly efficient enumeration . Unfortunately, this is currently possible only for rational numbers, see Appendix B. Case is equivalent to knowledge of unique tree enumeration algorithm. To achieve ( in Fig. 2, dashed diagonal line) algorithm must magically somehow know in advance all possible mathematical simplifications. Not only those trivial, like but anything mathematically imaginable, e.g, . We doubt this is possible at all. Some ideas, based on solved rational numbers enumeration, are presented in Appendix B. In practice, we achieved for calculators 1-4 (5), respectively. Large ratio for (5a

) is a result of very simple tree structure, in which only every nine-th odd-

RPN codes are valid. Without taking this into account, , i.e., the worst of four. Full calculator perform surprisingly well.

Fig. 2: Number of unique formulas versus number of valid codes . Solid lines show fits, with for calculators 1-4, respectively. Dashed diagonal line shows perfect case .

4 Quasi-convergence of enumerated formulas and constant recognition

Enumerated formulas provide sequence of elementary transcendental (EL) numbers. For (2) they are, in order of increasing Kolmogorov complexity :

(6)

where repeated numbers were omitted. Position within the same level of is unspecified, but follows from presumed enumeration loop and base- digits association with buttons. From sequence (6) we can extract subsequence of progressively better approximations for target number , e.g., using example from the introduction, . Analysis of sequences obtained this way is main goal of current section. We are interested in convergence properties and criteria for termination of sequence.

Let us assume for the moment, that we are able to provide on demand as many decimal digits for target number

as required. There are two mutually exclusive possibilities: (i) number is in class, defined e.g. by (2) or (3), or (ii) is true non-elementary transcendental constant outside EL class. In the former case, we expect that error eventually, at some finite , will drop to an infinitesimally small value, limited only by currently used precision. Search algorithm could then terminate and switch to some high-precision or symbolic algebra verification method. Unfortunately, equivalence problem is undecidable in general [24]. One cannot exclude, that both numbers deviate at some far more distant decimal place. In the case of (ii) i.e. truly transcendental number, however, sequence will converge indefinitely, by generation of approximations, with progressively more complicated formulas. Realistic convergence examples, computed using extended precision (long double in C) to prevent round-off errors, are presented in Figure 3. The example target from introduction was selected. Number is obviously in class. But is not explicite listed in any of (5), so it must be entered using appropriate sequence of buttons. For RPN calculator (5d) it is obvioulsy 2, SQRT, 3, SQRT, POW with . In Fig. 3 is is visible at , where decimal part show inner loop progress.

Fig. 3: Typical behavior of absolute error for subsequence of progressively better approximations in terms of formulas generated with ”calculators” (5), as a function of Kolmogorov complexity . For (5b) we have used . For (5c) (blue) and (5d) (black) sudden error drop-off to machine epsilon (zero) is observed for small , marking possible ”exact” formula match. Calculations were done using extended precision (long double in C) while horizontal line mark machine epsilon for IEEE 754-2008 (binary64).

Two types of behavior can be observed in Fig. 3. Vertical absolute error scale is adjusted to range , i.e., long double machine epsilon for this and related Figs. 4, 5. In typical situation absolute error for best approximation decreases exponentially with code length. However, if ”true” formula is encountered, error instantly drops off to limiting value. In example from Fig. 3 it is small multiple of machine epsilon for long double precision () or binary zero.

Criterion directly based on Kolmogorov complexity is inconvenient, if one deals with languages of various size, like our calculators 1-4, eq. (5). Convergence rate also depends on language size , cf. Fig. 3. One could use Kolmogorov complexity corrected for language size, or compile formulas generated in extended calculators (5c) or (5d) down to one of primitive forms given by (2) or (3). The former approach usually underestimate, and the latter heavily overestimates true complexity. Therefore we propose another criterion, independent of the language used to generate formulas. Instead of plotting error as a function complexity, we plot -th best approximation (Fig. 4). Now all curves nearly overlap, and observed lower error limit is .

Above considerations provide the first criterion for constant identification:

Criterion 1:

Identification candidate: if absolute error in subsequence of (6) progressively better approximations in terms of formulas deviates ”significantly” from estimated upper limit (drops to numerical ”zero”/machine epsilon in particular) we can stop search and return formula code, as possible identification candidate.

Failure of search: if absolute error in subsequence of progressively better approximations in terms of formulas follow and reach numerical precision limit, or computational resources are exhausted, search failed.

Candidate code must be then verified using symbolic methods, high-precision numerical confirmation test, and ultimately proved using standard mathematical techniques. Above criterion do not provide any numerical estimate for probability of successful identification. However, we point out, that even in case of possible misidentification, unexpected drop of error to machine epsilon marks stop of the search anyway. This is because finding better approximation would require formula with complexity already above threshold, given by intersection of dotted lines in Fig. 4. Beyond that, number of formulas with identical decimal expansion grows exponentially. This behavior mark search STOP criterion for cases, where ”smoking gun” feature from Fig. 3 was not encountered. In practice this still require a lot of computational resources, beyond capabilities of mid-range PC/laptop. That is why curves in Figs. 3-5 are still far from double epsilon, marked with DBL_EPSILON dotted line.

Fig. 4: Similar to Fig. 3, but has been replaced by number of subsequent best approximations found so far. This allows for direct comparison of languages of different length. Possible identification is marked by error for next approximation significantly below .

Expected decrease of approximation error for next best one is of previous. Therefore, we use -folding name for Criterion 1.

Fig. 5: Similar to Fig. 4, but has been replaced by number of valid RPN codes tested so far. Power-law is observed, with . Upper axis show estimated by power-law fit number of unique formulas.

Remarkable observation is provided by yet another Figure 5. Instead of or we used total number of valid RPN codes tested before encountering next approximation. Power law behavior is found, and band of data points (Fig. 5) can be roughly approximated as proportional to . Therefore, obtaining definite negative answer for constant known with machine precision of in terms of Criterion 1 require testing of codes. For double precision machine epsilon it is above . You need either hundreds of CPU cores, or a lot of patience (days of search). Positive identification can be much faster, of course, like for calculators 3, 4 in Figures 3, 4.

5 Statistical properties of the numbers

Observation from Fig. 5 and results of previous section lead to another view of constant identification problem. It can be described as a random process. Consider following numerical Monte Carlo experiment. We generate pseudo-random numbers

from exponential distribution with some scale

and Probability Distribution Function (PDF):

(7)

Next, we compare with our target number , and generate sequence of progressively better approximations. Surprisingly, observed behavior is similar to Fig. 4, but with larger fluctuations. Let’s further assume is known with numerical precision . In other words, true number is in the interval . Assuming (7) we can easily obtain probability of random hit into vicinity :

what gives average number of tries:

Above agrees qualitatively with result presented in Fig. 5 if .

However, statistical distribution of numbers is unknown, PDF (7) was chosen intuitively. Statistical properties of EL numbers should not depend on on language used to generate them, at least in the limiting case of large complexity. Numerical evidence, obtained by collecting real numbers generated according to procedure presented in Sect. 3, is presented in Fig. 6. Numbers with were discarded444See Appendix E in Supplemental Material for distribution of EL numbers on the complex plane. to simplify analysis and presentation.

Fig. 6:

Histogram of first 20k real numbers generated with use of calculators 1-4. Dashed line show exponentially transformed Cauchy distribution (

8).

None of the widespread distributions match numbers presented in Fig. 6. The most close are: Pareto, Levy and Cauchy distributions. Noteworthy, all of them have undefined mean and infinite dispersion. If we look at statistics of base 10 logarithm, it mimic t-Student, Cauchy and Laplace (double ) distributions. Probably this is indeed a brand new statistics, following some empirical distribution. True distribution of real EL numbers (Fig. 6) is, however, quite well described by transformed Cauchy distribution with PDF:

(8)

Therefore, later, statistical distribution of (real) EL numbers will be approximated by Cauchy-like distribution (8).

Let us assume statistical distribution of real EL constants has PDF empirically derived from histogram in Fig. 6 or is given by (8). Let the target number is know with precision , and statistical distribution of error for is .

might be normal or uniform distribution, for example. Then conditional probability

, that -th tested formula is not random match for target value is:

(9)

Then, replacing and integrating, likelihood for proper identification of target constant with i-th unique value is:

(10)

where , i.e. number of unique values tested so far. To understand (10), it might be approximated by

(11)

assuming and using Maclaurin expansion . Value of depends on error distribution , e.g., for uniform and

for normal distribution.

Likelihood (11) becomes zero for large , where search for is pointless anyway due to Criterion 1 (Sect. 4). The most intriguing application of (10) is when is still large compared with machine epsilon, and has uncertainty of statistical nature. This allow for selection of maximum likelihood formula(s) in marginal sense, i.e. for . Because calculation of requires storage of all previously computed values, it is reasonable to replace it with using power-law fit from Fig. 2:

(12)

Noteworthy, (12) as a function of , i.e. number of tested syntactically correct RPN codes, must have a maximum. Formula for with maximum probability is well defined. For limited precision maximum is very pronounced (Fig. 7), but for very small we might not be able to reach it at all, due to limited computational resources (e.g: Fig. 7, solid blue). Equations (10-12) can be understood as follows. If we find remarkably simple and elegant formula at the very beginning of the search, which is however quite far from target , e.g., several

’s in Gaussian distribution, we reject it on the basis of improbable error in measurement/calculation. Similarly, if we find formula well within error range, e.g.

, but after search covering billions of formulas, we reject it expecting random coincidence. Somewhere in the middle lays optimal formula, with maximum likelihood, estimated with use of (10). In practice, likelihood is very small, and -likelihood is more convenient:

Expanding and assuming is Gaussian, we have following useful approximation for -likelihood:

(13)

We note, that probability of the identification related to (10) could be, in principle, estimated directly, especially for single precision floats, as there are only of them. Using floats as bins for sequentially generated EL numbers, we can fill them with numbers associated with their complexity.

Fig. 7: Log-likelihood of subsequent best approximations, as estimated from (13). Red, green, blue and black lines are for calculators 1,2,3 and 4, respectively. Target values are approximations to our introduction example in form with: (dot-dashed, x), (dotted, o), (dashed, +), (solid, #).

Discussion above provides second criterion for constant identification:

Criterion 2:

Identification candidate: if likelihood given by (10) has reached maximum, formula has the highest probability, and should be returned.

One might return a few highest likelihood formulas near maximum as well, or one with largest value so far, if computational resources are limited. Noteworthy, likelihood value provide also relative quantitative estimate of identification probability. We may estimate likelihood (13), using calculator (5c), for (Fig. 7, solid blue line) to be identified as

compared to

for . From our Criterion 3, the former is more than times more probable than the latter. Without this sort of ”ranking”, discussion, if e.g. is ”simpler” or ”more elegant” compared to, e.g., might continue indefinitely, leading to nothing. Likelihood (7) or (13) provide numerical form of the Occam’s razor. It can be applied automatically, within software (CAS) environment.

6 Constant recognition as data compression

Sequence of progressively better approximations in form of RPN calculator codes can be viewed as a form of lossy compression. If exact formula is found, then compression becomes lossless. In the intermediate case, compression ratio provides measure, how good is some formula to recover decimal expansion. This process is illustrated in Table I. Calculator 3 defined in eq. (5c) has been used, because both decimal expansion of the target number and RPN code are strings of the same base-10 digits, i.e. 0123456789. Therefore, compression ratio is simply a number of correct digits divided by code length .

Numerical Code Compression Formula
value ratio
3.141592653589793238512 0 0.00
2.718281828459045235428 1 0.00
0.000000000000000000000 2 0.00
2.000000000000000000000 7 0.00
1.718281828459045235428 164 0.33
1.772453850905516027310 809 0.33
1.648721270700128146893 819 0.33
1.837877066409345483606 0043 0.50
1.837877066409345483606 7053 0.50
1.820796326794896619256 08485 0.60
1.820796326794896619256 80485 0.60
1.820796326794896619256 80845 0.60
1.820796326794896619256 88045 0.60
1.821126701185962651818 0338975 0.43
1.821126701185962651818 7033895 0.43
1.821126701185962651818 8303975 0.43
1.824360635350064073446 8819745 0.43
1.824360635350064073446 8781945 0.43
1.824360635350064073446 8197485 0.43
1.824360635350064073446 7819485 0.43
1.821126701185962651818 7830395 0.43
1.821662858741926632288 6091579 0.43
1.822361069544464599575 2298979 0.57
1.822413909696397869321 77408934 0.50
1.822722133555469366033 80790539 0.50
1.822690334737686312645 004377539 0.56
1.822634654966242214488 888854979 2.11
TABLE I:

In Table I, compression ratio of the target number in form of RPN calculator codes (5c) is shown.

Buttons/operations were assigned as follows:
0 Pi,
1 E,
2 I,
3 Log,
4 Plus,
5 Times,
6 -1,
7 2,
8 1/2,
9 Power.

In fact, last code in Table I 888854979, with pronounced compression ratio of , is equivalent to exact formula, with RPN sequence ½, ½, ½, ½, , , Power, 2, Power, i.e., .

In general case, compresion ratio is given by:

(14)

where is absolute precision of the approximation, - RPN code length, - number of calculator buttons.

Now we can formulate the last criterion for constant recognition:

Criterion 3:

Identification candidate: if compression ratio given by (14) is , or reaches maximum in the course of search, then formula should be returned.

Failure: if formula/code is unlikely match.

Criterion 3 has an advantage of being very simple. It do not require recorded history of search, like Criterion 1 (Sect. 4), of searching for maximum and knowledge of statistics, like Criterion 2 (Sect. 5). If indeed has a maximum, then it strenghten identification, but this is not required. You might ask for compresion ratio of any combination of decimal expansion and formula. The only required action is to compile formula to RPN code (5) or similar one. Therefore, it will work with any searching method, e.g: Monte Carlo, genetic or shortest path tree algorithms. It is weak compared to -folding or statistical criteria, but could easily exclude most of formulae produced in variety of software, especially very complicated ones, and those including large () integers.

7 Conclusions

Three criteria proposed in the article provide robust tool for decimal constant identification.

Fig. 8: Recognition indicators as a function of . Solid red line shows absolute -folding , while dotted red line relative -folding (upper panel). Green line show likelihood (13) (middle panel), and blue compression ratio (14) (bottom panel).

Process of recognition, applying all three criteria, is illustrated in Fig. 8, using blind-test target value of . Assuming all digits are correct, we adopted (gaussian) error . Using calculator (5c), we plot three indicators of matching quality. For -folding measure we use

and additionally:

where -th best approximation error is . For likelihood, we used (13), and for compression ratio (14). Clearly, three subsequent data points, for , stand out (Fig. 8).

From first criterion perspective (Sect. 4), error dropped nearly times compared to expected (Fig. 8, upper panel, solid red). Is also was times smaller compared to previously found , while statistically anticipated decrease was (Fig. 8, upper panel, dotted red). Using criterion 2 (Sect. 5) we found likelihood many many orders of magnitude larger compared to any other formulas. Three data points with nearly equal (Fig. 8, middle panel, green crosses), 16,17,18-th best approximations, are in fact the same formula typed using different RPN sequences. Differences are due to round-off errors. Shortest one is preferred in terms of criterion 3 (Sect. 6), see Fig. 8, lower panel. Constant is then unambiguously recognized by all three criteria as , with no other candidates in sight.

Above example show possible way to precise formulation of decimal constant recognition problem, as a task reverse to pushing button sequence, and its solution. Using three proposed criteria we are able to judge which formulas, matching given floating-point constant, are the most likely. In original formulation they require enumeration of all possible codes with growing complexity, but statistical and compression criteria are in fact independent of the method used to obtain expression. They require only formula to be compiled into RPN code. Then, Criterion 3 (14) can be applied directly, and Criterion 2 (13) indirectly, using code length to find upper limit on (Fig. 2).

Calculator used to enumerate codes can be arbitrary, as long as it includes ”irreducible” buttons equivalent to (2) or (3

). However, search results (required depth in particular) depend on calculator definition. This is especially visible in search of human-provided test cases, with strong bias towards decimal numerals, and repulsive reaction to power-towers and nested (exponential) function compositions. Unfortunately, the latter are strong at approximation, with multi-layer neutral networks being notable example using sigmoidal functions. Therefore, for numbers provided by humans, full calculator (

5d) is usually the best, for mathematical results (5c) and for randomly generated expressions (5a) or (5b).

Actually, while searching for constants, we are implicitly dealing with enumeration of all EL functions of one complex variable. This is especially visible in definition of calculator (2). Constant functions are just subset of them. Therefore, mathematical and statistical tools developed to solve constant identification problem, could be extended to identification of functions of one variable. This is a task for further research, with many more potential applications.

References

  • [1] David H. Bailey and Borwein Jonathan M. Exploratory Experimentation and Computation. Notices of the AMS, 58:1410–1419, November 2011.
  • [2] AMD EPYC™ 7002 Series Processors and DELL POWEREDGE™ R6525 Servers Set World Record on Industry Standard Decision Support Benchmark https://www.amd.com/system/files/documents/tpc-h-dell-3tb-exasol.pdf
  • [3]

    Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. Science, 324(5923):81–85, 2009

  • [4] Sohrab Towfighi, Symbolic regression by uniform random global search, Neural and Evolutionary Computing (cs.NE) Cite as: arXiv:1906.07848 [cs.NE]
  • [5] David H. Bailey and Simon Plouffe, Recognizing Numerical Constants, http://www.cecm.sfu.ca/organics/papers/bailey/paper/html/paper.html
  • [6] Maple (2017.3). Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario.
  • [7] Wolfram Research, Inc., Mathematica, 12.0.0 for Microsoft Windows (64-bit) (April 6, 2019), Champaign, IL (2019).
  • [8] Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103
  • [9] Robert Munafo, RIES - find algebraic equations, given their solution. https://mrob.com/pub/ries/index.html Accessed: 2019-11-25.
  • [10] Inverse Symbolic Calculator, http://wayback.cecm.sfu.ca/projects/ISC/ISCmain.html
  • [11] D. Hough, ”The IEEE Standard 754: One for the History Books” in Computer, vol. 52, no. 12, pp. 109-112, 2019. doi: 10.1109/MC.2019.2926614
  • [12] Ray J. Solomonof, INFORMATION AND CONTROL. Volume 7, No. 2, June 1964, pp. 224–254. Copyright by Academic Press Inc.
  • [13] J. M. Borwein and R. E. Crandall. Closed forms: What they are and why we care. Notices of the American Mathematical Society, 60:50–65, 2013.
  • [14] Timothy Y. Chow. What is a closed-form number? The American Mathematical Monthly, 106(5):440–448, 1999.
  • [15] Sandra Loosemore with Richard M. Stallman, Roland McGrath, Andrew Oram, and Ulrich Drepper, The GNU C Library Reference Manual for version 2.31, https://www.gnu.org/software/libc/manual/
  • [16] Donald Knuth, The Art of Computer Programming, Volume 4, Fascicle 2: Generating All Tuples and Permutations, Addison-Wesley Professional; 1 edition (February 24, 2005), ISBN-13: 978-0201853933, ISBN-10: 0201853930
  • [17] Alfred Tarski. A decision method for elementary algebra and geometry. U. S. Air Force Project Rand, R-109. Prepared for publication by J. C. C. McKinsey. Litho-printed. The Rand Corporation, Santa Monica, California, 1948
  • [18] Eli Maor, The story of a number , Princeton University Press, 1994, ISBN 0-691-03390-0
  • [19] Eric W. Weisstein, Constant, From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/Constant.html
  • [20] Denis Roegel. A reconstruction of the tables of Briggs’ Arithmetica logarithmica (1624). [Research Report] 2010. ffinria-00543939f
  • [21] Jan Lukasiewicz (1957). Aristotle’s Syllogistic from the Standpoint of Modern Formal Logic. Oxford University Press. (Reprinted by Garland Publishing in 1987. ISBN 0-8240-6924-2)
  • [22] C. L. Hamblin, Translation to and from Polish Notation, The Computer Journal, Volume 5, Issue 3, November 1962, Pages 210–213, https://doi.org/10.1093/comjnl/5.3.210
  • [23] L. Tychonievich L (2013) Enumerating Trees. URL https://www.cs.virginia.edu/ lat7h/blog/posts/434.html
  • [24] Richardson, Daniel (1968). ”Some undecidable problems involving elementary functions of a real variable”. Journal of Symbolic Logic. 33 (4). Association for Symbolic Logic. pp. 514–520. doi:10.2307/2271358. JSTOR 2271358.

Appendix A Enumeration example

Detailed example of formula enumeration algorithm. Top-down base system (see main text) has been used for simplicity. Ternary base digits were assigned to three RPN calculator buttons as: 0 E (), 1 LOG (), 2 POW (). After code length 3 invalid codes were omitted to save space.

Enum CODE Syntax RPN sequence formula
0 0 VALID E
1 1 INVALID
2 2 INVALID
3 00 INVALID
4 10 INVALID
5 20 INVALID
6 01 INVALID
7 11 INVALID
8 21 INVALID
9 02 INVALID
10 12 INVALID
11 22 INVALID
12 000 INVALID
13 100 INVALID
14 200 INVALID
15 010 INVALID
16 110 INVALID
17 210 INVALID
18 020 INVALID
19 120 INVALID
20 220 INVALID
21 001 VALID E, E, LOG
22 101 INVALID
23 201 INVALID
24 011 INVALID
25 111 INVALID
26 211 INVALID
27 021 INVALID
28 121 INVALID
29 221 INVALID
30 002 VALID E, E, POWER
31 102 INVALID
32 202 INVALID
33 012 INVALID
34 112 INVALID
35 212 INVALID
36 022 INVALID
37 122 INVALID
38 222 INVALID
210 00101 VALID E, E, LOG, E, LOG
219 00201 VALID E, E, POWER, E, LOG
228 00011 VALID E, E, E, LOG, LOG
255 00021 VALID E, E, E, POWER, LOG
291 00102 VALID E, E, LOG, E, POWER
300 00202 VALID E, E, POWER, E, POWER
309 00012 VALID E, E, E, LOG, POWER
336 00022 VALID E, E, E, POWER, POWER

Appendix B Enumeration of integer and rational numbers

If we restrict ourselves to simplest case of integer and rational numbers, enumeration procedure without repetition, i.e., one-to-one mapping of non-negative integers into integers is known:

Positive rationals can be enumerated by repeated composition of the function:

starting with zero.

Surprisingly, function is composition of two self-inverse functions:

Fig. 9: Self-inverse function, able to generate all integers and rationals without repetitions.

Illustrative example is provided as follows. Let us define two additional self-inverse functions:

(15)

and

(16)

Any integer and rational (including negative) can be obtained be repeated composition of functions and , cf. Fig. 9, starting with number (or a function) zero. Moreover, the appear in non-repetitive order:

Above properties are remarkable, and suggest possible way to non-repetitive generation of numbers by function composition, in unique order. However, it is unclear what kind of function(s) is could be. For example, self-inverse function related to exponentiation is:

where -1 can be replaced by any other constant. In fact, it could be generalized to:

with arbitrary . So far, our attempts to find failed. We only guess they are somehow related to and possibly .

Appendix C Proof of completeness of up-bottom base set

Goal of his section is to show, that all explicit elementary complex numbers can be reduced to three elements.

We start with symbols:

One can compute:

Reversing role of logarithm base and argument, we also get:

This way one can compute all natural numbers and their reciprocals (egyptian fractions).

Multiplication can be computed by:

while division is:

Doing addition is tricky, but possible:

Reciprocal is:

and sign change:

This complete basic 6 binary operations. We need only square root:

of -1:

to compute remaining trigonometric functions. In particular (see Fig. 10 for more readable form):

This completes the proof.

Fig. 10: Tree form of computed using primitive calculator CALC1.

Appendix D Proof of completeness of bottom-up base set

Goal of his section is to show, that all explicit elementary complex numbers can be reduced to calulator with four buttons.

We start with symbols:

One can calculate:

Now, we know how to compute:

Addition is:

One can change sign and compute reciprocal with:

So far we have:

succesor is:

Using succesor and reciprocal on can compute all integers and rationals. Multiplication and division using logarithms are well-known:

Binary exponentiation and logarithm are:

Let’s proceed to square root and :

Number is:

Now, calculating trigonometric functions is straightforward:

Above shows, that all constants, functions and binary operations from Sect. 2 can be computed using CALC2.

Appendix E Distribution of EL numbers on complex plane

Distrubution of the EL numbers generated by sequence on complex plane is presented in Figs. 11 and 12. Visually, it is far from random. However, it is not fractal, a because EL numbers include rationals, which are everywhere dense.

Fig. 11: Distribution of the EL numbers on complex plane for calculators 1-4. Calculator 1 (top-down) is shown in the upper-left panel. Only formulas with are shown, up to Kolmogorov complexity , what gives 1500k unique complex numbers. Upper-right panel is for CALC2 (down-top), up to ; lower-left shows CALC3 (Mathematica) for ; lower-right present full calculator CALC4 up to .
Fig. 12: Same as in Fig. 11, but his time distribution of the fractional part of the EL numbers on complex plane is shown.