FULL PRODUCT VERSION :
would be nice if 'java -version' would send output to a file

Java version "1.8.0_31"
Java (TM) SE Runtime environment (build 1.8.0_31-b13)
Java Hotspot(TM) 64-bit server VM (build 25.31-b07, mixed mode)

ADDITIONAL OS VERSION INFORMATION :
Windows 7 Home Premium, Service Pack 1

EXTRA RELEVANT SYSTEM CONFIGURATION :
Pentium Dual-core, E5700, 3GHz, 3GB RAM, 64-bit o/s

Running on indivual PC, as Java Application, not applet, with main() program.

Developed and running under DOS-type command window using javac.eve and java.exe



A DESCRIPTION OF THE PROBLEM :
I wrote a high-precision mathematical package using the BigDecimal class.The program is 5+ years old and is being enhanced constantly. This was probably a problem but never noticed. Of course it used earlier Java versions and was OK.

I suppose there is always the chance that this is my error but I cannot account for one after cutting away 90% of the program not involved in the failure. There is no polymorphism involved and no inheritance in the test that remains. No multitasking, no re-entrancy, no threads. Just straight synchronous execution.

When computations are done to up to about 1,000 decimal places, no problem. By accident while trying to increase efficiency, the precision of a few internal constants such as PI and ln(10) reached 1,500 places. When that happened, I found that the toBigInteger() function - which I use ONLY to format the final BigDecimal output myself - produced a value that is incorrect, in one case the value being exactly 1024 times the correct integer value. Another time it was 1024*256, and another was a greater power of 2. This happened in only ONE of several thousand benchmark tests I do.

I am unable to write a simple test program that will fail, however a scaled down version of my program will demonstrate it without anyone having to understand the math being done. If I try to copy/paste the code in this form, it may be too large... it would be better to email it to you. Or I can probably post it on a website I have access to. I've been chasing this bug for 7 days and only now realized my program was actually not failing, and only the function that formats the result is failing.

* Nothing here is copyright despite remaining comments; I'm sure the Mathematicians consulted to build Java know more than I do.

STEPS TO FOLLOW TO REPRODUCE THE PROBLEM :
In file Big.java, there is a constant called 'BLOCK' ( around line 60). It is the number of internal decimal places used by the program. If the value 1000 is used the program works, if 1500 is used, toBigInteger() in function print() produces wrong result.

The program should be recompiled and run with each value of BLOCK.

The batch file b_build.bat is the file that compiles the program.
Then do "java testGamma >outputFile" with a new outputFile for each test.

The line in the output file titled "x.toBigInt()" 2 lines up from the bottom will be different; it should match the value above it: it will be OK for BLOCK=1000 but not for BLOCK=1500 decimal places.

... the function toBigInteger() is used elsewhere and doesn't seem to be causing a problem with computations. That is still mysterious.

EXPECTED VERSUS ACTUAL BEHAVIOR :
EXPECTED -
Using BLOCK=1000 decimal places,

