Unsupervised Discovery of Formulas for Mathematical Constants

Authors: michael  shalyt
uri  seligmann
itay  beit
halachmi  ofir
david rotem elimelech
ido  kaminer
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:29:04

Abstract: Ongoing efforts that span over decades show a rise of AI methods for accelerating scientific discovery, yet accelerating discovery in mathematics remains a persistent challenge for AI.
Specifically, AI methods were not effective in creation of formulas for mathematical constants because each such formula must be correct for infinite digits of precision, with 'near-true' formulas providing no insight toward the correct ones. Consequently, formula discovery lacks a clear distance metric needed to guide automated discovery in this realm.

In this work, we propose a systematic methodology for categorization, characterization, and pattern identification of such formulas. The key to our methodology is introducing metrics based on the convergence dynamics of the formulas, rather than on the numerical value of the formula. These metrics enable the first automated clustering of mathematical formulas.
We demonstrate this methodology on Polynomial Continued Fraction formulas, which are ubiquitous in their intrinsic connections to mathematical constants, and generalize many mathematical functions and structures.
We test our methodology on a set of 1,768,900 such formulas, identifying many known formulas for mathematical constants, and discover previously unknown formulas for $\pi$, $\ln(2)$, Gauss', and Lemniscate's constants. The uncovered patterns enable a direct generalization of individual formulas to infinite families, unveiling rich mathematical structures. 
This success paves the way towards a generative model that creates formulas fulfilling specified mathematical properties, accelerating the rate of discovery of useful formulas.

Main paper: 
 expression identification through pslq analysis: 
 The PCF limit is computed and a targeted PSLQ search is used to identify an expression that matches the resulting digits and incorporates the constant C. 
introduction: 
 Historically, formulas of mathematical constants were a symbol of aesthetics and beauty. Continued fraction formulas such as those for the Golden Ratio ϕ and for tan(x)EQUATIONenable calculating infinitely many digits for these constants. Discovering such formulas often leads to profound revelations regarding the properties and underlying structure of fundamental constants. For example, the continued fraction formula for tan(x), shown in Eq. 1, was used by Johann Heinrich Lambert in the first proof of the irrationality of Pi [Lambert, 1768 ][Berggren et al., 2004] . Unfortunately, such formulas are notoriously hard to find on-demand, often relying on a mathematician's profound intuition. Part of the challenge is the lack of a well-defined 'distance', or a metric, between a formula and a given constant. i.e., there is no known way to tell whether a formula is nearly accurate. The formula either works, or it does not. In other fields of science, a prediction accurate to 1000 digits is precise enough for any practical need. However, in mathematics, if the 1001 st digit is wrong, the formula is incorrect and gives no insight regarding a correct formula. This lack of a metric is a substantial hurdle both for human efforts and for automated analysis, as many methods for optimization such as gradient descent become unsuitable.Recent efforts developed computer algorithms to discover a multitude of formula hypotheses for mathematical constants [Raayoni et al., 2021] , even implementing the first large-scale distributed computation for such discoveries [Elimelech et al., 2023] , but they relied mostly on exhaustive search methods. These approaches complements earlier applications of algorithms for automated theorem proving (ATP) (such as computer proofs of hypergeometric identities [Petkovsek et al., 1996] , Malarea [Urban, 2007] , and Flyspeck [Kaliszyk and Urban, 2012] ), and automated conjecture generation (ACG) (such as mechanical mathematics [Wang, 1960] , the Automated Mathematician [Lenat, 1982] , EURISKO [Lenat and Brown, 1983, Davis and Lenat, 1982] , and Graffiti [Fajtlowicz, 1988] ).Here we propose a fundamentally new methodology for automated investigation and discovery of formulas for mathematical constants. We constructed a large dataset of continued fractions, and enriched it with metrics based on their convergence dynamics, which are found to embody fundamental information about each continued fraction. These dynamical metrics enable the identification and generalization of patterns within the dataset. Using the metrics, we develop a process of categorization and clustering (Fig. 1 ) of continued fractions that share similar values of their dynamical metrics.Analyzing each automatically identified cluster of formulas, we find that all its members often relate to the same mathematical constant, showing the value of the dynamical metrics for the discovery of new formulas and the internal structure of families of such formulas. This novel method of formula discovery allowed us to identify both previously known and completely new formulas for constants such as π, ln(2), cot(1), the Golden Ratio, square roots of multiple integers, the Gauss constant, and the Lemniscate constant. ) If all anchors are in agreement, choose a random subset of other points in the cluster and use PSLQ to validate that they also relate to the same constant. If the validation is successful, the cluster is labeled. If not, the cluster should be split. (d.2) If the anchors relate to different constants, the cluster should be split -return to step c for finer clustering of the data. When focusing on a specific cluster, the best metrics could be different than those for the full dataset. (e) If no anchor is found in a certain cluster, attempt to label by (e.1) choosing a small subset of PCFs in the cluster and running a PSLQ search for each of them against a large set of potential constants. If a connection is found, the cluster now has an anchor -return to step d. (e.2) If an anchor is still not found, attempt to connect a sample of data points within the cluster using PSLQ. If successful, conclude that the cluster is correct, but has no identified constant. Define a new label for that cluster. If PSLQ failed to connect points within the cluster, return to step c for finer clustering.If no further refinement is appropriate, flag the cluster for further analytical investigation.As part of our analysis of metrics of continued fractions, we developed and applied the most complete classification of polynomial continued fractions known to date, detailed in Appendix B. This classification includes the prediction of whether the continued fraction converges directly based on its defining polynomials.Traditional clustering methods attempt to relate data points by calculating distance metrics based on the parameters of these data points, e.g., the coefficients of the defining polynomials. The most common approaches (like SVM) rely on linear classification, while more advanced methods rely on non-linear kernel transformations -but usually use various functions calculated directly on the data parameters. In our dataset, each point is a continued fraction formula defined by the polynomials used to construct it. However, we find that it is not sufficient to use the parameters of the polynomials, and not even the numerical limit of the continued fraction. Instead, we find that it is the dynamics of the continued fraction generated by these polynomials, rather than any direct function on their coefficients, which provides the most useful metrics for analysis. In other words, we find that the useful underlying metrics to extract from each data point are embedded within the intricate progression of the sequence created by the formula, rather than the explicit numerical value (limit) of that formula, or the coefficients defining it. Thus, in order to assess the distance between two polynomial continued fractions, and identify relations between such formulas, it is imperative to characterize the nuanced behavior of their sequences, analyzing trends in the convergence process of these sequences, spanning over numerous terms. Some of the metrics we extract, such as the irrationality measure, are well-known in the mathematical community, yet were never considered for a large-scale classification effort. The evaluation of the irrationality measure is technically challenging for formulas whose limit is not known in advance (which is the vast majority). This challenge made it impossible to extract the irrationality of formulas for a large dataset. Consequently, we develop a new algorithm -the Blind-δ algorithm (Section 3.4)to enable the extraction of the irrationality measure of a continued fraction without prior knowledge of its limit. This algorithm allowed us to extract the irrationality measure for the entire dataset.These advances provide the building blocks for our novel methodology for formula discovery. We cluster formulas by their 'closeness' to other formulas according to these new metrics, thus identifying promising formulas regardless of their numerical value (Fig. 1left ). Once a candidate formula is found, we numerically validate it by calculating its value to a large precision and then identifying its relation to a mathematical constant. The "generate ⇒ validate" approach is inspired by works in AI-driven code generation [Ridnik et al., 2024] and problem solving in geometry [Trinh et al., 2024] .3 Methodology for Data-Driven DiscoveryIn this work we chose to focus on polynomial continued fraction (PCF) formulas as our test case due to the combination of their simplicity and expressive power. PCFs relate to a wide range of mathematical fields, represent a variety of constants, are equivalent to infinite sums [Euler, 1748] , and cover mathematical functions such as Bessel functions, trigonometric functions, integral families, widely used Taylor series, and generalized hypergeometric functions [Cuyt et al., 2008] . Thus, studying PCFs can provide insight into a plethora of mathematical objects and applications.A PCF at depth n is defined as:EQUATIONwhere a n = a(n) and b n = b(n) are evaluations of polynomials with integer coefficients. The PCF value is the limit L = lim n→∞ p n q n (when it exists). The converging sequence of rational numbers pn qn provides an approximation of L, which is known as a Diophantine approximation.While irrational numbers cannot be expressed using a simple quotient of integers, they can be approximated by them. Moreover, some approximations are "better" than others, and one way to evaluate their quality is by a quantity called the irrationality measure [Hardy et al., 1979] .We define the irrationality measure of a sequence pn qn → L as the limit δ n → δ, withδ n = -log L - p n q n log |q n | -1 , qn = q n gcd(p n , q n ) (3)For every L ∈ R, the irrationality measure of L is defined as the supremum of all possible δ for which there is a sequence of distinct rational numbers pn qn → L; pn qn ̸ = L that satisfiesEQUATIONIt is known that for irrational numbers this measure is ≥ 1 (Dirichlet theorem for Diophantine approximations), and for rationals it is 0.Note that the irrationality measure of L is greater or equal to the irrationality measure of any specific sequence converging to the same L. While the irrationality measure of a sequence can be any number ≥ -1, the irrationality measure of its limit L is always either 0 or ≥ 1 [Church, 2019] .After a candidate expression has been identified, the PCF is computed to a greater depth. The resulting values are compared against the proposed expression to assess the accuracy of the match and verify convergence.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We discuss limitations and required assumptions of the Blind-δ Algorithm in section 3.4, the dynamical metric estimation via down-sampled curve fitting in Appendix A, and of PCF convergence and the delta-predictor formula (Eq.5) in Appendix F. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that are not acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community.Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
