static public BigDecimal sin(final BigDecimal x)
        {
                if ( x.compareTo(BigDecimal.ZERO) < 0)
                        return sin(x.negate()).negate() ;
                else if ( x.compareTo(BigDecimal.ZERO) == 0 )
                        return BigDecimal.ZERO ;
                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: sin(x)= - sin(x-pi)
                                */
                                return sin(subtractRound(res,p)) .negate() ;
                        }
                        else if ( res.multiply(new BigDecimal("2")).compareTo(p) > 0 )
                        {
                                /* pi/2<x<=pi: sin(x)= sin(pi-x)
                                */
                                return sin(subtractRound(p,res)) ;
                        }
                        else
                        {
                                /* for the range 0<=x<Pi/2 one could use sin(2x)=2sin(x)cos(x)
                                * to split this further. Here, use the sine up to pi/4 and the cosine higher up.
                                */
                                if ( res.multiply(new BigDecimal("4")).compareTo(p) > 0 )
                                {
                                        /* x>pi/4: sin(x) = cos(pi/2-x)
                                        */
                                        return cos( subtractRound(p.divide(new BigDecimal("2")),res) ) ;
                                }
                                else
                                {
                                        /* Simple Taylor expansion, sum_{i=1..infinity} (-1)^(..)res^(2i+1)/(2i+1)! */
                                        BigDecimal resul = res ;

                                        /* x^i */
                                        BigDecimal xpowi = res ;

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

                                        /* The error in the result is set by the error in x itself.
                                        */
                                        double xUlpDbl = res.ulp().doubleValue() ;

                                        /* The error in the result is set by the error in x itself.
                                        * We need at most k terms to squeeze x^(2k+1)/(2k+1)! below this value.
                                        * x^(2k+1) < x.ulp; (2k+1)*log10(x) < -x.precision; 2k*log10(x)< -x.precision;
                                        * 2k*(-log10(x)) > x.precision; 2k*log10(1/x) > x.precision
                                        */
                                        int k = (int)(res.precision()/ Math.log10(1.0 / res.doubleValue()))/2 ;
                                        MathContext mcTay = new MathContext( err2prec(res.doubleValue(),xUlpDbl/k) ) ;
                                        for(int i=1 ; ; i++)
                                        {
                                                /* TBD: at which precision will 2*i or 2*i+1 overflow? 
                                                */
                                                ifac = ifac.multiply(new BigInteger(""+(2*i) ) ) ;
                                                ifac = ifac.multiply( new BigInteger(""+(2*i+1)) ) ;
                                                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 set by the error in x itself.
                                        */
                                        mc = new MathContext(res.precision() ) ;
                                        return resul.round(mc) ;
                                }
                        }
                }
        } /* sin */