Big.print: var = 133697194382348382096607788465183338193565456046018844247569458972750894312486998024824325724736982326851035901784425686447370049387786067826320738076260459479550508892238045973237342344575773985256941526172255477095635860216762897905781270700308727774042935703399597077388908867710719254521182992345010067122532664398709404997934194880983562596833367715.16055351862835494830471966680500714677379932507670091998463150017256424061653979665143123279464197235913988149831580084826877330412482794367376691844058054222030534294681422995212574445159349829995415414703858298758123239415158255798114431318612025641643715381154114556666760352704924571922668118751259862678920840770160523030157600369321997816305387348253423618868989450392980097888212485812221827024255592193265113890702116088056809490739338023291975033964255691434246678280437248028358066225404761079657847575830923796458672560700372708577233099225155755307675956899669127488545144371417891874391412539793698933161651393444822597940152054167796341119220909092984471473051199565
     x = |var| = 133697194382348382096607788465183338193565456046018844247569458972750894312486998024824325724736982326851035901784425686447370049387786067826320738076260459479550508892238045973237342344575773985256941526172255477095635860216762897905781270700308727774042935703399597077388908867710719254521182992345010067122532664398709404997934194880983562596833367715.16055351862835494830471966680500714677379932507670091998463150017256424061653979665143123279464197235913988149831580084826877330412482794367376691844058054222030534294681422995212574445159349829995415414703858298758123239415158255798114431318612025641643715381154114556666760352704924571922668118751259862678920840770160523030157600369321997816305387348253423618868989450392980097888212485812221827024255592193265113890702116088056809490739338023291975033964255691434246678280437248028358066225404761079657847575830923796458672560700372708577233099225155755307675956899669127488545144371417891874391412539793698933161651393444822597940152054167796341119220909092984471473051199565
  x.toBigInt() = 133697194382348382096607788465183338193565456046018844247569458972750894312486998024824325724736982326851035901784425686447370049387786067826320738076260459479550508892238045973237342344575773985256941526172255477095635860216762897905781270700308727774042935703399597077388908867710719254521182992345010067122532664398709404997934194880983562596833367715
         fract = 0.16055351862835494830471966680500714677379932507670091998463150017256424061653979665143123279464197235913988149831580084826877330412482794367376691844058054222030534294681422995212574445159349829995415414703858298758123239415158255798114431318612025641643715381154114556666760352704924571922668118751259862678920840770160523030157600369321997816305387348253423618868989450392980097888212485812221827024255592193265113890702116088056809490739338023291975033964255691434246678280437248028358066225404761079657847575830923796458672560700372708577233099225155755307675956899669127488545144371417891874391412539793698933161651393444822597940152054167796341119220909092984471473051199565
gamma(191.5) = 133697194382348382096607788465183338193565456046018844247569458972750894312486998024824325724736982326851035901784425686447370049387786067826320738076260459479550508892238045973237342344575773985256941526172255477095635860216762897905781270700308727774042935703399597077388908867710719254521182992345010067122532664398709404997934194880983562596833367715.1605535186283549483047196668050071467737993250767009199846315001725642

ACTUAL -
Using BLOCK=1500 decimal places,

Big.print: var = 133697194382348382096607788465183338193565456046018844247569458972750894312486998024824325724736982326851035901784425686447370049387786067826320738076260459479550508892238045973237342344575773985256941526172255477095635860216762897905781270700308727774042935703399597077388908867710719254521182992345010067122532664398709404997934194880983562596833367715.1605535186283549483047196668050071467737993250767009199846315001725642406165397966514312368788366212835907221018006429865861935576037969689920847017491347684921760427268450103521543180431033216766352630610528625421933550878657809521731296329938415857688618490819163412729457134538944768496283846139030990974074454736971450786498287991276147257662277538418485175877317738907827654051704868340119187180819742111109203321567963552233287883215838474718307306258690368213756356472190948159600976822827872636574715641580559557972743259874324328380638560394954532125867008052074261694634579440709697433327190333921411468178548370048962295304744194143175800071297795593301187148355778549650500240776733194340373791813788411233127013175913361977517858603177367798960734841365806116683429648209693134999653555171297372610989885948877450288593568084527381016305368351571158266368749092737831257016303564518353053248881977975765003793655208366656682211306568092297891412205371815987203052261733294792617976011382077200358824303174106145381809843641541473794606620483417889940294784904376353207705112062485243357075027493962630290616547628175078072158738744008203004784710844245168307005430082
     x = |var| = 133697194382348382096607788465183338193565456046018844247569458972750894312486998024824325724736982326851035901784425686447370049387786067826320738076260459479550508892238045973237342344575773985256941526172255477095635860216762897905781270700308727774042935703399597077388908867710719254521182992345010067122532664398709404997934194880983562596833367715.1605535186283549483047196668050071467737993250767009199846315001725642406165397966514312368788366212835907221018006429865861935576037969689920847017491347684921760427268450103521543180431033216766352630610528625421933550878657809521731296329938415857688618490819163412729457134538944768496283846139030990974074454736971450786498287991276147257662277538418485175877317738907827654051704868340119187180819742111109203321567963552233287883215838474718307306258690368213756356472190948159600976822827872636574715641580559557972743259874324328380638560394954532125867008052074261694634579440709697433327190333921411468178548370048962295304744194143175800071297795593301187148355778549650500240776733194340373791813788411233127013175913361977517858603177367798960734841365806116683429648209693134999653555171297372610989885948877450288593568084527381016305368351571158266368749092737831257016303564518353053248881977975765003793655208366656682211306568092297891412205371815987203052261733294792617976011382077200358824303174106145381809843641541473794606620483417889940294784904376353207705112062485243357075027493962630290616547628175078072158738744008203004784710844245168307005430082
  x.toBigInt() = 2243066708746645393685321734362689344474497466222564090011830288188979868073765863054051074738268895685762429147992095177475799950509554621712833507984846200947666470594998420710933111780894188517836523483961583346416535484202437958951240027293550792550317525570046974480322440077898162448460863638098580418289228977672657808881821541304355522136594326164675081
         fract = 0.0615879346923771156694356347829677344174660279019808793587991475386996613410171385702634078239849988002978833763248416719337273441689255073315408118241032317788535077372003290591418417380863919615921050622886870320738231860432942652487384069542267213162067121514659251908900936793132149283985288966266097527204393215678451661259846038189607685313995224188491749206288228642798727744393741069077043870462375103713863183195945153206894896711459580831480200184641724563906945530518428061967576950044303504457375629104973235783455672111296903345675197495654609981363389136491410380345945787133575870905311647226508632570759371004183479169885995323768978516992662169844388941575653167047563260575818428625998959953554998255487144471183004257642964985874608823952319864233719102828818035019606347620276773020527061290379681941020996627964822129224854686793873261605054039178666626370514040423494345138615994703006652233886972860292406355302447935103189660567114559317129558843654388021127426989025775567719095272840401464407998026217700367878810454906880326357792552710014261423057961529376520424614212864472149744313816592732180970633470746195794327583122127331443545642984413658611712
