private PolynomialRoots()
    {
        // Do not instantiate
    }

    /**
     * Emulate the F77 sign function.
     * @param a double; the value to optionally sign invert
     * @param b double; the sign of which determines what to do
     * @return double; if b >= 0 then a; else -a
     */
    private static double sign(final double a, final double b)
    {
        return b >= 0 ? a : -a;
    }

    /**
     * LINEAR POLYNOMIAL ROOT SOLVER.
     * <p>
     * Calculates the root of the linear polynomial:<br>
     * q1 * x + q0<br>
     * Unlike the quadratic, cubic and quartic code, this is NOT derived from that Fortran90 code; it was added for completenes.
     * @param q1 double; coefficient of the x term
     * @param q0 double; independent coefficient
     * @return Complex[]; the roots of the equation
     */
    public static Complex[] linearRoots(final double q1, final double q0)
    {
        if (q1 == 0)
        {
            return new Complex[] {}; // No roots; return empty array
        }
        return linearRoots(q0 / q1);
    }

    /**
     * LINEAR POLYNOMIAL ROOT SOLVER.
     * <p>
     * Calculates the root of the linear polynomial:<br>
     * x + q0<br>
     * Unlike the quadratic, cubic and quartic code, this is NOT derived from that Fortran90 code; it was added for completenes.
     * @param q0 double; independent coefficient
     * @return Complex[]; the roots of the equation
     */
    public static Complex[] linearRoots(final double q0)
    {
        return new Complex[] { new Complex(-q0, 0) };
    }

    /**
     * QUADRATIC POLYNOMIAL ROOT SOLVER
     * <p>
     * Calculates all real + complex roots of the quadratic polynomial:<br>
     * q2 * x^2 + q1 * x + q0<br>
     * The code checks internally if rescaling of the coefficients is needed to avoid overflow.
     * <p>
     * The order of the roots is as follows:<br>
     * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
     * negative last).<br>
     * 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
     * q1 : coefficient of x term q0 : independent coefficient
     * @param q2 double; coefficient of the quadratic term
     * @param q1 double; coefficient of the x term
     * @param q0 double; independent coefficient
     * @return Complex[]; the roots of the equation
     */
    public static Complex[] quadraticRoots(final double q2, final double q1, final double q0)
    {
        if (q2 == 0)
        {
            return linearRoots(q1, q0);
        }
        return quadraticRoots(q1 / q2, q0 / q2);
    }

    /**
     * QUADRATIC POLYNOMIAL ROOT SOLVER
     * <p>
     * Calculates all real + complex roots of the quadratic polynomial:<br>
     * x^2 + q1 * x + q0<br>
     * The code checks internally if rescaling of the coefficients is needed to avoid overflow.
     * <p>
     * The order of the roots is as follows:<br>
     * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
     * negative last).<br>
     * 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
     * q1 : coefficient of x term q0 : independent coefficient
     * @param q1 double; coefficient of the x term
     * @param q0 double; independent coefficient
     * @return Complex[]; the roots of the equation
     */
    public static Complex[] quadraticRoots(final double q1, final double q0)
    {
        boolean rescale;

        double a0, a1;
        double k = 0, x, y, z;

        // Handle special cases.
        if (q0 == 0.0 && q1 == 0.0)
        {
            // Two real roots at 0,0
            return new Complex[] { Complex.ZERO, Complex.ZERO };
        }
        else if (q0 == 0.0)
        {
            // Two real roots; one of these is 0,0
            // x^2 + q1 * x == x * (x + q1)
            Complex nonZeroRoot = new Complex(-q1);
            return new Complex[] { q1 > 0 ? Complex.ZERO : nonZeroRoot, q1 <= 0 ? nonZeroRoot : Complex.ZERO };
        }
        else if (q1 == 0.0)
        {
            x = Math.sqrt(Math.abs(q0));

            if (q0 < 0.0)
            {
                // Two real roots, symmetrically around 0
                return new Complex[] { new Complex(x, 0), new Complex(-x, 0) };
            }
            else
            {
                // Two complex roots, symmetrically around 0
                return new Complex[] { new Complex(0, x), new Complex(0, -x) };
            }
        }
        else
        {
            // The general case. Do rescaling, if either squaring of q1/2 or evaluation of
            // (q1/2)^2 - q0 will lead to overflow. This is better than to have the solver
            // crashed. Note, that rescaling might lead to loss of accuracy, so we only
            // invoke it when absolutely necessary.
            final double sqrtLPN = Math.sqrt(Double.MAX_VALUE); // Square root of the Largest Positive Number
            rescale = (q1 > sqrtLPN + sqrtLPN); // this detects overflow of (q1/2)^2

            if (!rescale)
            {
                x = q1 * 0.5; // we are sure here that x*x will not overflow
                rescale = (q0 < x * x - Double.MAX_VALUE); // this detects overflow of (q1/2)^2 - q0
            }

            if (rescale)
            {
                x = Math.abs(q1);
                y = Math.sqrt(Math.abs(q0));

                if (x > y)
                {
                    k = x;
                    z = 1.0 / x;
                    a1 = sign(1.0, q1);
                    a0 = (q0 * z) * z;
                }
                else
                {
                    k = y;
                    a1 = q1 / y;
                    a0 = sign(1.0, q0);
                }
            }
            else
            {
                a1 = q1;
                a0 = q0;
            }
            // Determine the roots of the quadratic. Note, that either a1 or a0 might
            // have become equal to zero due to underflow. But both cannot be zero.
            x = a1 * 0.5;
            y = x * x - a0;

            if (y >= 0.0)
            {
                // Two real roots
                y = Math.sqrt(y);

                if (x > 0.0)
                {
                    y = -x - y;
                }
                else
                {
                    y = -x + y;
                }

                if (rescale)
                {
                    y = y * k; // very important to convert to original
                    z = q0 / y; // root first, otherwise complete loss of
                }
                else // root due to possible a0 = 0 underflow
                {
                    z = a0 / y;
                }
                return new Complex[] { new Complex(Math.max(y, z), 0), new Complex(Math.min(y, z), 0) };
            }
            else
            {
                // Two complex roots (zero real roots)
                y = Math.sqrt(-y);

                if (rescale)
                {
                    x *= k;
                    y *= k;
                }
                return new Complex[] { new Complex(-x, y), new Complex(-x, -y) };
            }
        }
    }

    /**
     * CUBIC POLYNOMIAL ROOT SOLVER.
     * <p>
     * Calculates all (real and complex) roots of the cubic polynomial:<br>
     * c3 * x^3 + c2 * x^2 + c1 * x + c0<br>
     * The first real root (which always exists) is obtained using an optimized Newton-Raphson scheme. The other remaining roots
     * are obtained through composite deflation into a quadratic.
     * <P>
     * The cubic root solver can handle any size of cubic coefficients and there is no danger of overflow due to proper
     * rescaling of the cubic polynomial. The order of the roots is as follows: 1) For real roots, the order is according to
     * their algebraic value on the number scale (largest positive first, largest negative last). 2) Since there can be only one
     * complex conjugate pair root, no order is necessary. 3) All real roots precede the complex ones.
     * @param c3 double; coefficient of the cubic term
     * @param c2 double; coefficient of the quadratic term
     * @param c1 double; coefficient of the linear term
     * @param c0 double; coefficient of the independent term
     * @return Complex[]; array of Complex with all the roots
     */
    public static Complex[] cubicRoots(final double c3, final double c2, final double c1, final double c0)

        EventProducingIterator<T> iterator = new EventProducingIterator<T>(this.parent.iterator(), this.sourceIdProvider);
        // WEAK reference as an iterator is usually local and should be eligible for garbage collection
        iterator.addListener(this, EventProducingIterator.OBJECT_REMOVED_EVENT, ReferenceType.WEAK);
        return iterator;
    }

    /** {@inheritDoc} */
    @Override
    public void notify(final EventInterface event) throws RemoteException
    {
        // pass through the OBJECT_REMOVED_EVENT from the iterator
        if (event.getType().equals(EventProducingIterator.OBJECT_REMOVED_EVENT))
        {
            this.fireEvent(OBJECT_REMOVED_EVENT, this.parent.size());
        }
    }

    /** {@inheritDoc} */
    @Override
    public boolean remove(final Object o)
    {
        boolean changed = this.parent.remove(o);
        if (changed)
        {
            this.fireEvent(OBJECT_REMOVED_EVENT, this.parent.size());
        }
        return changed;
    }

    /** {@inheritDoc} */
    @Override
    public boolean removeAll(final Collection<?> c)
    {
        boolean changed = this.parent.removeAll(c);
        if (changed)
        {
            this.fireEvent(OBJECT_REMOVED_EVENT, this.parent.size());
        }
        return changed;
    }

    /** {@inheritDoc} */
    @Override
    public boolean retainAll(final Collection<?> c)
    {
        boolean changed = this.parent.retainAll(c);
        if (changed)
        {
            this.fireEvent(OBJECT_REMOVED_EVENT, this.parent.size());
        }
        return changed;
    }

    /** {@inheritDoc} */
    @Override
    public Object[] toArray()
    {
        return this.parent.toArray();
    }

    /** {@inheritDoc} */
    @Override
    public <E> E[] toArray(final E[] a)

    public EventProducingCollection(final Collection<T> parent, final IdProvider sourceIdProvider)
    {
        Throw.whenNull(parent, "parent cannot be null");
        Throw.whenNull(sourceIdProvider, "sourceIdprovider cannot be null");
        this.parent = parent;
        this.sourceIdProvider = sourceIdProvider;
    }

    /** {@inheritDoc} */
    @Override
    public Serializable getSourceId()
    {
        return this.sourceIdProvider.id();
    }

    /** {@inheritDoc} */
    @Override
    public int size()
    {
        return this.parent.size();
    }

    /** {@inheritDoc} */
    @Override
    public boolean isEmpty()
    {
        return this.parent.isEmpty();
    }

    /** {@inheritDoc} */
    @Override
    public void clear()
    {
        int nr = this.parent.size();
        this.parent.clear();
        if (nr != this.parent.size())
        {
            this.fireEvent(OBJECT_REMOVED_EVENT, this.parent.size());
        }
    }

    /** {@inheritDoc} */
    @Override
    public boolean add(final T o)

    public EventProducingList(final List<E> parent, final IdProvider sourceIdProvider)
    {
        Throw.whenNull(parent, "parent cannot be null");
        Throw.whenNull(sourceIdProvider, "sourceIdprovider cannot be null");
        this.parent = parent;
        this.sourceIdProvider = sourceIdProvider;
    }

    /** {@inheritDoc} */
    @Override
    public Serializable getSourceId()
    {
        return this.sourceIdProvider.id();
    }

    /** {@inheritDoc} */
    @Override
    public int size()
    {
        return this.parent.size();
    }

    /** {@inheritDoc} */
    @Override
    public boolean isEmpty()
    {
        return this.parent.isEmpty();
    }

    /** {@inheritDoc} */
    @Override
    public void clear()
    {
        int nr = this.parent.size();
        this.parent.clear();
        if (nr != this.parent.size())
        {
            this.fireEvent(OBJECT_REMOVED_EVENT, this.parent.size());
        }
    }

    /** {@inheritDoc} */
    @Override
    public void add(final int index, final E element)

    public EventProducingCollection(final Collection<T> parent, final IdProvider sourceIdProvider)
    {
        Throw.whenNull(parent, "parent cannot be null");
        Throw.whenNull(sourceIdProvider, "sourceIdprovider cannot be null");
        this.parent = parent;
        this.sourceIdProvider = sourceIdProvider;
    }

    /** {@inheritDoc} */
    @Override
    public Serializable getSourceId()
    {
        return this.sourceIdProvider.id();
    }

    /** {@inheritDoc} */
    @Override
    public int size()
    {
        return this.parent.size();
    }

    /** {@inheritDoc} */
    @Override
    public boolean isEmpty()
    {
        return this.parent.isEmpty();
    }

    /** {@inheritDoc} */
    @Override
    public void clear()
    {
        int nr = this.parent.size();
        this.parent.clear();
        if (nr != this.parent.size())
        {
            this.fireEvent(OBJECT_REMOVED_EVENT, this.parent.size());
        }
    }

    /** {@inheritDoc} */
    @Override
    public boolean add(final T o)

        String format = String.format("%%%ds", fieldsPerLine + (fieldsPerLine - 1) / this.extraSpaceAfterEvery);
        this.prototypeLine = String.format(format, "");
    }

    /** String builder for current output line. */
    private StringBuilder buffer = new StringBuilder();

    /** {@inheritDoc} */
    @Override
    public String getResult()
    {
        String result = this.buffer.toString();
        this.buffer.setLength(0);
        return result;
    }

    /** {@inheritDoc} */
    @Override
    public int getMaximumWidth()
    {
        return this.prototypeLine.length();
    }

    /** {@inheritDoc} */
    @Override
    public boolean append(final int address, final byte theByte) throws IOException
    {
        if (this.buffer.length() == 0)
        {
            this.buffer.append(this.prototypeLine);
        }
        int lineByte = address % this.fieldsPerLine;
        int index = lineByte + lineByte / this.extraSpaceAfterEvery;

        StackTraceElement[] ste = steList.toArray(new StackTraceElement[steList.size()]);
        String where = ste[0].getClassName() + "." + ste[0].getMethodName() + " (" + ste[0].getLineNumber() + "): ";
        Object[] args = argList.toArray();
        String formattedMessage;
        try
        {
            formattedMessage = where + String.format(message, args);
        }
        catch (IllegalFormatException exception)
        {
            formattedMessage = where + message + " [FormatException; args=" + argList + "]";
        }

        // throw all other exceptions through reflection
        T exception;
        try
        {
            Constructor<T> constructor = ClassUtil.resolveConstructor(throwableClass, new Class<?>[] {String.class});

    public boolean addAll(final Collection<? extends T> c)
    {
        boolean changed = this.parent.addAll(c);
        if (changed)
        {
            this.fireEvent(OBJECT_ADDED_EVENT, this.parent.size());
        }
        else
        {
            if (!c.isEmpty())
            {
                this.fireEvent(OBJECT_CHANGED_EVENT, this.parent.size());
            }
        }
        return changed;
    }

    /** {@inheritDoc} */
    @Override
    public boolean contains(final Object o)
    {
        return this.parent.contains(o);
    }

    /** {@inheritDoc} */
    @Override
    public boolean containsAll(final Collection<?> c)
    {
        return this.parent.containsAll(c);
    }

    /** {@inheritDoc} */
    @Override
    public EventProducingIterator<T> iterator()

    public static Complex[] rootsDurandKerner(final Complex[] a)
    {
        int n = a.length - 1;
        double radius = 1 + maxAbs(a);

        // initialize the initial values, not as a real number and not as a root of unity
        Complex[] p = new Complex[n];
        p[0] = new Complex(Math.sqrt(radius), Math.cbrt(radius));
        double rot = 350.123 / n;
        for (int i = 1; i < n; i++)
        {
            p[i] = p[0].rotate(rot * i);
        }

        double maxError = 1.0;
        int count = 0;
        while (maxError > 0 && count < MAX_STEPS_DURAND_KERNER)