Graeffe's root squaring method

[Albert Chan]

Theoretically, we could use Graeffe's root squaring method to solve [VA] SRC #012b ROOT
Practically, we get hit with overflow, way before 10-12 digits accuracy mandatory requirement.

The Graeffe Root-Squaring Method for Computing the Zeros of a Polynomial

Cleve Moler, July 11, 2016 Wrote:As far as I know, no one has ever made a serious piece of mathematical software out of the Graeffe Root-squaring method. If you know of one, please contribute a comment.

Still, I think root squaring is interesting topic by itself.
I like to share the PM's I sent (total 9).

[Albert Chan]

11-18-2022, 12:42 PM

Werner Wrote:- the Graeffe squaring process quickly overflows; in this case, you obtained 12+ digits of the correct answer - but you cannot know that unless you performed one or more additional steps, which would overflow. I tried it out with smaller degrees (in excel..) and I could only confirm about 6 digits before overflow. And then remains the Newton-Rhapson step to refine the value, but all I have is the abs value, not the estimate of the root, so I have n possible complex roots to verfiy. The link to Cleve Moler's article says just that. (thanks for that link BTW!)

P(x) = 2 + 3x + 5x^2 + 7x^3 + 11x^4 + ...

All P roots inside complex unit circle.

We can assume min root real part closer to -1 than +1, to keep imag part small.
Root squaring process gives min abs root ≈ 0.8, so I would first try x = 0.8*cis(3/4*pi ≈ 2.36)

CAS> p := makelist(ithprime, 148+1,1,-1) :;
CAS> f := horner(p, x) :;
CAS> newton(f, x = 0.8*exp (3/4*pi*i))

−0.645758096347+0.483177676218*i

CAS> polar(Ans)

[0.806513599261, 2.49922322693]

Magnitude matched root squaring produced min --> this (and its conjugate) is the root to seek.

[Albert Chan]

11-20-2022, 06:13 AM

Werner Wrote:I see what you mean - roots will be closer to -1 so there's no need to try out the complex roots to the right of the origin, you just took one 3/4 of the way, and it happens to produce the right answer. What if it didn't?
Werner

Guessing way is again bisect phase angle: 7/8*pi, 5/8*pi

(11-16-2022 01:28 PM)Albert Chan Wrote:  If we consider polynomial coefficients with geometric progression:

f(x) = 1 + (r*x) + (r*x)^2 + ... + (r*x)^n = ((r*x)^(n+1) - 1) / ((r*x) - 1)

From RHS numerator, all f roots abs = 1/r

P(x) = 2 + 3x + 5x^2 + 7x^3 + 11x^4 + ...

P roots, sorted in abs: (-0.6458±0.4832i), (0.4472±0.7248i), (-0.2853±0.8292i), ...
With corresponding abs: 0.8065, 0.8517, 0.8769, ...

Primes does not grow as fast as geometric progression. (sum of reciprocal primes diverges)

P min abs root is due to ratio, (5/2=2.5) > (11/5=2.2)
If the ratios were about the same, we expected f(x) roots pattern, with similar sized roots.

Example, R(x) = P(x)+0.5, 5/(2+0.5) = 2.

R roots, sorted in abs: (-0.6811±0.5122i), (0.4692±0.7114i), (-0.2905±0.8666i), ...
With corresponding abs: 0.8522, 0.8522, 0.9140, ...

If guess close enough, we can use Newton's method to zeroed in P min abs root.

More direct way is to try smaller degree polynomial, say, all primes under 100.
Roots in polar form, sorted in abs, with root that "belongs" to 2 and 2.5 bolded.
(2nd proot needed, to confirm location of root that really "belongs" to 2)

proot(2.0 + 3x + 5x^2 + ... + 97*x^24) = 0.788∠±2.50, 0.810∠±1.03, ...
proot(2.5 + 3x + 5x^2 + ... + 97*x^24) = 0.815∠±1.01, 0.818∠±2.50, ...

--> P(x) min abs root around 0.8∠±2.5

[Albert Chan]

11-20-2022, 10:52 AM

Albert Chan Wrote:proot(2.0 + 3x + 5x^2 + ... + 97*x^24) = 0.788∠±2.50, 0.810∠±1.03, ...
proot(2.5 + 3x + 5x^2 + ... + 97*x^24) = 0.815∠±1.01, 0.818∠±2.50, ...

--> P(x) min abs root around 0.8∠±2.5

To create a minimum degree polynomial, we need more accurate guess.

We know primes does not grow as fast as geometric series.
Assume it does not grow at all, we have:

P remaining terms = 101*x^25 + 103*x^26 + ... ≈ 102*x^25 * (x^9999-1)/(x-1)

XCas> f := r2e(makelist(ithprime, 25,1,-1)) :; // above quoted polynomial
XCas> c := 102*x^25 * (x^9999-1) / (x-1) :; // correction
XCas> polar(newton(f+c, x=0.8*exp(2.5*i)))