gamma(191.5) = 2243066708746645393685321734362689344474497466222564090011830288188979868073765863054051074738268895685762429147992095177475799950509554621712833507984846200947666470594998420710933111780894188517836523483961583346416535484202437958951240027293550792550317525570046974480322440077898162448460863638098580418289228977672657808881821541304355522136594326164675081.0615879346923771156694356347829677344174660279019808793587991475386996


REPRODUCIBILITY :
This bug can be reproduced always.

---------- BEGIN SOURCE ----------
Four files are below... the build file and 3 programs. unless I email them to you, you must put each into the named files... I removed all package names; can all go in one folder.

File b_build.bat:

del *.class
javac Big.java
javac gamma.java
javac testGamma.java

-------------

File testGamma.java :

import java.math.BigDecimal;
// (c) Copyright Bruce R. 2003-infinity

public class testGamma
{
    int CALC = 70; // Calculation precision
    int MAXP = 88; // Extra precision
    int PLACES = 70; // Decimal places to print

    // -----------------------------------------------------------------------------------------

    public static void main(String args[])
    {
        new testGamma().test();
    }

    private void test()
    {
        Big.SetPrecision(MAXP,CALC);

        System.out.println("\r\n----------------------------------------------------");
        System.out.println("Test of large z integers + a fraction ");

        int x = 191; // calculating gamma(191.5)
        BigDecimal bigx = new BigDecimal(x+0.5);

        System.out.println("\r\nNext x value = " + bigx );
        System.out.println( " precision = " + Big.DEC_PLACES + ", " + Big.MAX_PLACES);

        BigDecimal val = gamma.gamma(bigx);

        System.out.println("gamma(" + Big.print(bigx,10) + ") = " + val );
        System.out.println("gamma(" + Big.print(bigx,10) + ") = " + Big.print(val,PLACES) );
    }
}

--------------

File gamma.java :

// gamma.java - Gamma functions
// (c) Copyright Bruce R. 2003-infinity

/// package bigCalc;

import java.math.BigDecimal;
import java.math.BigInteger;

public class gamma
{
    private static double ACCURACY_FACTOR = 1.0 / Math.log10(2*Math.PI);

    private static int STEP = 25; // arbitrary precision increment