δ-predictor formula: 
 The classification of a large number of continued fraction formulas requires an efficient and accurate calculation of the irrationality measure δ for each formula. This calculation is challenging because it depends on the asymptotic behavior of the converging sequence, and because δ appears as an exponent of a large basis number. The δ-Predictor formula that we present here provides a way around this challenge -requiring no specific knowledge about the convergence rate and trajectory, or even about the sequence limit itself:EQUATIONwhere λ 1 (n) and λ 2 (n) are the eigenvalues of the matrix0 b n 1 a n , |λ 1 (n)| > |λ 2 (n)|.This formula extends a hypothesis made in a previous work [David et al., 2021] , which was limited to PCFs with deg(B) = 2deg(A) and with a qn that grows exponentially. As we found in this work, Eq.5 works for any converging PCF. It was validated numerically and proven for the deg(B) = 2deg(A) case in Appendix F. This formula helps estimate the irrationality measure, a critical dynamical metric for our work. Specifically, the asymptotic behavior of qn and λ1 /λ2 are still required for finding δ predicted , but they are usually easier to derive.Each PCF formula is defined by the polynomials that generate it. This work focuses on polynomials up to 2nd degree:a n = A 2 n 2 + A 1 n + A 0 , b n = B 2 n 2 + B 1 n + B 0 , with integer coefficients in the domain -5 ≤ A i ≤ 5, -5 ≤ B i ≤ 5.We removed the a = 0 and b = 0 cases, as they break the PCF structure, leaving us with 1,768,900 formulas. Some of these PCFs do not converge to a single limit, rendering their measured metrics meaningless (see Appendix B for the classification method we developed to predict PCF convergence). We filtered out all formulas that do not converge, providing the final filtered dataset of 1,543,926 formulas.The conventional classification of the PCFs is by the coefficients of their polynomials a n , b n , (A 2 , A 1 , A 0 , B 2 , B 1 , B 0 ), and by their numerical limit L, which we evaluate here at depth n = 2000.Going beyond these conventional classification, our methodology relies on dynamics-based metrics calculated for each formula:• The irrationality measure: for each PCF, we calculate δ predicted (Eq.5) and compute δ directly using the Blind-δ algorithm (presented in section 3.4) at depth n = 1000. Fig. 2a presents example δ evaluations.• The convergence rate dynamics, comprised of three parameters: for each PCF, we fit the approximation error ϵ(n) := pn qn -L , which scales as ϵ(n) ∼ n! η • e γn • n β for large n. We store the fitted η (factorial coefficient), γ (exponential coefficient), and β (polynomial coefficient). The process is detailed in Appendix A.• The growth rate of qn , comprised of three parameters: For each PCF, we fit the denominator to qn ∼ n! η ′ • e γ ′ n • n β ′ for large n. We store the fitted parameters η ′ , γ ′ , β ′ .Based on this set of metrics, we applied unsupervised clustering for unlabeled data (using the densitybased OPTICS algorithm [Ankerst et al., 1999] ). The clustering is the key component in the algorithm we developed (Fig. 1 ), leading to the discovery of a variety of formulas and data patterns that exposed formula families (see sections 4.1, 4.2, 4.3, and 4.4 for selected results). 1) and the Silver Ratio. The δ of these constants is known to be 1 (green dashed lines). The blue dots show the numerical convergence of δ (Eq.3) to the correct value. The red dots show the evaluated δ-Predictor formula (Eq.5), following the numerical δ very closely in the Silver Ratio formula, while taking a completely different (and much slower) trajectory in the cot(1) formula; yet both converge to the correct value δ = 1. For the purposes of clustering, δ predicted was evaluated at n = 10 9 , providing an accurate estimation for δ. (b) δ n (n = 1000) vs. the limit value for PCFs in our dataset. While δ values seem to follow a pattern, the limit value distribution does not contain relevant information (the higher density of PCFs near the Y axis arises from the small coefficients of the polynomials in our dataset). Our dataset contains 913, 056 irrationality-proving formulas, most of which are not yet linked to any known constant. (c) Exponential growth coefficients of qn and of ϵ(n) for PCFs with deg(B) = 2deg(A). Note the surprising "band-structures" that this view reveals. A few of the clusters have been identified, but the reason for the appearance of these "bands" and the properties of most clusters remain as open questions. (d) Example PCFs in the dataset that converge to a value close to the constant cot(1) (±10 -5 ) and yet are not related to it, showcasing the challenge of mathematical formula discovery. For visual clarity, error bars not shown. See Appendix A for a discussion regarding measurement errors.The irrationality measure δ of a PCF is of mathematical interest, and (as we will see in section 4) is a powerful dynamical metric. Unfortunately, evaluating δ using Eq.3 requires knowing the series limit L, making its estimation for a large set of unlabeled PCFs impractical.The Blind-δ algorithm was created in order to circumvent this limitation. Instead of inspecting the convergence behavior of pn qn → L, we inspect the convergence behavior of pn qn → pm qm for specific m > n. Given a rational approximation pn qn → L, we approximate the error rate ϵ(n) withp n q n - p m q m = ϵ(n) • 1 - ϵ(m) ϵ(n) . If 0 < s < 1 -ϵ(m) ϵ(n)< S is bounded away from zero and infinity for all n large enough, then the approximation of the Blind-δ algorithm has the same convergence rate as Eq.3, bypassing the need to evaluate L. Intuitively, this condition holds whenever |ϵ(n)| → 0 fast enough, which is true for the vast majority of PCFs (see Appendix F for details).Note that m has to grow with n. In practice, our implementation of the algorithm in this work uses m = 2n, so in order to study δ up to n = 1000, we use m = 2000.As part of the automated formula discovery flow we choose the best metric (for each step), in terms of representation power, which is measured by applying the Davies-Bouldin Index [Davies and Bouldin, 1979] on clustering using each metric individually. Table 1 shows results for a randomly chosen sample of 25K converging PCFs. Note the extremely poor performance of the PCF limit L, in agreement with Fig. 2b, d . This dimensionality reduction is important for efficiency during the clustering step and for better explainability. The former is because the dataset size grows exponentially with the PCF degree and with the magnitude of the polynomial coefficients. Other metrics were tested. Some have been shown to have little to no representation power (e.g. p n and q n , as defined in Eq.2, modulo various primes, their sign, their GCD etc.) while others show potential and are left for future study (e.g. the leading Fourier coefficients of the "noise" around the fit of qn ). A relatively small number of metrics were measured and used, which helped keep the results mathematically explainable. Nevertheless, the clustering using the metrics in Table 1 showcases the strength of our dynamical metrics approach.Once the expression has been verified, an automated analytical proof is attempted via the following steps.Given two polynomials a(n) and b(n) defining a PCF, Euler's formula [Euler, 1748] states that when there exist polynomials h1 (n), h 2 (n), f (n) such that b(n) = -h 1 (n)h 2 (n), f (n)a(n) = f (n -1)h 1 (n) + f (n + 1)h 2 (n + 1), the limit of the PCF equals:f (-1)h 1 (0) f (0) + f (1)h 2 (1) f (0) • ∞ k=0 f (0)f (1) f (k)f (k + 1) k i=1 h 1 (i) h 2 (i + 1) -1 (10)The proof is then completed using known identities of infinite sums.If the PCF's polynomials do not satisfy Euler's formula's conditions, the PCF is flagged for further examination.For example, one may take the PCF defined by the polynomials a n = -2 and b n = 4n 2 + 4n + 1 which was labeled as related to π. The PCF was numerically identified and validated as converging to1 -3 3-3π 4Next, the polynomials h 1 , h 2 and f were identified as 2n + 1, -(2n + 1) and 1 respectively, which corresponds to the sum 1 -3 ∞ k=0 k n=1 -2n-1 2n+3. This infinite sum is known and was indeed proven to converge to 1 -3 3-3π 4 .For the entire list of proven formulas, see table 5 in appendix G.Not all PCFs converge. Clearly, if pn qn does not have a well defined limit, then some of our numerically measured metrics lose their meaning. Though we had algorithmic safeguards to detect such cases and remove them from the analyzed set, it was valuable to identify a pattern and formulate a rule-set that predicts the convergence of a PCF.For that purpose we turned to the matrix representation of a continued fraction to depth n (see Appendix F.1 for details):p N -1 p N q N -1 q N = N -1 n=1 0 b n 1 a n = = N -1 n=1 a n-1 1 0 0 1 a 0 N -1 n=1 0 bn an a n-1 1 1 1 0 0 a N -1Assuming a n ̸ = 0 for n ≥ 0.Analyzing the eigenvalues of the matrices within the matrix product as n → ∞ allows for examining the asymptotic behavior of the continued fraction. Their characteristic polynomial is λ 2 -λ -bn anan-1 , and we propose observing the discriminant of this polynomial -more specifically its dominant power of n:EQUATIONHere we assume C s ̸ = 0 and s is some integer. Based on the data, we compiled table 3 as a summary of the conjectured behavior of any polynomial continued fraction based on s and C s .Table 3 : Summary of PCF behavior characterized by s and C s as defined in Eq.11.ConvergenceC s > 0 C s < 0 ✗ s ≥ 3 s ≥ 0 ✓ s ≤ 2 s ≤ -1We can further elaborate on the converging cases by discussing the conjectured rate of convergence. Usually, a PCF is expected to converge at a sub-exponential rate, but in the case of s = 0, C 0 > 0 it is expected to converge faster:• If C 0 ̸ = 1 then the PCF will converge at an exponential rate, and the exact rate of convergence increases monotonically as C 0 → 1, with a vertical asymptote at C 0 = 1. The convergence rate is identical for C 0 and 1 C0 . • If C 0 = 1 then the PCF will converge at a factorial rate. More specifically, if we find the second most dominant power∆ n = C 0 + Ct n t + O( 1 n t+1) for some C t ̸ = 0 and integer t > 0 then the precision will grow at a rate of O(n! t ).We used these rules (in conjunction with the measurements mentioned in section 3.3) to validate that all PCFs we analyze and cluster do converge and their measured metrics are well defined.Polynomial continued fractions use two polynomials a n = a(n) and b n = b(n) to generate a sequence of rationals p n /q n . However, the same sequences with identical behaviour can be generated using more then one set of polynomials. By identifying transformations under which the dynamics of p n /q n remains invariant, we can formally prove equivalence between data points, validating the clustering power of the chosen metrics.By rearranging the continued fraction definition, we can see how equivalent a n and b n series can arise:EQUATIONIndeed, by defining a new pair of polynomials a ′ n = a n c n ; b ′ n = b n c n c n-1 we get an equivalent continued fraction which converges to c0pn qn . Clearly, since the resulting sequencep ′ n q ′n is identical to the original one, it exhibits the same dynamics. We call this process "Inflation by c n ".The metrics we are interested in are mostly not affected by a finite number of elements in the sequence. For example, both the convergence rate and δ discuss an overall trend as n grows. Consequently, we can initiate the sequence at different values of n ̸ = 0 without changing the latent parameters. When expressing these transformations as modification to the continued fraction definition, we see that the limit of the continued fraction might change due to this shift in sequence initiation, but only by a rational fractional transform.a 0 + b 1 a 1 + b 2 a 2 + b 3 . . . + b n a n + . . . = p n q n ⇒ a 1 + b 2 a 2 + b 3 . . . + b n a n + . . . = b 1 pn qn -a 0For example, we consider the cluster of formulas related to the golden ratio shown in figure 3b . A large portion of these PCFs stem from transforming the known formula for the golden ratio shown in Eq.1 via the methods aforementioned. The exact transformations are depicted in Table 4 .Table 4 : Continued fractions converging to linear fractional transformations of the Golden Ratio ϕ, found using the top left cluster of Figure 3b . Numerous data points in this cluster exhibit identical sequence dynamics and are equivalent under the inflation and index indentation transformations. The equivalent data points create families of continued fractions in the cluster. Discrepancies between the calculated irrationality measure within the same family is ascribed to numerical inaccuracies, typically on the order 0.001. However, when comparing families, discrepancies in the irrationality measure rise to a magnitude of 0.04, suggesting potential deeper distinctions among these PCFs. EQUATIONThe methodology can be scaled to larger data sets in multiple ways. For example, by extending the range of PCF coefficients from [-5, 5] to [-10, 10] , the size of the data set increases by approximately 50 times. Additionally, by considering polynomials of higher degrees -such as advancing to thirddegree a n and fourth degree b n with coefficients in [-5, 5] -the dataset size can be amplified by approximately 1,333 times.To manage the substantial computational demands associated with these expansions, a three-fold solution is proposed.First, a dynamically chosen measurement depth is implemented during evaluation, aiming for a fixed precision across all PCFs rather than maintaining a constant depth for all computations. This approach optimizes computational efficiency by adjusting the measurement depth based on the specific convergence rate of each PCF.Second, using known equivalences of PCFs such as inflation (see Appendix C), it is possible to substantially decrease the effective size of the dataset that needs to be measured. In particular, when c n = -1, we observe that the sign of a n does not affect the dynamics of the sequence -only flips the sign of the limit to -L. For every PCF its inflation by -1 is also contained in the current data set, and clearly will have the same dynamics-based metrics. This equivalence single handedly de-facto cuts the size of the current data set by half (to 771,963 converging formulas).Third, the inherently parallelizable nature of computing metrics for each formula is leveraged.The algorithm has been adapted and deployed within the Berkeley Open Infrastructure for Network Computing [Anderson, 2004] , enabling parallel computation across thousands of volunteer computers. Assuming a typical contribution of approximately 1,000 BOINC volunteer cores, processing the expanded data set requires about one month of computational time.The approaches described enable the methodology to handle larger datasets efficiently, facilitating the analysis of more complex polynomial continued fractions within practical computational limits. Fig. 6 and Fig. 7 show alternate clustering and 2D visualization approaches (in addition to Fig. 4 ). The same set of anchors was used for all 3 versions. The tSNE algorithm (Fig. 6 ) provides a visualization that separates well between clusters, but loses out on explainability. In Fig. 7 the opposite approach was taken -using only 2 dynamical metrics for clustering allows each cluster to be defined very clearly, but information from other dynamical metrics (which can be used for better cluster separation) is lost along the way.Despite their differences, all 3 techniques point to similar results and conclusions, providing additional validation for our conclusions.The growth rate for simple continued fraction or equivalently for constant linear recurrences is well understood, and usually boils down to the matrix defining the recurrence, and its eigenvalues. In our case, the coefficient in the recurrence also depend on n, so their study is more involved, however the ideas are similar, which we now describe To find whether or not the sequence pn qn converges and if so what is its convergence rate, we note the continued fraction formula b (1)a (1) + b (2) a (2) + b (3) . . . + b(n-1) a(n-1)+0 = p n q n ,can be rewritten in matrix form asp n-1 p n q n-1 q n = n-1 1 0 b (k) 1 a (k) .In particular this implies that both p n and q n satisfy the same linear recurrence:u n+1 = a (n) u n + b (n) u n-1 , with initial conditions p 0 p 1 q 0 q 1 = 1 0 0 1 .Trying to determine if there is convergence, we use the Cauchy condition. For any m ≥ n we have thatp m q m - p n q n = m-1 n p k+1 q k+1 - p k q k = m-1 n p k+1 q k -q k+1 p k q k q k+1 = - m-1 n det p k p k+1 q k q k+1 q k q k+1 = - m-1 n (-1) k k j=1 b (j) q k q k+1 . The sequence K ∞ 1 b(n) a(n) converges if and only if ∞ 1 k j=1 b(j)q k q k+1 converges, and to the same limit L. More over, the convergence rate isϵ (n) := p n q n -L = ∞ n (-1) k k j=1 b (j) q k q k+1 .This suggests that we should understand the growth rate of both q k and k j=1 b (j). Note that the convergence and its rate might depend on the sign of(-1) k k j=1 b(j) q k q k+1 . 1. Suppose that (-1) k k j=1 b(j) q k q k+1 = 1 k d .If the signs do not alternate, then∞ n (-1) k k j=1 b(j) q k q k+1 = ∞ n 1 k d . This diverge if d = 1and has order of magnitude 1 k d-1 for d > 1. However, with alternating signs we get the smaller bound∞ 2n (-1) k k d = ∞ n 1 (2k) d - 1 (2k + 1) d = ∞ n (2k + 1) d -(2k) d (2k) d (2k + 1) d ∼ ∞ n d (2k) d-1 4k 2d ∼ 1 n d .Thus, it always converges and with better rates.2. However, for faster converging sequences we do not expect alternating sign to affect the convergence rate. For example, ifm-1 1 (-b(k)) qmqm-1= λ m for some 0 < λ < 1, then with only positive signs the limit will be λ n 1+λ while for alternating signs it will be (-λ) n 1+λ , so in any case the convergence rate is exponential.F.2 Growth rate of m-1 k=1 |b (k)|Let b (x) be a polynomial of degree d, with leading coefficient of absolute value B. Then there exists a constant C > 0 such that for any integer N we have(N e) -C ≤ N k=1 b (k) Bk d ≤ (N e) C .Proof. We may assume that the leading coefficient of b is positive. Writing b (x) = (k) = 1 + d-1 0 b j B 1 k d-j .Hence, we can find an integer constant C 0 ≥ 1 such that for all k ≥ 1 we have1 - C 0 k ≤ b (k) ≤ 1 + C 0 k .For all k large enough, all the expression above are positive, so we getln 1 - C 0 k ≤ ln b (k) ≤ ln 1 + C 0 k .With the goal of summing up these expressions from 1 to infinity, we claim that there is some constant M > 0 such that for any C ′ ∈ R, and 2 |C ′ | ≤ n < N we have thatEQUATION)Given this claim we conclude that-(C 0 ln (N ) + [M -C 0 ln (2C 0 )]) ≤ N k=2C0+1 ln b (k) Bk d ≤ C 0 ln (N ) + [M -C 0 ln (2C 0 )] .For another C large enough (independent of N ), we can start the summation from k = 1 to get-C (ln (N ) + 1) ≤ N k=1 ln b (k) ≤ C (ln (N ) + 1) .Finally, exponenting it back we get the result we wanted:(N e) -C ≤ N k=1 b (k) ≤ (N e) C .We are left to prove Equation ( 13).Using the Taylor expansion of ln (1 + x) for |x| ≤ 1 2 , we know that there is some large enough0 < M 0 such that |ln (1 + x) -x| ≤ M 0 x 2 . It follows that for 2 |C ′ | ≤ n < N we have N k=n ln 1 + C ′ k - C ′ k ≤ M 0 C ′2 N n 1 k 2 ≤ M 0 C ′2 ζ (2) .In addition, we have thatN n 1 k - N n-1 1 x dx ≤ 1, and N n-1 1 x dx = ln N n -1 . Therefore N n ln 1 + C ′ k -C ′ ln N n -1 ≤ |C ′ | + M 0 C ′2 ζ (2)is uniformly bounded.The sequence q n satisfies the linear recurrenceq n+1 = a (n) q n + b (n) q n-1 ,or in matrix form(q n , q n+1 ) = (q n-1 , q n ) M (n) 0 b (n) 1 a (n) .If both a (x) , b (x) are constant, and therefore M = M (n) is a constant matrix, then this problem reduces to simply (q n , q n+1 ) = (q 0 , q 1 ) M n . Its a standard exercise to approximate q n using the eigenvectors decomposition of M . However, in general not only M (n) is non-constant, its entries have different orders of magnitude.Thus, we would like to move to an "equivalent" system where at the very least M (n) converges to some matrix M ∞ , and then hope to show that the behavior of q n can be read from the system with M n ∞ . This equivalent system will be built in two steps: first we "balance" the matrix, so its coordinates growth rate are the same, and then taking it outside as a scalar, the remaining sequence of matrices will converge.This balancing is split into two cases according to the degrees of (n) = a(n) (n+1) d and b (n) = b(n) (n(n+1)) d, by our choice of d we see that the coefficient or the recurrence converge, and not both to zero:lim n→∞ ã (n) = A lim n→∞ b (n) = B.Here too we can also write it in a matrix form, namely(q n , qn+1 ) = (q n-1 , qn ) 0 b (n) 1 ã (n) .We now have a limit matrix, and the dynamics of such a matrix is well known. If both eigenvalues are real which are distinct in absolute value, then we expect exponential convergence. If both are non real, and therefore complex conjugate we expect it to behave like a rotation, and therefore will not converge. In both of these cases, since the eigenvalues are distinct in the limit, this holds for almost all n, so this behavior should hold in general.In the discriminant zero, the situation is much more delicate, since we can converge to zero in many ways. For example, the discriminant along the way can be negative, positive or zero. In this notes we will restrict the study only to the two real eigenvalues with different absolute values.The main goal of this section is to approximate the growth rate of a solution u n to the recurrenceu n+1 = a n u n + b n u n-1 ,where both a n , b n converge (and not both to zero) or in matrix form( vn vn+1 ) = ( vn-1 vn ) M n , M n = 0 bn 1 an , where M n → M := 0 b 1 a .The first step is the standard conjugation to a simpler matrix. Indeed, if D = P M P -1 is simpler, e.g. diagonal, then D n := P M n P -1 → D , and n 1 M i = P -1 n 1 D i P , so we more or less need to understand n 1 D i .In the constant diagonal caseD n = λ 1 0 0 λ 2 with λ 1 > |λ 2 |, we expect that for almost every initial condition (α 1 , β 1 ) D k ∼ λ k 1 .This is true as long as the initial vector is not in R • e 2 , and we have similar behaviour for other type of matrices. When the D n are not constant, we need to take a little bit more care. The image you should have in mind is the following: Instead of the two eigenvectors being on the X and Y axes, they only converge to it, so we only know that they are somewhere inside the red and blue regions. Thus, to understand this system we first need a separation condition saying that these regions are disjoint. Assuming the X-axis is the pulling axis (larger eigenvalue), we will need at least one point outside the error region around the Y axis, which we call the initial condition. Once both these conditions hold, a standard investigation of diagonalizable product will show that the point's orbit converge towards the eigenvector in the X-region. As this region shrinks to X in the limit, we see that the limit of the orbit is there as well. Suppose thatD n → D where D = λ+ 0 0 λ- with 0 ≤ λ- λ+ < 1 and let κ n = 1 |λ+| max k≥n ∥D k -D∥ ∞ .Fix some initial z 1 and letz k = (z 1 ) k-1 1 D n .Assuming that for some n we have• Separation condition: λ- λ+ + 4κ n < 1 and• Initial condition: |z n | < µn+ √ (µn-2κn)(µn+2κn) 2κn , µ n = 1 -λ- λ+ -2κ n .Then limk→∞ |z k | = 0. Note that µn+ √ (µn-2κn)(µn+2κn) 2κn∼ 1-|λ| κn → ∞ as κ n → 0, so this initial condition becomes easier to satisfy as n → ∞.Proof. First, proving the claim for 1 λ+ D n instead of D n , we may assume that the limit is D = 1 0 0 λ where λ = λ- λ+ .Next, note that wheneverµ := 1 -|λ| -2κ > 2κ > 0, we have that (µ -2κ) (µ + 2κ) = µ 2 -4κ 2 < µ. Setting ν ± (κ) = µ κ ± (µ κ -2κ) (µ κ + 2κ) 2κ ,we get that 0 < ν -(κ) < ν + (κ) are real numbers, and the condition in the assumption is|z n | < ν + (κ n ).Our main goal is to prove our process satisfies:1. |z k | < ν + (κ n ) for all k ≥ n and, 2. We have lim sup k |z k | ≤ ν -(κ n ) .Assuming these two steps are true, the full proof is not too far behind. Indeed, since D n → D, the sequence κ n := sup k≥n ∥D n -D∥ converges to zero, and note that asκ n → 0 we get that ν + (κ n ) ↗ ∞ and ν -(κ n ) ↘ 0. Assuming step (1), for k ≥ m ≥ n we have |z k | < ν + (κ n ) ≤ ν + (κ m ), and by step (2) we get that lim sup k |z k | ≤ ν -(κ m ) → 0.For the remaining of the proof, without loss of generality we may assume that n = 1 and just write κ, µ instead of κ n , µ n .To prove these two steps, consider the change fromz k to z k+1 . Writing D k = 1+ε1,1 ε1,2 ε2,1 λ+ε2,2 , since z k+1 = (z k ) D k and ∥D -D k ∥ ∞ ≤ κ, we get that |z k+1 | = ε 1,2 + z k (λ + ε 2,2 ) (1 + ε 1,1 ) + z k ε 2,1 ≤ κ + |z k | (|λ| + κ) 1 -κ -κ |z k | .Note that the final denominator is positive, so that the last inequality is valid. Indeed, using the conditions of the claim we get1 -κ (1 + |z k |) ≥ 1 -κ 1 + µ + (µ -2κ) (µ + 2κ) 2κ > 1 -(κ + µ) = |λ| + κ > 0.Thus, we can rewrite the inequality asEQUATIONThe goal now is to show that if |z k | is "large", then |z k+1 | is much smaller, and if |z k | is small, then |z k+1 | cannot increase too much.A simple computations shows that the eigenvalues of this matrix areγ ± = |λ| + 1 ± (µ + 2κ) (µ -2κ) 2 ,and since(µ + 2κ) (µ -2κ) ≤ µ ≤ 1 -|λ|, we get that γ + > γ -> 0.Finally, the corresponding (right) eigenvectors arev ± = ν ∓ 1 .To simplify the notations, let us conjugate by the matrixT = ν + ν - 1 1 to obtain T -1 M ε T = γ -0 0 γ + .Note that the Mobius mapT -1 (z) := 1 ν + -ν - 1 -ν - -1 ν + (z) = - z -ν - z -ν + = -1 + ν --ν + z -ν + sends ν -→ 0, ν + → ∞ and 0 → -ν- ν+ < 0.In particular, it is monotone increasing on [0, ν + ), so that our two steps from above are equivalent to1. T -1 (|z k |) ∈ [-ν- ν+ , ∞), 2. lim sup k T -1 (|z k |) ∈ -ν- ν+ , 0 ,and the claim's original assumption is thatT -1 (|z 1 |) ∈ [-ν- ν+ , ∞).However, now this claim is simple, since in these notations we get thatT -1 (M ε (|z k |)) = T -1 M ε T T -1 (|z k |) = γ - γ + • T -1 (|z k |) , and 0 < γ- γ+ < 1. Thus, if T -1 (|z k |) ∈ [-ν- ν+ , ∞), then so is T -1 (M ε (|z k |)) ∈ [-ν- ν+ , ∞), so by Equation ( 14) and the monotonicity of T , we obtain thatT -1 (|z k+1 |) ≤ γ - γ + • T -1 (|z k |) ,which implies the two steps.Returning back to the recursion, we get the following.Suppose that we have a solution to the recurrencev n+1 = a n v n + b n v n-1 , where a n → a, b n → b and suppose that λ ± are the roots of x 2 = ax + b with 0 ≤ |λ -| < λ + . Writing κ ′ n = 1 |λ+| max k≥n max {|a k -a| , |b k -b|} and C (λ ± ) := 1+|λ+| |λ+-λ-|, Assume that for some n we have• Separation condition: λ- λ+ + 4C (λ ± ) κ ′ n < 1 and • Initial condition: λ --vn vn-1 ≥ C (λ ± ) κ ′ n |λ+-λ-| 1- λ - λ + -4C(λ±)κ ′ n , Then v n v n-1 → λ + .Proof. Set M n = 0 bn 1 an and M = ( 0 b 1 a ) as in the beginning of this section. With P =1 λ+ 1 λ- and P -1 = 1 λ--λ+ λ--λ+ -1 1we have that D = P M P -1 = λ+ 0 0 λ-. We would like to apply Lemma F.3.2 to the matrices D n = P M n P -1 .For the separation condition on the infinity norm, we haveP M n P -1 -D ∞ = P (M n -M ) P -1 ∞ = 1 |λ --λ + | 1 λ+ 1 λ- 0 bn-b 0 an-a λ--λ+ -1 1 ∞ = 1 |λ --λ + | 1 λ+ 1 λ- b-bn bn-b a-an an-a ∞ ≤ C(λ±) 1 + |λ + | |λ + -λ -| ∥M n -M ∥ ∞ .Thus, the separation condition of this theorem implies the separation condition of Lemma F.3.2:λ - λ + + 4κ n ≤ λ - λ + + 4C (λ ± ) 1 |λ + | max k≥n ∥M n -M ∥ ∞ < 1Next, for the initial condition , setting(v k-1 v k ) := (v 0 v 1 ) k-1 1 M n = (v 0 v 1 ) P -1 k-1 1 D n Pwe have(α k , β k ) = (v 0 v 1 ) P -1 k-1 1 D n = (v k-1 , v k ) P -1 .Setting z n = βn αn , we get that|z n | = β n α n = -λ + v n-1 + v n λ -v n-1 -v n = 1 + λ + -λ - λ --vn vn-1 ≤ 1 + λ + -λ - λ --vn vn-1 = ( * ) .Using the assumption that λ--vn vn-1 ≥ C (λ ± ) κ ′ n |λ+-λ-| 1- λ - λ + -4C(λ±)κ ′ n ≥ κ n |λ+-λ-| 1- λ - λ + -4κn, we see that the expression above is( * ) ≤ 1+ |λ + -λ -| κ n |λ+-λ-| 1- λ - λ + -4κn = 1+ 1 -λ- λ+ -4κ n κ n = 2µ n -2κ n 2κ n < µ n + (µ n -2κ n ) (µ n + 2κ n ) 2κ n .This was the second condition needed for Lemma F.3.2 , so we can now conclude that1 + λ + -λ - λ --vn vn-1 = β n α n → 0which implies that vn vn-1 → λ + .We return now to the original problem withα = K ∞ 1 b(n) a(n)and assume that a (x) , b (x) have degrees d a , d b . As mentioned before, we split our study into two cases:The balanced case Assume that d b = 2d a = 2d, and let A, B be the leading coefficients of a (x) , b (x) respectively.In this case the limit matrix is M ∞ = 0 B 1 A , and we assume that the roots λ ± of x 2 = Ax + B satisfy 0 < |λ -| < λ + . Using Theorem F.3.2 once the two conditions hold, we obtainq n+1 q n (n + 1) d = q n+1 / (n + 1)! d q n /n! d → λ + , implying that q n = n! d λ n + exp (o (n)).As for the product of the b (k), using Claim F.2 we have thatN k=1 |b (k)| = exp (o (N )) • B N • N ! 2d .Putting them together as in the error rate expression, we get :m-1 k=1 |b (k)| |q m-1 q m | = |B| m-1 • (m -1)! 2d (m -1)! d (m)! d λ 2m-1 + exp (o (m)) = ( * ) . Note that |B| = |det (M ∞ )| = |λ -λ + |, so that the expression above is ( * ) = |λ -/λ + | m • exp (o (m)) = exp (m log |λ -/λ + | + o (m)) .Thus, for given ε > 0 where λ- λ+ + ε < 1, and for any m large enough we see that ( * )≤ λ- λ+ + ε m-1. We conclude that the error rate for all n large enough is bounded from above byp n q n -α ≤ ∞ m=n+1 m-1 k=1 |b (k)| |q m-1 q m | ≤ ∞ m=n+1 λ - λ + + ε m-1 = λ - λ + + ε n 1 1 -λ- λ+ + ε . It follows that ln p n q n -α ≤ n ln λ - λ + + ε -ln 1 - λ - λ + + ε ∼ n ln λ - λ + .The unbalanced caseEQUATION(n (n + 1)) da = 0. This time the two roots of x 2 = Ax + 0 are λ = 0, A. If needed, we can use a simple continued fraction inflation K ∞ 1 (-1) 2 b(n) (-1)a(n) and assume that A > 0. Using Theorem F.3.2, if the two conditions hold, we obtainq n = n! d A n exp (o (n)).Letting B be the leading coefficient of b (x) in absolute value, Claim F.2 implies thatN k=1 |b (k)| = exp (o (N )) • BN • N ! d b .Again, together we obtain thatm-1 k=1 |b (k)| |q m-1 q m | = Bm-1 • (m -1)! d b (m -1)! d m! d A 2m-1 exp (o (m)) = 1 (m -1)! 2da-d b • B A 2 m exp (o (m))Similarly to the previous case, given ε > 0, and using the fact that 2d a -d b ≥ 1, for all n large enough we obtainp n q n -α ≤ ∞ m=n+1 m-1 k=1 |b (k)| |q m-1 q m | ≤ ∞ m=n+1 1 (m -1)! 2da-d b • B A 2 + ε m-1 = 1 n! 2da-d b B A 2 + ε n ∞ m=0 n! (n + m)! 2da-d b • B A 2 + ε m ≤ 1 n! 2da-d b B A 2 + ε n ∞ m=0 1 m! 2da-d b • B A 2 + ε m .The infinite sum in the last exression converges to some finite limit C, so we conclude thatln p n q n -α ≤ (d b -2d a ) ln (n!) + n ln B A 2 + ε + ln C ∼ (d b -2d a ) n • ln |n| .-4n 2 + 5n + 4 n + 1 -n -4 1 1 - 5 335 3 -160 log (2) -4 4n 2 + 4n 2n -2n -2 1 - 4 -2+4 log (2) -3 n 2 n -n -2 1 - 3 3-3 log (2) -3 n 2 + 2n n -n -2 1 - 3 -15 2 +12 log (2) -3 n 2 + 2n + 1 n + 1 -n -1 n + 3 2 1 3 - 10 3(-20+30 log (2)) -3 n 2 + 4n + 3 n + 1 -n -3 1 1 - 4 34-48 log (2) -3 n 2 + 4n + 4 n + 2 n -2 n + 5 2 6 5 - 21 5( 147 2 -105 log (2)) -2 n 2 + n n -n -1 1 - 2 -2+4 log (2) -2 n 2 + 3n + 2 n + 1 -n -2 1 1 - 3 9-12 log (2) -2 4n 2 2n -2n 1 -1 log (2) -1 n 2 n -n 1 -1 log (2) -1 n 2 + 2n + 1 n + 1 -n -1 1 1 - 2 2-2 log (2) -1 n 2 + 4n + 4 -n -2 n + 2 1 2 - 3 -3 2 +3 log (2) -4 4n 2 -1 2n -1 -2n -1 1 -3 -3 2 + 3π 4 -1 -2 4n 2 -4n + 1 2n -1 -2n + 1 1 -4 π -1 -2 4n 2 + 4n + 1 2n + 1 -2n -1 1 1 -3 3-3π 4 -5 4n 2 + 2n -2 2n + 1 -2n -2 1 - 4 -20 3 + 16 √ 2 3 -1 -5 4n 2 + 2n 2n + 1 -2n n + 3 4 - 32Γ( 11 4 ) 9 - -56Γ( 5 4 )Γ( 7 4 ) π -14 Γ( 3 4 ) -1 3 -3 4n 2 -2n 2n -1 -2n 1 -1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? Answer: [Yes] Justification: Yes, the main results, motivations and aspirational goals are included in the abstract (section 1) and introduction (section 2). Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? 
