diff --git a/src/Math-Numerical/PMNewtonZeroFinder.class.st b/src/Math-Numerical/PMNewtonZeroFinder.class.st
index a253bdd3..9f5f2303 100644
--- a/src/Math-Numerical/PMNewtonZeroFinder.class.st
+++ b/src/Math-Numerical/PMNewtonZeroFinder.class.st
@@ -39,6 +39,15 @@ PMNewtonZeroFinder class >> function: aBlock1 derivative: aBlock2 [
 	^(self new) setFunction: aBlock1; setDerivative: aBlock2; yourself
 ]
 
+{ #category : #information }
+PMNewtonZeroFinder class >> with: precision [
+
+	| rootFinder |
+	rootFinder := self new.
+	rootFinder desiredPrecision: precision.
+	^ rootFinder
+]
+
 { #category : #operation }
 PMNewtonZeroFinder >> computeInitialValues [
 	"Private - If no derivative has been defined, take an ad-hoc definition.
diff --git a/src/Math-Polynomials/PMPolynomial.class.st b/src/Math-Polynomials/PMPolynomial.class.st
index d77820df..d941313e 100644
--- a/src/Math-Polynomials/PMPolynomial.class.st
+++ b/src/Math-Polynomials/PMPolynomial.class.st
@@ -233,27 +233,27 @@ PMPolynomial >> reciprocal [
 
 { #category : #information }
 PMPolynomial >> roots [
-
-	^ self roots: Float defaultComparisonPrecision
+		^ (self roots: Float defaultComparisonPrecision) asSortedCollection asArray
 ]
 
 { #category : #information }
 PMPolynomial >> roots: aNumber [
 
-	| pol roots x rootFinder |
-	rootFinder := PMNewtonZeroFinder new.
-	rootFinder desiredPrecision: aNumber.
-	pol := self class coefficients: ( coefficients reverse collect: [ :each | each asFloat]).
+	| pol roots root rootFinder |
+	rootFinder := PMNewtonZeroFinder with: aNumber.
+	pol := self class coefficients:
+		       (coefficients reverse collect: [ :each | each asFloat ]).
 	roots := OrderedCollection new: self degree.
-	[ rootFinder setFunction: pol; setDerivative: pol derivative.
-	  x := rootFinder evaluate.
-	  rootFinder hasConverged
-		] whileTrue: [ roots add: x.
-					   pol := pol deflatedAt: x.
-					   pol degree > 0
-						 ifFalse: [ ^roots].
-					 ].
-	^roots
+	[
+	rootFinder
+		setFunction: pol;
+		setDerivative: pol derivative.
+	root := rootFinder evaluate.
+	rootFinder hasConverged ] whileTrue: [
+		roots add: root.
+		pol := pol deflatedAt: root.
+		pol degree > 0 ifFalse: [ ^ roots ] ].
+	^ roots
 ]
 
 { #category : #'double dispatching' }
diff --git a/src/Math-Tests-Polynomials/PMPolynomialTest.class.st b/src/Math-Tests-Polynomials/PMPolynomialTest.class.st
index 3b5a2fc9..aa20322a 100644
--- a/src/Math-Tests-Polynomials/PMPolynomialTest.class.st
+++ b/src/Math-Tests-Polynomials/PMPolynomialTest.class.st
@@ -235,13 +235,36 @@ PMPolynomialTest >> testPolynomialRoots [
 
 	| polynomial roots |
 	polynomial := PMPolynomial coefficients: #( -10 -13 -2 1 ).
-	roots := polynomial roots asSortedCollection asArray.
+	roots := polynomial roots .
 	self assert: roots size equals: 3.
 	self assert: (roots at: 1) + 2 closeTo: 0.
 	self assert: (roots at: 2) + 1 closeTo: 0.
 	self assert: (roots at: 3) - 5 closeTo: 0
 ]
 
+{ #category : #'iterative algorithms' }
+PMPolynomialTest >> testPolynomialRootsConstantsHaveNoRoots [
+
+	| constant |
+	"Here, compute the roots of the constant C = 1"
+	constant := PMPolynomial coefficients: #( 1 ).
+	self
+		should: [ constant roots ]
+		raise: Error
+		description: 'Function''s derivative seems to be zero everywhere'
+]
+
+{ #category : #'iterative algorithms' }
+PMPolynomialTest >> testPolynomialRootsForLinear [
+
+	| linearPolynomial roots |
+	"Here, compute the roots of the linear (2x + 1)"
+	linearPolynomial := PMPolynomial coefficients: #( 1 2 ).
+	roots := linearPolynomial roots.
+	self assert: roots size equals: 1.
+	self assert: (roots at: 1) closeTo: -0.5
+]
+
 { #category : #'function evaluation' }
 PMPolynomialTest >> testPolynomialSubtraction [
 	| polynomial |
@@ -253,3 +276,14 @@ PMPolynomialTest >> testPolynomialSubtraction [
 	self assert: (polynomial at: 3) equals: -1.
 	self assert: (polynomial at: 4) equals: 0
 ]
+
+{ #category : #'iterative algorithms' }
+PMPolynomialTest >> testPolynomialWithRepeatedRoots [
+	| polynomialWithRepeatedRoots roots |
+	"Here, compute the roots of the quadratic (2x + 1)^2 = 4 x^2 + 4 x + 1"
+	polynomialWithRepeatedRoots := PMPolynomial coefficients: #(1 4 4).
+	roots := polynomialWithRepeatedRoots roots .
+	self assert: roots size equals: 2.
+	self assert: (roots at: 1) closeTo: -0.5 .
+	self assert: (roots at: 2) closeTo: -0.5 .
+]