    public static BigDecimal gamma(BigDecimal realzPlus1)
    {
        BigDecimal realZ = realzPlus1.subtract(Big.ONE);

        // Get internal precision.
        Big origPrec = Big.GetPrecision();

        int precision = origPrec.dec;

        BigDecimal sum = null;
        BigDecimal prev = null;

        int stepFactor = ((int)realzPlus1.doubleValue()) / 10;
        if (stepFactor < 1)
            stepFactor = 1;

        for (boolean again = true; again; precision += STEP * stepFactor, stepFactor = 1)
        {
            int a = (int)Math.ceil(precision * ACCURACY_FACTOR);
            BigDecimal bigA = new BigDecimal(a);

            System.out.println("Gamma. a = " + a + ". Setting precision = " + (precision+STEP) +
                                ", " + (precision+STEP+10));

           Big.SetPrecisionAndExponent(precision+STEP+10,precision+STEP,Big.NO_MAX_POWER);

            BigDecimal realZplusA = realZ.add(bigA);
            BigDecimal realZplusHalf = realZ.add(Big.HALF);

            BigDecimal lnBase = Big.ln(realZplusA); // termNumer = (z+a)^(z+1/2)
            lnBase = Big.mult(lnBase,realZplusHalf);
            BigDecimal termNumer = Big.exp(lnBase);
            if (termNumer == null)
            {
                System.out.println("****** UNEXPECTED ERROR *******");
                continue;
            }

            BigDecimal bigTermDenom = Big.exp(realZ);
            BigDecimal termDenom = bigTermDenom.negate();

            // Calculate the series... series is valid for real(z) >= -1/2.

            // Term k=0
            BigDecimal sqrt2pi = Big.sqrt(Big.TWO_PI); // sqrt(2pi)

            BigDecimal expZA = Big.exp( realZ.add(bigA) ); // e^(z+a)

            BigDecimal temp = Big.mult(termNumer,sqrt2pi); // sum = (z+a)^(z+1/2)*sqrt(2pi)/e^(z+a)
            temp = Big.div(temp,expZA);
            sum = temp;

            // Terms k=1 thru k = a-1
            BigDecimal fact = Big.ONE; // factorial number
            BigDecimal zPlusK = realZ; // will be each z+k
            for (int k = 1; k <= a-1; k++)
            {
                BigInteger iTerm = new BigInteger(""+(a-k)).pow(k-1);
                BigDecimal term = new BigDecimal(iTerm);
                term = Big.mult(term,Big.sqrt( new BigDecimal(a-k) ));

                if (k >= 3)
                    fact = fact.add(Big.ONE);

                termDenom = Big.mult(termDenom,fact);
                termDenom = Big.mult(termDenom,Big.E);
                termDenom = termDenom.negate();

                BigDecimal cTerm = Big.mult(term,termNumer);
                cTerm = Big.div(cTerm,termDenom);

                zPlusK = zPlusK.add(Big.ONE);
                cTerm = Big.div(cTerm,zPlusK);
                sum = sum.add(cTerm);
            }

            if (prev != null)
            {
                System.out.println(" Prev sum cmp = " + prev);
                System.out.println(" Next sum cmp = " + sum);
                System.out.println(" CONSTANTS precision: DEC = " + Big.CONSTANT_DEC_PLACES +
                                    ", CONSTANT_MAX = " + Big.CONSTANT_MAX_PLACES);

                again = Big.cmp(prev,sum,origPrec.dec) != 0;
            }
            prev = sum;
            System.out.println(" Loop again = " + again + ", places compared = " + origPrec.dec);
        }

        Big.ResetPrecision(origPrec);
        return sum;
    }
}

-----------

File Big.java :


// Big.java - Implementation of BigDecimal functions.
// (c) Copyright Bruce R. 2003-infinity

/// package bigCalc;

import java.math.BigDecimal;
import java.math.BigInteger;

public class Big
{
    public static int ROUND_DOWN = BigDecimal.ROUND_DOWN;

    // -------------------------

    // Numeric constants

    public static BigDecimal ZERO = new BigDecimal("0");
    public static BigDecimal ONE = new BigDecimal("1");
    public static BigDecimal TWO = new BigDecimal("2");
    public static BigDecimal FOUR = new BigDecimal("4");
    public static BigDecimal SEVEN = new BigDecimal("7");
    public static BigDecimal TEN = new BigDecimal("10");
    public static BigDecimal ONE_HUND = new BigDecimal("100");
    public static BigDecimal NEG_ONE = new BigDecimal("-1");
    public static BigDecimal HALF = new BigDecimal("0.5");
    public static BigDecimal NEG_HALF = new BigDecimal("-0.5");


