Need a generalization of RootApproximant to recognize linear combinations over algebraic numbers

Question

RootApproximant does a very good job when I need to recognize an algebraic number and when enough of its digits are known (or even when an unlimited number of digits can be obtained from a numerical computation given enough time).

But often I need to recognize linear combinations of the form $\alpha+\beta\cdot\tau$ where $\alpha$ and $\beta$ are unknown algebraic numbers, and $\tau$ is a known fixed transcendental number (e.g. Pi, E or Log[2]).

Some simple cases can be solved using WolframAlpha lookup:

WolframAlpha["2.421441469079183123",
    IncludePods -> "PossibleClosedForm",
    TimeConstraint -> ∞]

Is there a way to solve this problem in general in Mathematica?

Answer 1

I can offer some ideas but it would take work to put together into a package. I'll illustrate with the example sqrt(3)+sqrt(5)pi.

val = Sqrt[3] + Sqrt[5]*Pi;

Set up an vector of values involving powers of both pi and val. We will work at fairly high precision (300 digits). I chose to go to degree 6 in each, which is overkill for this example but of course one would not know that if one did not know the number under consideration. (Likewise one does not know what precision might be needed. That sort of thing involves trial and error unless a priori bounds are known for both coefficients and degree.)

e = 300;
vals = Table[{val, Pi}^j, {j, 0, 6}];
tvals = Transpose[vals];
prods = Union[Flatten[Outer[Times, tvals[[1]], tvals[[2]]]]];
n = Length[prods];

Now set up an integer lattice with this vector as column 1 and an identity matrix tacked on to the right. When we reduce it, the matrtix on the right in effect records the linear operations used in that reduction. That is to say, if a linear relation appears between the powers we set up, it will be recorded in the corresponding row elements to the right of the first column. (If this is not making sense, recall how one is taught to invert matrices. This is the same idea.)

LatticeReduce wants exact values for input so we multiply by 10^precision and round. This is also useful because we want the elements to the right of column 1 to be small relative to column 1, so that their magnitudes do not affect the result of the reduction.

iden = IdentityMatrix[n];
lat = Transpose[Prepend[iden, Round[10^e*prods]]];
redlat = LatticeReduce[lat];

Now things proceed by hand. We check the first column in the reduced lattice to see if any elements are small.

redlat[[All, 1]]

Out[280]= {0, 3, 19, 20, 11, 3, 10, 0, -3, 8418239, 1665689, \
12765195, -9350569, 8425940, 1294210, 7162793, 1256878, 4894036, \
-357934, 421674, 9371963, -16456134, 2330713, 14461320, 7543239, \
14075640, 12607973, -3613106, 6864536, 2542139, -15683973, 14355643, \
-18398453, 1357865, -3346617, -5348079, 17039740, -7005888, 20053768, \
-6467532, -6804288, -12349164, 2944867, -11953635, -3433782, 3208561, \
-4995208, -10462168, 10265253}

Jackpot. Not only are several small, one is indeed zero. But we need to determine whether the true value is likely to be zero, since this could be the result of roundoff error.

I will verify that, for this example, we really have a zero relation.

Most[redlat[[1]]].prods // Expand

Out[282]= 0

In an actual problem we only know those values in prods to finite precision so we would just use Dot and see whether the result was near zero. For example, in this instance row 8 column 1 is also zero. But this is a "spurious" relation resulting, I would suspect, from roundoff error.

In[284]:= Most[redlat[[8]]].prods // N

Out[284]= -5.773318038*10^6

While this is not encouraging, something I fond was that many nonzero small values in column 1 corresponded to "near relations" (in-laws?) that were off by exact integers. Once the exact integer is found by numeric evaluation, one can construct a polynomial pretty much as I do below, but subtracting that integer to obtain zero.

Anyway, back to the found relation.

Most[redlat[[1]]].prods

9 \[Pi] - 30 \[Pi]^3 + 25 \[Pi]^5 - 
 6 \[Pi] (Sqrt[3] + Sqrt[5] \[Pi])^2 - 
 10 \[Pi]^3 (Sqrt[3] + Sqrt[5] \[Pi])^2 + \[Pi] (Sqrt[3] + 
    Sqrt[5] \[Pi])^4

This next part could be automated but I'll just go forward by hand. One can write this as a polynomial in y, over Q(pi). (This could be done in automated fashion from knowing what power products were in each position of the prods vector. I use the fact that I know, symbolically, what val is and how to change it to y.)

The polynomial is

In[316]:= pol = Most[redlat[[1]]].prods /. val -> y

Out[316]= 
9 \[Pi] - 30 \[Pi]^3 + 25 \[Pi]^5 - 6 \[Pi] y^2 - 
 10 \[Pi]^3 y^2 + \[Pi] y^4

So we are interested in roots of y as functions of Pi. We'll set this to zero and solve for y. Then we take the solution that numerically matches the value of interest.

solns = y /. Solve[pol == 0, y];
soln = Select[solns, Abs[# - val] < 10^(-20) &]

During evaluation of In[330]:= N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -Sqrt[3]-Sqrt[5] \[Pi]+Sqrt[3+2 Sqrt[15] \[Pi]+5 \[Pi]^2]. >>

Out[331]= {Sqrt[3 + 2 Sqrt[15] \[Pi] + 5 \[Pi]^2]}

Unfortunately this is not readily simplified. There are some ideas here for how one might take that next step. I'll just show the computation.

fax = 
 Select[FactorList[pol, 
    Extension -> {Sqrt[2], Sqrt[3], Sqrt[5]}][[All, 1]], ! 
    FreeQ[#, y] &]

Out[338]= {Sqrt[3] - Sqrt[5] \[Pi] + y, Sqrt[3] + Sqrt[5] \[Pi] - y, 
 Sqrt[3] + Sqrt[5] \[Pi] + y, Sqrt[3] - Sqrt[5] \[Pi] - y}

One can easily read off the correct root from this list.