discovered formulas for mathematical constants: 
 The first step in validating the dynamical metrics approach is using basic heuristics on the metric space to find PCFs related to mathematical constants. There are some PCFs in the dataset that have a known irrational limit (like the examples in Eq.1 and the PCF familyB A + B A + . . . = 2B A + √ A 2 + 4Bfor constant A and B), so we expected to find some of them. Through this test, we also found previously unknown PCF formulas related to mathematical constants. Note that known mathematical formulas are both the anchors for labeling and a test set in our method.A natural heuristic is inspecting PCFs with δ ≈ 1, marking their limits as irrational. Another heuristic we used is focusing on PCFs with η ′ ≈ 0, as it was a very strong indicator for mathematical constant (a) Projecting the data on the δ vs. η ′ (q n factorial coefficient) plane, it is easy to see the emerging subsets. We focus on PCFs with η ′ ≈ 0, as a previous work [Elimelech et al., 2023] indicated this as an important property. (b) Clustering in the δ vs. γ ′ (q n exponential coefficient) plane shows examples of common properties within a cluster, like rationality or convergence to a specific constant (up to a linear fractional transformation). Focusing further on the deg(B) > 2deg(A) cluster (as it is a clear anomaly in the η ′ ≈ 0 subset), we used a PSLQ algorithm to identify links between these formulas and mathematical constants. This identification was feasible since a preliminary step identified a promising subset ∼ 5, 000 times smaller than the initial dataset. (c) The result of this clustering and identification procedure is a structured arrangement of formulas that reveal a range of novel formulas related to constants such as π, ln(2), √ 2, Gauss' constant, and Lemniscate's constant. (d) Keeping only PCFs with B 2 = 1 we are left with a highly symmetrical "checkerboard pattern" of formulas for π and ln(2), which was generalized into infinite formula families hypotheses (see section 4.3). Error bars not shown for visual clarity, see Appendix A for a discussion regarding measurement errors.formulas in a previous work [Elimelech et al., 2023] . Combining the two gives a subset (see Fig. 3a top left) that contains PCFs such as:EQUATIONRemoving the limitation on qn growth rate, one can find the cot(1) formula shown in Fig. 2a :EQUATIONOn the other hand, relaxing the limitation on δ, focusing only on η ′ ≈ 0, a rich structure emerges (Fig. 3b ). Diving deeper into the B-dominated subset, we find formulas (Fig. 3c ) for the Gauss constant G GA [Finch, 2003] :EQUATIONLemniscate constant L Lemniscate [Finch, 2003] :4 + 2 . . . + 4n 2 -2n 4 + . . . = -6 L Lemniscate -4 (9)As well as for second order roots, π and ln(2) (see section 4.3). Note that unlike the formulas in Eq.6 and Eq.7, which are analytically proven, the formulas in Eq.8 and Eq.9 are (to the best of the authors' knowledge) novel. Their limits were numerically validated to a large precision, yet formal proofs for these formula hypotheses remain an open challenge.It should be noted that usually in number theory, a bigger δ is considered "good", whereas a smaller (often negative) δ is considered "bad". We use δ as a metric, without "judgment". These novel formulas (Eq.8, Eq.9 and the infinite family of formulas shown in section 4.3), which have the "bad" δ ≈ -1, are a demonstration that our "non-judgmental" approach is successful.This section shows that clusters in the latent space of dynamics-based metrics are successful in grouping together different formulas in a way that exposes their shared properties, such as the mathematical constant to which they relate.Looking at the top left cluster in Fig. 3b (defined by qn exponential coefficient < 0.6 and δ > 0.9), we recognize the canonical form of the Golden Ratio PCF (shown in Eq.1). This cluster also contains 21 additional PCFs, with different generating polynomials, some of higher degree. As it turns out, all of them are linear fractional transformations of √ 5 (see Appendix C), which were labeled automatically by the formula discovery algorithm (Fig. 1 ). Another example of property conservation within clusters is the rational cluster marked in green on Fig. 3b . The limits of the PCFs in this subset are varied, and its spread is real (i.e., not only due to numerical imperfections). Yet, all the PCFs in this cluster converge to rationals -which is not directly measured by any of the latent space dimensions. Note the multi-anchor clusters of e and e 2 , as well as the algebraic roots: these clusters failed to single out a specific constant, yet relate to constants of similar nature -suggesting meaningful clustering nevertheless. For the sake of visualization the algorithm stopped after the second iteration. In a standard run these multi-anchor clusters would have been separated via additional metrics.As mentioned in section 4.1, the deg(B) > 2deg(A), η ′ ≈ 0 cluster, contains many formulas of interest (see Fig. 3c and d ). They were discovered via a PSLQ algorithm, identifying linear fractional relations between the limit values of PCFs in the subset and notable mathematical constants (such as π or e). This is a computationally heavy operation, and it would be challenging to run it on all 1.5M formulas in the data set. Yet by first identifying the promising clusters, we reduce the search space ∼ 5, 000 times, allowing for a deeper inspection of each PCF.Once the "checkerboard" pattern in Fig. 3d was discovered, the hypothesis was expanded into 2 infinite families of PCFs with sub-exponential convergence relating to π and ln(2):• a n = i + 2j + 1, b n = n 2 + (i + k)n,with integers i, j ≥ 0, and k ∈ {0, 1}. This is expected to be related to π if k = 1, and to ln(2) if k = 0 (in fact, this pattern can be generalized even further, into a novel 3-dimensional Conservative Matrix Field, provided in Appendix C. See [Elimelech et al., 2023] for the definition of Conservative Matrix Fields).Another formula family was discovered via clustering in the γ vs. γ ′ space. The algebraic roots subset (marked by a green circle in Fig. 2c ) was generalized into:• a n = -2n + j -1, b n = -n 2 + jn + k for integer j, k such that b n has real roots that are not positive integers. This is expected to converge to a root of b n .These are novel experimental results and mathematical hypotheses -awaiting proof. 3b . Formulas that are equivalent to an existing element in the original dataset are marked with "o", while "x" marks truly novel formulas. The higher-degree PCFs fit, using the same dynamical metrics, into clusters trained on lower-degree PCFs, showing the ability to identify novel formulas despite being of mathematical forms not seen during training.1 6 𝜋𝜋 2 = 1 + -1 5+ -16 13+ -𝑛𝑛 4 2𝑛𝑛 2 +𝑛𝑛+1+⋱The final validation for the dynamical metrics approach is to show its potential even on test sets of a different mathematical structure than the training set. For that purpose, the clusters created based on the original dataset were treated as a classifier, and a test set of higher-degree PCFs (up to 3 rd degree A n and 4 th degree B n , coefficient range [-2,2]) was created, measured and classified.Naturally, not all high-degree PCFs fit neatly into the existing clusters (as they represent constants that were not present in the original dataset), but some were correctly identified and labeled, discovering novel formulas in the process (see Fig. 5 for a sample of the results). 
discussion and outlook: 
 This work marks an important step toward the vision of automated on-demand formula creation in mathematics. Going beyond all previous algorithms in this field, we connect the challenge of formula creation to robust ML methods. This methodology provides a wide variety of automatically generated formulas, including both previously known and previously unknown ones, exposing their underlying mathematical structure and enabling new proofs.The next research step directly building on our methodology could help to finally reveal the complete intricate mathematical structure of PCFs. For example, starting with the "band-structure" found in Fig. 2c , or with clusters of formulas with various structures representing the same mathematical constant. Further exploration of our conjectures from section 4 could have more impact on mathematics, perhaps achieving further generalizations and prescriptive formula generation.The technique presented here can be applied to a larger scope of continued fractions and for completely different types of formulas. For more general continued fractions, dynamical metrics such as the numerical trajectories and the corresponding sequences of δ (in addition to its asymptotic value) hold valuable information even in continued fractions that do not converge at all. We expect these dynamical metrics to provide a "fingerprint" for wider families of formulas and perhaps even for the mathematical constants themselves. This approach was directly applied in this work to higher polynomial degrees, larger polynomial coefficients, and can be expanded to continued fractions not based on polynomials. Looking beyond continued fractions, metrics that are derived from the dynamics of a numerical calculation of certain formulas are an especially good fit for automated computer-assisted investigations. Such metrics can be measured for a variety of mathematical structures, including ones whose evaluation is iterative or recursive, that are defined via an infinite sum, or any other process which produces rational approximants. Any such mathematical structure can be measured, clustered, and identified using the proposed method -treating the generating functions as a black box. We believe that such dynamical metrics can unveil patterns and underlying structures in broad fields of mathematics and in other areas of science.To exemplify this universal concept, we looked into higher depth recursion relations, which are a promising research direction because little is known about their global structure, yet they are involved in several important conjectures. For example, the best rational approximation formula known for Euler's gamma constant is constructed via such a recursion relation [Aptekarev, 2009] . This family of formulas is broader than continued fractions, yet they can be described by the same metrics as PCFs. Another type of mathematical structure successfully analyzed using the same method is hypergeometric functions, showing the applicability of our measurement-clustering-generation approach to a bigger family of mathematical functions. This generalization can be useful in a wide variety of contexts, such as investigations of integral formulas (e.g., Beukers-type integrals [Beukers, 1979 , Dougherty-Bliss et al., 2022 , Brown and Zudilin, 2022] ).Our work was based on a limited-size dataset and on a small set of metrics. It would be intriguing to test the extracted conjectures on larger datasets, which can help reveal additional, more intricate, phenomena. Considering the success we had using a relatively small set of metrics, we would like to use an order-of-magnitude larger set of metrics and find what new predictions can be recovered. In fact, the creation and evaluation of the metrics themselves can be automated.Taking a broader perspective, the methodology presented in this work can be seen as a general prescription for tackling scientific discovery challenges in mathematics and theoretical physics that rely on numerical evaluations and generalizations. Such an advance is especially exciting for such challenges that were considered in the past to require intuitive leaps of creativity.When characterizing PCFs, we use several metrics extracted from the dynamic behavior of the formula:• The growth coefficients η, γ, β (of the form n! η • e γn • n β ) of the convergence rate ϵ(n).• qn (as defined in Eq.3) growth coefficients: η ′ , γ ′ (of the form n! η ′ • e γ ′ n ).• The δ (as defined in Eq.3) calculated using the Blind-δ Algorithm described in section 3.4.To measure the growth coefficients of qn and ϵ(n), the values of log (ϵ(n)) (see section 3.4) and of log (q n ) were evaluated up to depth 1000.The most resource-intensive values that are generated are p n , q n and gcd(p n , q n ) -all other values are calculated from them (and require less precision). For the worst case PCF this requires 36MB of memory (without optimizations) and ∼ 1.9 seconds of run time on a single core of a basic workstation, which translates to an upper cap of ∼ 900 hours for the whole data set. In practice we used a high power cluster with 64 cores, which ran each iteration of the measurements in ∼ 8.5 hours.Once these values are calculated, using scipy [Virtanen et al., 2020] and numpy [Harris et al., 2020] a fit of the form log n! η • e γn • n β was calculated for qn and ϵ(n), producing the dynamical metrics.A curve fit using 1000 points is a fairly heavy operation, unsuited for large scale investigations. Instead, we used an extreme down-sampling. Specifically, only 5 points were used for the fit. One may justifiably wonder if 5 data points are sufficient to fit accurately enough the desired metrics.A test comparing between a 5 data point fit and a 1000 data point fit was done. As the test set, 50 PCFs were randomly chosen out of each of 9 categories (450 total test cases). The categories were all combinations of deg(a) = 0, 1, 2 and deg(b) = 0, 1, 2. Focusing on the dominant coefficients (γ and η), for each case, a full (1000 point) fit was performed (producing γ f , η f ), and compared to the down sampled fit of 5 points (producing γ p , η p ). We tested 2 methods of choosing the 5 points, even (i = 6, 206, 406, 606, 806) and logarithmic (i = 6, 125, 250, 500, 1000). The relative error was then calculated (|γp-γ f | |γ f | and |ηp-η f | |η f |) for ϵ(n) and qn . The relative errors were then averaged over the test set (results summarised in table 2) -showing the 5-point fit to be almost as good as the full 1000-point fit. In our measurements we use the logarithmic point distribution as it gives better results for most metrics. After a PCF has been identified as potentially connected to a known constant C, it undergoes a verification process in order to establish whether it indeed converges to a known mathematical expression involving C. This process consists of several step:Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: The git repository will is linked in this camera-ready version. The data set is mathematical formulas and constants -so freely available to all.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: In this work we are using unsupervised recurring clustering as the main machine learning method, and the results are mostly validated using the (test) subset of known PCF formulas for mathematical constants. This subset was included in the full data set but was not given any special treatment during clustering (as mentioned in section 4.1).Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: The main errors in this work originate in the dynamical metrics measurement stage. Error bars are not shown (for visual clarity), but these errors are discussed in Appendix A, and mentioned in every figure that shows measured dynamical metrics.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: There are no special or extreme resource requirements, so we did not focus on this question, but a brief discussion of required memory and runtime is in Appendix A.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that did not make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: This work does not involve human participants, special data sets, does not have implications for the broader public (only for scientists) and cannot be used to do harm.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: Although our work is applied ML, its impact is on mathematical research, not society at large.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The methods in this paper do not pose such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: The only existing assets that were used are open source Python packages, that were properly credited in Appendix A.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.-2+2 √2 -1 
 
 Appendix:  


Decision: This paper presents an original contribution to automatized discovery of mathematical formulas. While the machine learning methods it uses are rather basic, it does an extensive job of evaluating a large number of continued fractions. It discovers previously unknown (at least to the authors and the reviewers) formulas. This is a valuable contribution to the field and will be of interest to the NeurIPS community. 

The reviews contain many points that improve the presentation of the paper. The authors are strongly encouraged to carefully take them into account in order to increase the impact of their interesting work.

Reviews: review: paper_summary: paper_summary: The paper addresses a long-standing challenge in number theory by proposing a new methodology for the categorization, characterization, and pattern identification of mathematical formulas, specifically Polynomial Continued Fraction (PCF) formulas. The authors introduce metrics based on the convergence dynamics of these formulas, enabling the first automated clustering of mathematical formulas. The methodology is demonstrated on a dataset of 1.7M PCF formulas, leading to the identification of both known and previously unknown formulas for significant mathematical constants such as π, ln(2), Gauss, and Lemniscate constants. The uncovered patterns allow for the generalization of individual formulas to infinite families, revealing rich mathematical structures. This work sets the stage for a generative model capable of creating continued fractions with specified mathematical properties, potentially accelerating the discovery of useful formulas.
strength_weakness: strength_weakness: The paper has some weaknesses that should be addressed to improve its quality:


While the addressed task and the methodology is unique, it may not be the best fit for a machine learning conference like NeurIPS. The focus on mathematical discovery might be better suited for a specialized conference or journal in mathematics or computational mathematics.

The paper’s writing can be further improved for clarity and conciseness. The abstract could be more concise, and the figures are currently blurry and not well-organized, detracting from the overall presentation.

The study is based on a limited-size dataset and a small set of metrics. Expanding the dataset and incorporating a broader range of metrics could enhance the robustness and applicability of the findings.

The newly identified formulas for significant constants need to be further verified. Ensuring their correctness and utility is crucial for the validity of the contributions.
questions: questions: How do the authors plan to further verify the newly discovered formulas for significant constants such as π and ln(2)?

Can the authors provide more details on how their methodology can be scaled to larger datasets and more complex formulas?

How can the approach be adapted or extended to other types of mathematical formulas beyond Polynomial Continued Fractions?
limitations: limitations: The paper could benefit from a more detailed discussion on the impact of the limited dataset size and the potential benefits of using larger and more diverse datasets.

The choice of metrics is crucial for the clustering and characterization of formulas. Providing a rationale for the selected metrics and discussing potential additional metrics would strengthen the paper.

Further elaboration on the methods for verifying the newly discovered formulas would enhance the credibility of the findings.
paper_summary: paper_summary: They generate continued fraction formulas and test if they evaluate to mathematical constants. 
They introduce a distance metric to compare formulas.
They discover novel formulas for known constants.
strength_weakness: strength_weakness: It generates formula hypotheses, so we do not always know if the formulas are actually correct

Rather limited structure of the formulas
questions: questions: all formulas are continued fraction formula of quadratic polynomials? is that not rather limited?

in (1) you give the example tan(x), but in the definition of a/b (p4 118), there is only n and no x.  So you cannot actually get a formula for tan(x), can you? 

is there a list of all 1,543,926 formulas? how many are equal to pi or e? is there pi^pi or e^e among them? or ln(pi)?
do you know how many are irrational? or even how many are transcendent?

how many are proven to be correct and how many remain hypotheses?

have you consulted mathematicians if they are going to use these formulas for anything?

>(4)

is it necessary to introduce q̃_n and not use q_n without tilde there?
limitations: limitations: yes
paper_summary: paper_summary: The paper presents a classification of 1.5 million polynomial continued fractions (PCF), continued fractions having as coefficients the integer values of two polynomials, $A(n)$ and $B(n)$, with $A$ $B$ of degree two, with integer coefficients in $[-5,5]$. PCF are classified according to the asymptotic properties of the sequence of approximation errors (difference between the convergents and the limit), the asymptotic properties of the denominator of the convergents (in simplest terms), and the measure of irrationality, which is the difference between the logs of the two previous metrics.

Along these metrics, the authors observe that groups of PCF with the same limit tend to cluster together, and that by "anchoring" these clusters on known formulas, one can discover new PCF decompositions of mathematical constants.
strength_weakness: strength_weakness: The link with AI or machine learning is not completely clear to me. Whereas the clustering techniques used in the paper, and the t-SNE projections used to represent the results, are commonly used by ML practitionners, most of the analyses conducted in this paper amount to descriptive statistical analysis of a specific mathematical dataset. To demonstrate a possible link with Machine Learning, it would be useful that the authors discuss (and perhaps demonstrate) how their approach scales to large sets of PCF, by letting the coefficients and degrees of $A$ and $B$ grow larger. 

I lean towards acceptance because of the potential interest of such approaches in AI for Science.
questions: questions: * Figure 1 caption is very long, in the interest of clarity, it might be worth describing your methodology in a specific section (3.2?)
* l.96: couldn't we assume that the convergents $p_n/q_n$ are always in simplest terms? this would avoid having to introduce of $\tilde q_n$. Besides, I believe it is assumed in the usual definition of the measure of irrationality.
* in section 3.1, you explain that the irrationality measure is either 0 or larger than 1. You then claim that the blind-$\delta$ method provides a good estimator of $\delta$, and figure 2.a seems to support this, yet in figure 2.b most estimates of $\delta$ are below $1$, and some are even negative. What happens? Doesn't this compromise the use of blind-$\delta$ as an estimator of the irrationality measure?
* Section 3.5, can you elaborate on "representation power", what it is? why is it important?
* Table 1: The exponential factor of the growth coefficient seems useless, why is it? Also, could the larger value of the Davies Boulding Index the irrationality measure be a sign of the problem with its estimation?
* Figure 4: can you provide a description of the axes in the tSNE graph? (maybe switch to PCA for more explainability)
* l.29: shouldn't Lambert's original paper be quoted, instead of a modern compilation?
* l. 100 "equal" (typo)
* l.135: the error rate $\epsilon$ is used before it is defined (in section 3.4)
* Figure 4, constant $C1$ is not defined anywhere, lemniscate?
limitations: limitations: yes
paper_summary: paper_summary: The authors propose an algorithm to filter, cluster and label polynomial continued fractions. It is based on several features linked to the asymptotic behaviour of their rational approximations. The authors detail how these are computed. They apply this algorithm to a large set of formulae ; they discuss the properties of the obtained clusters and highlight some novel PCFs linked to known mathematical constants.
strength_weakness: strength_weakness: The claim that "we connect the challenge of formula creation to modern approaches in AI for Science" seems a bit bold. To my understanding this work does not involve modern AI in the sense that it requires manual and careful feature extraction and that nothing is learnt.
questions: questions: - Part 3.4 was a bit confusing. The Blind-δ Algorithm is presented as a proxy for delta eq. 4. Is it not rather to estimate the approximation error epsilon ?

- There is an inconsistency : part 3.2 it is n^beta vs table 1 P(n) and no beta ; is it because beta is hard to estimate ? is it usefull for clustering ?

- Could the authors comment on the fact that the predicted delta eq. 5 can be quite far from the actual delta (fig. 2a, fig. 2b) ? Yet it seems it provides relevant information for clustering ; what does it capture ?

- More generally, did the authors try to extract more features, other metrics, not necessary having a mathematical interpretation, in an automated way ?

l. 100 "equval"

l. 152 "but how is it related to the actual series delta?" is not clear

l. 197 "We’ll start" informal
limitations: limitations: -

Comments: comment: I thank the reviewers for their replies, the precisions and the additional results they give.
Thank you very much for your replies, which clarify a number of my questions. I will keep my rating, and believe this paper is a good fit for NeurIPS.