    // These do NOT need to be re-calculated when the precision is changed:

    private static BigDecimal ONE_TENTH = new BigDecimal("0.1"); // Minimum arg for ln(x)



    // These must be re-calculated when the precision is changed:

    public static BigDecimal LN_10 = ln(TEN);
    public static BigDecimal E = taylor(ONE,0,1,false);
    public static BigDecimal PI = ComputePi();

    public static BigDecimal TWO_PI = mult(TWO,PI);
    public static BigDecimal PI_HALF = div(PI,TWO);
    public static BigDecimal PI_QTR = div(PI,FOUR);


    // -------------------------

    // Precision parameters

    // This default maximum is currently experimental:
    public static final int NO_MAX_POWER = -1;

    public static int MAX_EXPONENT = NO_MAX_POWER;

    public static int DEC_PLACES; // # decimal places for calculations.
    public static int MAX_PLACES; // extra decimal places for some calculations.
    private static int SLACK = 15; // Some extra digits after DEC_PLACES.


    private static final int BLOCK = 1000; // xxxxxxxx THIS WORKS
// private static final int BLOCK = 1500; // THIS CAUSES THE ERROR
    public static int CONSTANT_DEC_PLACES = 0;
    public static int CONSTANT_MAX_PLACES = 0+SLACK;


    private static int saved_DEC_PLACES = -111;
    private static int saved_MAX_PLACES = -111;


    // Set the default precision
    private static int dummy = SetPrecisionAndExponent(28,20,NO_MAX_POWER);


    // --------------------
    // For use as the precision object:
    public int dec; // DEC_PLACES
    public int max; // MAX_PLACES
    public int exp; // MAX_EXPONENT


    // =============================================================================================
    // PRECISION FUNCTIONS

    public static int SetPrecisionAndExponent(int maxPlaces, int decPlaces, int maxExponent)
    {
        SetPrecision(maxPlaces,decPlaces);
        MAX_EXPONENT = maxExponent; // Set the max exponent for exp(), pow() and similar funcs.
        return 0; // need return value for static call at top
    }

    // -------------------------

    public static int SetPrecision(int maxPlaces, int decPlaces)
    {
        if (decPlaces < 1)
            decPlaces = 1;
        if (maxPlaces < decPlaces)
            maxPlaces = decPlaces + SLACK;

        if (decPlaces > CONSTANT_DEC_PLACES || maxPlaces > CONSTANT_MAX_PLACES)
        {
System.out.println(">>>> NEED TO RECALC. CONSTANT_DEC = " + CONSTANT_DEC_PLACES + ", CONSTANT_MAX = " + CONSTANT_MAX_PLACES);
            // Advance the constant precision to the next 'BLOCK' precision size.
            // For 1 to 500, const prec = 500; for 501 to 1000, const prec = 1000; etc.

            int prec = Math.max(decPlaces,maxPlaces);
            int size = ((prec + BLOCK-1)/BLOCK) * BLOCK; // size = 500*k that is >= precision
            CONSTANT_DEC_PLACES = size;
            CONSTANT_MAX_PLACES = size+SLACK;
System.out.println(" New constant sizes are: DEC = " + CONSTANT_DEC_PLACES + ", CONSTANT_MAX = " + CONSTANT_MAX_PLACES);

            DEC_PLACES = size;
            MAX_PLACES = size+SLACK;

            CalcConstants(); // re-calculate the mathematical constants.
        }
        DEC_PLACES = decPlaces;
        MAX_PLACES = maxPlaces;
        return 0; // need return value for static call at top
    }

    // -------------------------

    public static Big GetPrecision()
    {
        Big out = new Big();
        out.dec = DEC_PLACES;
        out.max = MAX_PLACES;
        out.exp = MAX_EXPONENT;
        return out;
    }

    public static void ResetPrecision(Big orig)
    {
        DEC_PLACES = orig.dec;
        MAX_PLACES = orig.max;
        MAX_EXPONENT = orig.exp;
    }

    // =============================================================================================
    // MATHEMATICAL FUNCTIONS

    // recip() computes the reciprocal of a BigDecimal number.

    public static BigDecimal recip(BigDecimal x)
    {
        return div(ONE,x);
    }

    // -------------------------