[0.806496, 2.50034]

We now have 4 digits accuracy.
For reference, rounded to 6 digits, P min abs root = 0.806514 ∠ 2.49922

[Albert Chan]

11-20-2022, 11:31 AM

Albert Chan Wrote:XCas> c := 102*x^25 * (x^9999-1) / (x-1) :; // correction

It is simpler to assume correction with infinite degress. We will get same result.

XCas> c := 102*x^25 / (1-x) :; // correction, infinite degrees.

[Albert Chan]

11-20-2022, 02:15 PM

We can again apply correction term to minimum degree polynomial.
For 12 digits accuracy, this reduced primes needed from 149, down to 119

n degree f + correction = f + p(n+2)*x^(n+1) / (1-x)

To simplify, assume prime gap = 0 --> p(n+2) = p(n+1)
Multiply by (1-x), x^(n+1) term disappears.

We get back n degree polynomial, with prime gap coefficents.
(if we assume 0-th prime = 0)

XCas> p := makelist(ithprime,119,1,-1);
XCas> g := expand(r2e(p - p[1:])) :;
XCas> polar(newton(g, x=0.8*exp(2.5*i)))

[0.806513599261, 2.49922322692]

[Albert Chan]

11-23-2022, 10:13 AM

Werner Wrote:- the Graeffe squaring process quickly overflows; in this case, you obtained 12+ digits of the correct answer - but you cannot know that unless you performed one or more additional steps, which would overflow. I tried it out with smaller degrees (in excel..) and I could only confirm about 6 digits before overflow. And then remains the Newton-Rhapson step to refine the value, but all I have is the abs value, not the estimate of the root, so I have n possible complex roots to verfiy. The link to Cleve Moler's article says just that. (thanks for that link BTW!)

We may get answer directly with plain root squaring.
This post may be relevant.

(11-20-2022 09:11 PM)Albert Chan Wrote:  I think the bottleneck is PROOT. What is its time complexity?
For example, below trick will reduce primes needed. How much time will it save?

Primes grow slower than any geometric series. Asymptotically, p_n ≈ n*ln(n)

Let M = minimum absolute root function

M(2+3x+5x^2+...+887*x^153) ≈ M((2+3x+5x^2+...+683*x^123) + 683*x^124/(1-x))

Last term assumed primes stop growing after p₁₂₄ = 683, thus the factor 1/(1-x)
Multiply by (1-x), we have another polynomial that contains the same min abs root, but with less terms.

2 + (3-2)*x + (5-3)*x^2 + (7-5)*x^3 + ... + (683-677)*x^123

Relative to primes, prime gap is tiny, allowing for more squarings.
To delay overflow a bit longer, work with halved (or more) of the prime gap polynomial.


Correction: to get 12 digits accuracy, we need 1 more prime.

2 + (3-2)*x + (5-3)*x^2 + (7-5)*x^3 + ... + (683-677)*x^123 + (691-683)*x^124

[Albert Chan]

11-23-2022, 11:10 AM

The trick does not work with all polynomial in general, only for VA ROOT P(x)
P(x) roots are all inside the unit disc., with P min abs root much smaller than 1

Let M = minimum abs root function

M(P(x)) = M(P(x)*(1-x)) = M(2 + (3-2)*x + (5-3)*x^2 + (7-5)*x^3 + ...)

We can reduce size of coefficients by 2 (or more), with 2nd order forward difference (see my next PM)

M(P(x)) = M(P(x)*(1-x)²) = M(2 + ((3-2)-2)*x + ((5-3)-(3-2))*x^2 + ((7-5)-(5-3))*x^3 + ...)

[Albert Chan]

11-23-2022, 11:17 AM

If we don't mind negative coefficients, we can assume prime gap don't grow, and do gap of gap.
Maximum size of coefficients will drop by 2, or more.

XCas> p := makelist(ithprime,124+1,1,-1)      → [691,683,677,673,661,659,...]
XCas> p -= p[1:]      → [8,6,4,12,2,6,6,4,2,10,12,...]
XCas> p -= p[1:]      → [2,2,-8,10,-4,0,2,2,-8,-2,10,...]
XCas> min(abs(proot(p)))      → 0.806513599261


Assumption is not as good though (last prime gap we picked may be huge)
Above example, assumption is very good.

p125-p124 = 691-683 = 8 (which, we assumed the same for all future prime gaps)
p126-p125 = 701-691 = 10
p127-p126 = 709-701 = 8
p128-p127 = 719-709 = 10
p129-p128 = 727-719 = 8
...

[Albert Chan]

11-23-2022, 01:26 PM

Werner Wrote:

Albert Chan Wrote:We may get answer directly with plain root squaring.
This post may be relevant.

Oh no, it is certainly relevant, and may *aid* in determining the answer. It is the only method that will give us an idea of what the real min abs root is. You got 12 correct digits after all - but only because you knew the answer already? To verify the correctness, you would have to square further once or twice, but you'd run into overflow.

