8000 ENH: Improve the computation of polynomials from roots. · numpy-buildbot/numpy@ef767a5 · GitHub
[go: up one dir, main page]
Skip to content < 8000 script type="application/json" data-target="react-partial.embeddedData">{"props":{"docsUrl":"https://docs.github.com/get-started/accessibility/keyboard-shortcuts"}}

Commit ef767a5

Browse files
charrischarris
committed
ENH: Improve the computation of polynomials from roots.
The original method was overly sensitive to roundoff. Of the two approaches considered, gauss integration or binary subdivision of the roots, the latter is more compatible with using other number representations such as mpmath. No method is going to be suitable for large numbers of arbitrary zeros but the current method is a significant improvement.
1 parent 13010c7 commit ef767a5

File tree

6 files changed

+66
-24
lines changed

6 files changed

+66
-24
lines changed

‎numpy/polynomial/chebyshev.py

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -540,10 +540,17 @@ def chebfromroots(roots) :
540540
return np.ones(1)
541541
else :
542542
[roots] = pu.as_series([roots], trim=False)
543-
prd = np.array([1], dtype=roots.dtype)
544-
8000 for r in roots:
545-
prd = chebsub(chebmulx(prd), r*prd)
546-
return prd
543+
roots.sort()
544+
n = len(roots)
545+
p = [chebline(-r, 1) for r in roots]
546+
while n > 1:
547+
m = n//2
548+
tmp = [chebmul(p[i], p[i+m]) for i in range(m)]
549+
if n%2:
550+
tmp[0] = chebmul(tmp[0], p[-1])
551+
p = tmp
552+
n = m
553+
return p[0]
547554

548555

549556
def chebadd(c1, c2):

‎numpy/polynomial/hermite.py

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -287,10 +287,17 @@ def hermfromroots(roots) :
287287
return np.ones(1)
288288
else :
289289
[roots] = pu.as_series([roots], trim=False)
290-
prd = np.array([1], dtype=roots.dtype)
291-
for r in roots:
292-
prd = hermsub(hermmulx(prd), r*prd)
293-
return prd
290+
roots.sort()
291+
n = len(roots)
292+
p = [hermline(-r, 1) for r in roots]
293+
while n > 1:
294+
m = n//2
295+
tmp = [hermmul(p[i], p[i+m]) for i in range(m)]
296+
if n%2:
297+
tmp[0] = hermmul(tmp[0], p[-1])
298+
p = tmp
299+
n = m
300+
return p[0]
294301

295302

296303
def hermadd(c1, c2):

‎numpy/polynomial/hermite_e.py

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -287,10 +287,17 @@ def hermefromroots(roots) :
287287
return np.ones(1)
288288
else :
289289
[roots] = pu.as_series([roots], trim=False)
290-
prd = np.array([1], dtype=roots.dtype)
291-
for r in roots:
292-
prd = hermesub(hermemulx(prd), r*prd)
293-
return prd
290+
roots.sort()
291+
n = len(roots)
292+
p = [hermeline(-r, 1) for r in roots]
293+
while n > 1:
294+
m = n//2
295+
tmp = [hermemul(p[i], p[i+m]) for i in range(m)]
296+
if n%2:
297+
tmp[0] = hermemul(tmp[0], p[-1])
298+
p = tmp
299+
n = m
300+
return p[0]
294301

295302

296303
def hermeadd(c1, c2):

‎numpy/polynomial/laguerre.py

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -283,10 +283,17 @@ def lagfromroots(roots) :
283283
return np.ones(1)
284284
else :
285285
[roots] = pu.as_series([roots], trim=False)
286-
prd = np.array([1], dtype=roots.dtype)
287-
for r in roots:
288-
prd = lagsub(lagmulx(prd), r*prd)
289-
return prd
286+
roots.sort()
287+
n = len(roots)
288+
p = [lagline(-r, 1) for r in roots]
289+
while n > 1:
290+
m = n//2
291+
tmp = [lagmul(p[i], p[i+m]) for i in range(m)]
292+
if n%2:
293+
tmp[0] = lagmul(tmp[0], p[-1])
294+
p = tmp
295+
n = m
296+
return p[0]
290297

291298

292299
def lagadd(c1, c2):

‎numpy/polynomial/legendre.py

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -314,10 +314,17 @@ def legfromroots(roots) :
314314
return np.ones(1)
315315
else :
316316
[roots] = pu.as_series([roots], trim=False)
317-
prd = np.array([1], dtype=roots.dtype)
318-
for r in roots:
319-
prd = legsub(legmulx(prd), r*prd)
320-
return prd
317+
roots.sort()
318+
n = len(roots)
319+
p = [legline(-r, 1) for r in roots]
320+
while n > 1:
321+
m = n//2
322+
tmp = [legmul(p[i], p[i+m]) for i in range(m)]
323+
if n%2:
324+
tmp[0] = legmul(tmp[0], p[-1])
325+
p = tmp
326+
n = m
327+
return p[0]
321328

322329

323330
def legadd(c1, c2):

‎numpy/polynomial/polynomial.py

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -186,10 +186,17 @@ def polyfromroots(roots) :
186186
return np.ones(1)
187187
else :
188188
[roots] = pu.as_series([roots], trim=False)
189-
prd = np.array([1], dtype=roots.dtype)
190-
for r in roots:
191-
prd = polysub(polymulx(prd), r*prd)
192-
return prd
189+
roots.sort()
190+
n = len(roots)
191+
p = [polyline(-r, 1) for r in roots]
192+
while n > 1:
193+
m = n//2
194+
tmp = [polymul(p[i], p[i+m]) for i in range(m)]
195+
if n%2:
196+
tmp[0] = polymul(tmp[0], p[-1])
197+
p = tmp
198+
n = m
199+
return p[0]
193200

194201

195202
def polyadd(c1, c2):

0 commit comments

Comments
 (0)
0