    public static BigDecimal mult(BigDecimal x1, BigDecimal x2)
    {
        BigDecimal prod = x1.multiply(x2);
        return prod.setScale(MAX_PLACES,ROUND_DOWN);
    }

    // -------------------------

    public static BigDecimal div(BigDecimal x1, BigDecimal x2)
    {
        if (x2.compareTo(ZERO) == 0)
            return null;

        BigDecimal dividend = x1.setScale(MAX_PLACES,ROUND_DOWN);
        return dividend.divide(x2,ROUND_DOWN);
    }

    // -------------------------

    public static BigDecimal[] mod(BigDecimal x1, BigDecimal x2)
    {
        BigDecimal quot,prod;
        BigDecimal qr[] = new BigDecimal[2];

        if (x2.compareTo(ZERO) == 0)
            return null;

        BigDecimal x2Abs = x2.abs();
        quot = div(x1,x2Abs);
        qr[0] = new BigDecimal(quot.toBigInteger());
        prod = mult(x2Abs,qr[0]);
        qr[1] = x1.subtract(prod);
        return qr;
    }

    // =============================================================================================

    public static int cmp(BigDecimal x, BigDecimal y)
    {
        return cmp(x,y,DEC_PLACES);
    }

    // -------------------------

    public static int cmp(BigDecimal x, BigDecimal y, int places)
    {
        BigInteger ix = x.movePointRight(places).toBigInteger();
        BigInteger iy = y.movePointRight(places).toBigInteger();
        return ix.compareTo(iy);
    }

    // =============================================================================================

    public static String print(BigDecimal var, int decPlaces)
    {
System.out.println("\r\nBig.print: var = " + var);
        boolean negative = var.signum() < 0;

        BigDecimal x = var.abs();
System.out.println(" x = |var| = " + x);

        BigInteger whole = x.toBigInteger();
System.out.println(" x.toBigInt() = " + whole);
        BigDecimal fract = x.subtract( new BigDecimal(whole) );
System.out.println(" fract = " + fract);

        // Format the integer part.
        String sWhole = "";
        while (whole.compareTo(BigInteger.ONE) >= 0)
        {
            BigInteger qr[] = whole.divideAndRemainder(BigInteger.TEN);
            sWhole = qr[1] + sWhole;
            whole = qr[0];
        }
        if (sWhole.equals(""))
            sWhole = "0";

        // Format the fraction part.
        int trailZeroes = 0;
        String sFract = ".";
        for (int i = Math.abs(decPlaces); i > 0; i--)
        {
            fract = fract.movePointRight(1);
            whole = fract.toBigInteger();
            sFract += whole;
            if (whole.compareTo(BigInteger.ZERO) == 0)
                trailZeroes++;
            else
                trailZeroes = 0;
            fract = fract.subtract( new BigDecimal(whole) );
        }

        // Remove trailing zeros from the fraction string if not required to keep them.
        if (decPlaces >= 0 && trailZeroes > 0)
        {
            int len = sFract.length() - trailZeroes;
            sFract = sFract.substring(0,len);
        }
        if (sFract.length() == 1)
            sFract = "";

        String sOut = negative ? "-" : "";
        sOut += sWhole + sFract;

        return sOut;
    }


    // ======================================================================

    public static void CalcConstants()
    {
        LN_10 = ln(TEN);
        E = taylor(ONE,0,1,false);
        PI = ComputePi();

        TWO_PI = mult(TWO,PI);
        PI_HALF = div(PI,TWO);
        PI_QTR = div(PI,FOUR);
    }

    // -------------------------

    // pi = 16*atan(1/5) - 4*atan(1/239) )

    private static BigDecimal ComputePi()
    {
        BigDecimal one5th = new BigDecimal("0.2");
        BigDecimal atan5 = mult(atan(one5th),FOUR);
        BigDecimal two39 = new BigDecimal("239");
        BigDecimal atan239 = atan(recip(two39));
        atan5 = atan5.subtract(atan239);
        atan5 = mult(atan5,FOUR);

        return atan5.setScale(MAX_PLACES,ROUND_DOWN);
    }

    // =================================================================================
    // Taylor series
    //
    // Input: startTerm = first term's value of series
    // startIndex = index (i.e. exponent or divisor) of first term
    // delta = increment of exponent and factorial divisor from term to term
    // alternating = true if alternating series, false if positive series

