static public BigDecimal cos(final BigDecimal x)
        {
                if ( x.compareTo(BigDecimal.ZERO) < 0)
                        return cos(x.negate());
                else if ( x.compareTo(BigDecimal.ZERO) == 0 )
                        return BigDecimal.ONE ;
                else
                {
                        /* reduce modulo 2pi
                        */
                        BigDecimal res = mod2pi(x) ;
                        double errpi = 0.5* Math.abs(x.ulp().doubleValue()) ;
                        MathContext mc = new MathContext( 2+err2prec(3.14159,errpi) ) ;
                        BigDecimal p= pi(mc) ;
                        mc = new MathContext( x.precision() ) ;
                        if ( res.compareTo(p) > 0 )
                        {
                                /* pi<x<=2pi: cos(x)= - cos(x-pi)
                                */
                                return cos( subtractRound(res,p)) .negate() ;
                        }
                        else if ( res.multiply(new BigDecimal("2")).compareTo(p) > 0 )
                        {
                                /* pi/2<x<=pi: cos(x)= -cos(pi-x)
                                */
                                return cos( subtractRound(p,res)).negate() ;
                        }
                        else
                        {
                                /* for the range 0<=x<Pi/2 one could use cos(2x)= 1-2*sin^2(x)
                                * to split this further, or use the cos up to pi/4 and the sine higher up.
                                        throw new ProviderException("Not implemented: cosine ") ;
                                */
                                if ( res.multiply(new BigDecimal("4")).compareTo(p) > 0 )
                                {
                                        /* x>pi/4: cos(x) = sin(pi/2-x)
                                        */
                                        return sin( subtractRound(p.divide(new BigDecimal("2")),res) ) ;
                                }
                                else
                                {
                                        /* Simple Taylor expansion, sum_{i=0..infinity} (-1)^(..)res^(2i)/(2i)! */
                                        BigDecimal resul = BigDecimal.ONE ;

                                        /* x^i */
                                        BigDecimal xpowi = BigDecimal.ONE ;

                                        /* 2i factorial */
                                        BigInteger ifac = BigInteger.ONE ;

                                        /* The absolute error in the result is the error in x^2/2 which is x times the error in x.
                                        */
                                        double xUlpDbl = 0.5*res.ulp().doubleValue()*res.doubleValue() ;

                                        /* The error in the result is set by the error in x^2/2 itself, xUlpDbl.
                                        * We need at most k terms to push x^(2k+1)/(2k+1)! below this value.
                                        * x^(2k) < xUlpDbl; (2k)*log(x) < log(xUlpDbl);
                                        */
                                        int k = (int)(Math.log(xUlpDbl)/ Math.log(res.doubleValue()) )/2 ;
                                        MathContext mcTay = new MathContext( err2prec(1.,xUlpDbl/k) ) ;
                                        for(int i=1 ; ; i++)
                                        {
                                                /* TBD: at which precision will 2*i-1 or 2*i overflow? 
                                                */
                                                ifac = ifac.multiply(new BigInteger(""+(2*i-1) ) ) ;
                                                ifac = ifac.multiply( new BigInteger(""+(2*i)) ) ;
                                                xpowi = xpowi.multiply(res).multiply(res).negate() ;
                                                BigDecimal corr = xpowi.divide(new BigDecimal(ifac),mcTay) ;
                                                resul = resul.add( corr ) ;
                                                if ( corr.abs().doubleValue() < 0.5*xUlpDbl ) 
                                                        break ;
                                        }
                                        /* The error in the result is governed by the error in x itself.
                                        */
                                        mc = new MathContext( err2prec(resul.doubleValue(),xUlpDbl) ) ;
                                        return resul.round(mc) ;
                                }
                        }
                }
        } /* BigDecimalMath.cos */