Polynomial square root

Jane Sullivan

unread,

Aug 31, 2012, 8:49:09 AM8/31/12

to

I have this function to extract the square root of a polynomial. A
polynomial is represented as a vector of the coefficients. If the
polynomial does not have a square root, a value is returned, which may
or may not be useful, because it is up to the user to ensure that the
input is a square.

∇ z←SQRT x;i;n;q;r
⍝ Find the square root of a Polynomial
⎕SIGNAL(1≠2|⍴x)/11 ⍝ Must be an odd number of terms
z←,x[0]*0.5 ⍝ First term of the result
r←1↓x ⍝ First Remainder
q←2×z
:For i :In 1+⍳0.5×¯1+⍴x
z,←n←r[0]÷q[0] ⍝ Next term in the result
r←1↓r-(⍴r)↑n×q←q,n ⍝ Next remainder
q←q+(-⍴q)↑n
:EndFor
∇

SQRT 1 ¯2 ¯1 0 3 2 1
1 ¯1 ¯1 ¯1
SQRT 1 ¯4 10 ¯16 19 ¯16 10 ¯4 1
1 ¯2 3 ¯2 1
SQRT 9 6 ¯11 ¯4 4
3 1 ¯2

Can anyone think of any improvements that can be made to it?

--
Jane

James J. Weinkam

unread,

Aug 31, 2012, 5:04:39 PM8/31/12

to

Jane Sullivan wrote:
> I have this function to extract the square root of a polynomial. A polynomial is represented as a vector of the

> coefficients. If the polynomial does not have a square root, a value isreturned, which may or may not be useful,

> because it is up to the user to ensure that the input is a square.
>
> ∇ z←SQRT x;i;n;q;r
> ⍝ Find the square root of a Polynomial
> ⎕SIGNAL(1≠2|⍴x)/11 ⍝ Must be an odd number of terms
> z←,x[0]*0.5 ⍝ First term of the result
> r←1↓x ⍝ First Remainder
> q←2×z
> :For i :In 1+⍳0.5×¯1+⍴x
> z,←n←r[0]÷q[0] ⍝ Next term in the result
> r←1↓r-(⍴r)↑n×q←q,n ⍝Next remainder
> q←q+(-⍴q)↑n
> :EndFor
> ∇
>
> SQRT 1 ¯2 ¯1 0 3 2 1
> 1 ¯1 ¯1 ¯1
> SQRT 1 ¯4 10 ¯16 19 ¯16 10 ¯4 1
> 1 ¯2 3 ¯2 1
> SQRT 9 6 ¯11 ¯4 4
> 3 1 ¯2
>
> Can anyone think of any improvements that can be made to it?
>

Here's a version in APL2 that uses a different approach: Consider a polynomial, B, of degree 2k with coefficients b. If
B is a perfect square, its square root is a polynomial A of degree k with unknown coefficients a. If we symbolically
square A, we obtain 2k+1 simultaneous second degree equations in the k+1 variables a. The equation for b[0] only
involves a[0]; that for b[1] only involves a[0] and a[1], etc up to k. (A similar situation holds true working backwards
from 2k+1 to k.)

My version gives results identical to yours. No surprise. Timing comparisons are inconclusive.

[0] a←psqrt b;c;d;i;n;⎕IO
[1] ⎕IO←0 ◊ a←(1+n←⌊(⍴b)÷2)⍴0 ◊ d←2×a[0]←b[0]*0.5 ◊ →(n=0)/0 ◊ i←1
[2] lp:a[i]←(b[i]-+/(⌽c)×c←1↓i↑a)÷d ◊ i←i+1 ◊ →(i≤n)/lp

James J. Weinkam

unread,

Aug 31, 2012, 5:06:48 PM8/31/12

to

That should have been 2k.

Roger Hui

unread,

Aug 31, 2012, 7:32:37 PM8/31/12

to

The obligatory one-liner (in Dyalog APL):

psqrt2←{⍵{⍵,(2×1↑⍵)÷⍨⍺[⍴⍵]-c+.×⌽c←1↓⍵}⍣(⌊2÷⍨⍬⍴⍴⍵)⊢0.5*⍨1↑⍵}

psqrt2 1 ¯2 ¯1 0 3 2 1
1 ¯1 ¯1 ¯1
psqrt2 1 ¯4 10 ¯16 19 ¯16 10 ¯4 1

1 ¯2 3 ¯2 1

psqrt2 9 6 ¯11 ¯4 4
3 1 ¯2
psqrt2 ,17
4.12311

psq←{+⌿(-⍳⍴⍵)⌽⍵∘.×⍵,1↓(⍴⍵)⍴0}

psqrt2 psq 3 1 4 ¯1 5 9.26
3 1 4 ¯1 5 9.26