    private static BigDecimal taylor(BigDecimal x, int first, int delta, boolean alternating)
    {
        x = x.setScale(MAX_PLACES,ROUND_DOWN);
        BigDecimal sum = ZERO;
        int sign = 1; // first term is positive

        // Generate first term = x^first
        BigDecimal xToN = ONE;
        BigDecimal nFact = ONE;
        BigDecimal n = ZERO;
        for (int i = 1; i <= first; i++)
        {
            xToN = mult(xToN,x);
            n = n.add(ONE);
            nFact = mult(nFact,n); // integer
        }

        BigDecimal term;
        do
        {
            // Calculate and add current term.
            term = div(xToN,nFact);
            if (sign > 0)
                sum = sum.add(term);
            else
                sum = sum.subtract(term);

            // Prepare for next term.
            if (alternating)
                sign = -sign;

            for (int i = 1; i <= delta; i++)
            {
                xToN = mult(xToN,x);
                n = n.add(ONE);
                nFact = mult(nFact,n); // integer
            }
        }
        while (cmp(term,ZERO) != 0);

        return sum.setScale(MAX_PLACES,ROUND_DOWN);
    }

    // =================================================================================

    // exp(x)
    //
    // Output: e^x, OR null if x is too large

    public static BigDecimal exp(BigDecimal x)
    {
        BigDecimal lowLimit = mult(LN_10,new BigDecimal(-MAX_PLACES));
        if (cmp(x,lowLimit) < 0)
            return ZERO;

        // If x is too large, return error.
        if (MAX_EXPONENT > 0 && x.doubleValue() > MAX_EXPONENT)
            return null;

        // If |x| <= 100, do Taylor series normally... converges quickly for arg <= 100.
        if (cmp(x.abs(),ONE_HUND) <= 0)
            return taylor(x,0,1,false);

        // For |x| > 100 (large x), convert +/-x to a modulo-100 integer and remainder, then
        // compute e^x = e^(x%100) * (e^+/-100)^[x/100]. Although virtually impossible for
        // x/100 to be beyond integer precision, treat it as +/-infinity if it is.

        BigDecimal qr[] = mod(x,ONE_HUND); // quot & rmdr both have same sign as x

        Integer count = qr[0].abs().intValue();

        BigDecimal ex = taylor(qr[1],0,1,false);

        BigDecimal arg = x.signum() >= 0 ? ONE_HUND : ONE_HUND.negate();
        BigDecimal e100 = taylor(arg,0,1,false);
        for (int i = 1; i <= count; i++)
        {
            ex = mult(ex,e100);
            if (ex.compareTo(ZERO) == 0) // If x < 0 and |x| is so large that the
                return ZERO; // result becomes zero, return the zero result.
        }
        return ex;
    }

    // ======================================================================
    // ln(x)
    // valid only for x > 0. ln(x) range is all x.
    //
    // Output: ln x if x > 0, OR
    // null if x <= 0

    public static BigDecimal ln(BigDecimal x)
    {
        if (cmp(x,ZERO) <= 0)
            return null;

        if (cmp(x,ONE) == 0)
            return ZERO;

        int count = 0;
        while (cmp(x,ONE_HUND) > 0)
        {
            x = x.movePointLeft(1);
            count++;
        }
        while (cmp(x,ONE_TENTH) < 0)
        {
            x = x.movePointRight(1);
            count--;
        }

        BigDecimal numer = x.subtract(ONE);
        BigDecimal denom = x.add(ONE);
        x = div(numer,denom);

        BigDecimal xterm = x;

        BigDecimal term = xterm;
        BigDecimal xSquared = mult(x,x);

        BigDecimal logx = ZERO;
        BigDecimal div = ONE;

        do
        {
            logx = logx.add(term);
            xterm = mult(xterm,xSquared);
            div = div.add(TWO);
            term = xterm.divide(div,ROUND_DOWN);
        }
        while (cmp(term,ZERO) != 0);

        logx = mult(logx,TWO);
        if (count != 0)
        {
            int extra = 1;
            int c = Math.abs(count);
            if (c > 0)
                extra = String.valueOf(c).length() + 1;

            Big origPrec = GetPrecision();
            SetPrecision(origPrec.max+extra,origPrec.dec+extra);

            BigDecimal prod = mult(LN_10, new BigDecimal(count));
            logx = logx.add(prod);
            ResetPrecision(origPrec);
        }
        return logx.setScale(MAX_PLACES,ROUND_DOWN);
    }

