public Vector<BigComplex> roots(int digits)
        {
                RatPoly mon = monic() ;

                Random rand = new Random() ;
                MathContext mc = new MathContext(digits+3, RoundingMode.DOWN) ;

                Vector<BigComplex> res =new Vector<BigComplex>() ;

                final int d = mon.degree() ;
                double randRad =0. ;
                for(int i=0 ; i <= d ; i++)
                {
                        /* scale coefficient at maximum degree */
                        double absi = Math.abs(mon.at(i).doubleValue()) ;
                        if ( absi > randRad) 
                                randRad = absi ;
                }
                randRad += 1.0 ;

                /* initial values randomly in radius 1+randRad
                */
                for(int i=0 ; i < d ; i++)
                {
                        double rad = randRad*rand.nextDouble() ;
                        double phi = 2.0*3.14159*rand.nextDouble() ;
                        res.add(i, new BigComplex(rad* Math.cos(phi),rad* Math.sin(phi)) ) ;
                }

                /* iterate until convr indicates that all values changed by less than the digits 
                * precision indicates.
                */
                boolean convr = false ;
                for(int itr =0 ; ! convr ; itr++)
                {
                        convr = true ;
                        Vector<BigComplex> resPlus =new Vector<BigComplex>() ;
                        for(int v=0 ; v < d; v++)
                        {
                                /* evaluate f(x)/(x-root1)/(x-root2)/... (x-rootdegr), Newton method
                                */
                                BigComplex thisx = res.elementAt(v) ;
                                BigComplex nv = mon.valueOf(thisx,mc) ;
                                for(int j=0; j < d ; j++)
                                {
                                        if ( j != v )
                                                nv = nv.divide(thisx.subtract(res.elementAt(j)),mc) ;
                                }

                                /* is this value converged ? */
                                if ( nv.abs(mc).doubleValue() >  thisx.abs(mc).doubleValue()* Math.pow(10.0, -digits) )
                                        convr =false;
        
                                thisx = thisx.subtract(nv) ;

                                /* If unstable, start over */
                                if ( thisx.abs(MathContext.DECIMAL32).doubleValue() > randRad )
                                        return roots(digits) ;

                                resPlus.add(thisx) ;
                        }
                        res = resPlus ;

                }
                return res;
        } /* roots */