Nice to know prime gap trick work

This is a quick and dirty way to do root powering. (symbolically, thus slow)
BTW, I had trouble translate a[::n] to HP Prime Cas.

Code:

rootpow(a,n) := { local w;
     w := exp(2*pi/n*i);
     a := symb2poly(product(horner(a,x*w^k),k=0..n-1));
     return normal(a[::n])
}

Below is for illustration purpose only.
We could use integer coefficients without worrying about overflow.
Or, numerically by XCas mpfr float, with exponents range into billions.

XCas> p := makelist(ithprime,124+1,1,-1)      → [691,683,677,673,661,...]
XCas> k, q := 2, reverse(p - p[1:]) .* .25       → 2, [0.5,0.25,0.5,0.5,1.0,...]
XCas> q := rootpow(q, 2) :; k*=2; surd(q[0]/q[2],k), surd(q[2]/q[4],k)

Repeat last expression until k=1024, we have:

Code:

4   , 0.707106781187, 0.816496580928
8   , 0.725699929618, 0.874520449757
16  , 0.812448915479, 0.827434355334
32  , 0.79583362162 , 0.901285129669
64  , 0.80643033679 , 0.852976123366
128 , 0.805938612598, 0.850214544854
256 , 0.806514291756, 0.851660093256
512 , 0.80651360183 , 0.851658832026
1024, 0.806513599261, 0.851662531880

Final check for min abs root does not require full root squaring.
For next P 2 minimum roots, we need 4*2+1 = 9 columns
(same idea for the next after, we need 8*2+1 = 17 columns)

XCas> q := q[:9]
XCas> q := rootpow(q, 2) :; k*=2; surd(q[0]/q[2],k), surd(q[2]/q[4],k)

2048, 0.806513599261, 0.85166253289      # P min abs root, and 2nd min

For rough estimate, P 2nd min is probably more important than actual min.
If we get a root with abs below its 2nd min, we've locked to min root.
All is left is to improve the root, to required accuracy.

[EdS2]

Thanks for sharing - an interesting investigation!

[Werner]

In Excel, we can verify that the result is correct to 14 digits after 9 iterations, by doing the 10th ;-)
(barely avoiding overflow - needed to scale the coefficients by 2 to achieve that btw)
   

Seeing that the 42S' range is slightly larger than Excel's, we can thus use Graeffe's method there, successively doubling the number of primes used, for instance.

We can save time on overflow checks by scaling the domain of the polynomial by 1.2 or so:

1 + 1.5x + 2.5x^2 + 3.5x^3 ... = 0

take y=x*s, s=1.2 then:

1 + 1.5*(y/s) + 2.5*(y/s)^2 + 3.5*(y/s)^3 + ... = 0

The answer we seek after n iterations is then

(b0/b2)^(1/(2^(n+1))/s

We now get underflows instead of overflows, and can disregard them, and I can even do an 11th iteration. Ok, I did a bit of trial and error there ;-)

Cheers, Werner

[Albert Chan]

(11-23-2022 07:10 PM)Albert Chan Wrote:  P(x) = 2 + 3x + 5x^2 + 7x^3 + 11x^4 + ...

All P roots inside complex unit circle.

We can assume min root real part closer to -1 than +1, to keep imag part small.
Root squaring process gives min abs root ≈ 0.8, so I would first try x = 0.8*cis(3/4*pi ≈ 2.36)

Taylor series:

P(x+y*i) = P(x) + P'(x)*(y*i) + P''(x)*(y*i)^2/2! + P'''(x)*(y*i)^3/3! + ...

P with positive coefficients, if x>0, we have |P(x)| > |P(-x)|

My intution is that if x>0, it take bigger imaginery part "correction" to get P(x+y*i) = 0
This implied min abs root is in Q2 (or Q3 for conjugate root)

This is an educated guess.

[Albert Chan]

Another way to estimate location of min abs root, start with something simple.

2 + 3x = 0 --> x = -2/3, which suggested P min abs closer to -2/3 than 2/3

(2+3x)*(1+2.5*x^2) = 2 + 3x + 5x^2 + 7.5x^3 = 0 --> x = [-2/3, ±1√(-2.5)] ≈ [-0.667, ±0.632i]

7.5 ≈ 7, which suggested P min root and 2nd min root phase angle differs by about pi/2.
The 2 min roots should be far apart, in the complex plane.

Both suggestions turned out to be true, even with small degree polynomial

(11-23-2022 07:12 PM)Albert Chan Wrote:  proot(2.0 + 3x + 5x^2 + ... + 97*x^24) = 0.788∠±2.50, 0.810∠±1.03, ...

And, P with infinite degrees:

proot(2 + 3x + 5x^2 + ... + 97*x^24 + ...) = 0.807∠±2.50, 0.852∠±1.02, ...