    // ======================================================================
    // sqrt(x)
    // valid only for x >= 0. sqrt(x) range is >= 0.
    //
    public static BigDecimal sqrt(BigDecimal x)
    {
        int n = 2; // n=2 --> square root

        if (cmp(x,ZERO) == 0)
            return ZERO;
        if (cmp(x,ONE) == 0)
            return ONE;

        // The process may need two tries if convergence fails. On second try, increase
        // the internal MAX precision but keep the DEC precision as is.

        Big origPrec = GetPrecision();

        BigDecimal root = null;
        int count = 0;
        for (int tries = 2; tries > 0; tries--)
        {
            // Set up constants.
            BigDecimal bigN = new BigDecimal(n);
            BigDecimal factor = div( bigN.subtract(ONE), bigN);
            BigDecimal numDivN = div(x,bigN);

            // starting point for the calculation.
            root = ONE;
            if (cmp(x,ONE) > 0)
            {
                int digits = x.toBigInteger().toString().length();
                digits = (digits + n-1) / n;
                if (digits > 0)
                    root = new BigDecimal("5").movePointRight(digits-1);
            }

            // The iteration formula is root' = (n-1)/n * root + (x/n / root^(n-1))
            count = 1000; // Iteration limit
            BigDecimal prior;
            do
            {
                prior = root;

                BigDecimal a = numDivN;
                for (int i = 1; i < n; i++)
                    a = a.divide(root,ROUND_DOWN);
                root = mult(root,factor).add(a);
            }
            while (cmp(root,prior) != 0 && --count > 0);

            if (count > 0) // If didn't get stuck,
                break; // the root was found.

            SetPrecision(origPrec.max+origPrec.dec/2,origPrec.dec);
        }

        ResetPrecision(origPrec);
        root = root.setScale(MAX_PLACES,ROUND_DOWN);

        return root;
    }

    // -------------------------

    // atan(x)
    // valid for all x. atan(x) range is -pi/2 <= atan(x) <= +pi/2.

    public static BigDecimal atan(BigDecimal x)
    {
        BigDecimal atanx;

        // If the result of tan(x) was null, then x was 90 degrees.
        if (x == null)
            return PI_HALF;

        if (cmp(x,ONE) == 0)
            return PI_QTR;

        if (cmp(x,NEG_ONE) == 0)
            return PI_QTR.negate();

        BigDecimal absx = x.abs();
        if (cmp(absx,ONE) <= 0)
            atanx = atanSeries(absx);
        else
        {
            BigDecimal arg = recip(absx);
            atanx = atanSeries(arg);
            atanx = PI_HALF.subtract(atanx);
        }

        if (cmp(x,ZERO) < 0)
            atanx = atanx.negate();

        return atanx;
    }

    // -------------------------

    // Euler's Arctan series (converges faster than the logarithm-type series)

    private static BigDecimal atanSeries(BigDecimal x)
    {
        x = x.setScale(MAX_PLACES,ROUND_DOWN);

        BigDecimal xsq = mult(x,x);
        BigDecimal div = xsq.add(ONE);

        BigDecimal term = x.divide(div,ROUND_DOWN);
        BigDecimal xfactor = xsq.divide(div,ROUND_DOWN);

        BigDecimal sum = ZERO;
        BigDecimal n = ZERO;
        div = ONE;

        do
        {
            sum = sum.add(term);

            n = n.add(ONE);
            BigDecimal numer = mult(mult(n,n),FOUR); // numer is an integer
            term = mult(term,numer);
            term = mult(term,xfactor);

            div = div.add(ONE); // div is an integer
            term = term.divide(div,ROUND_DOWN);

            div = div.add(ONE);
            term = term.divide(div,ROUND_DOWN);
        }
        while (cmp(term,ZERO) != 0);

        sum = sum.setScale(MAX_PLACES,ROUND_DOWN);
        return sum;
    }
}
---------- END SOURCE ----------

CUSTOMER SUBMITTED WORKAROUND :
have not tried yet but probably can find one