/*  RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
 *  CONDITION ABS(X) < 1/B, ERROR OTHERWISE.
 *  USES NEWTONS METHOD TO SOLVE THE EQUATION
 *  EXP1(-Y) = X, THEN REVERSES SIGN OF Y.
 */
static void
mp_lns(const MPNumber *x, MPNumber *z)
{
    int t, it0;
    MPNumber t1, t2, t3;

    /* ln(1+0) = 0 */
    if (mp_is_zero(x)) {
        mp_set_from_integer(0, z);
        return;
    }

    /* Get starting approximation -ln(1+x) ~= -x + x^2/2 - x^3/3 + x^4/4 */
    mp_set_from_mp(x, &t2);
    mp_divide_integer(x, 4, &t1);
    mp_add_fraction(&t1, -1, 3, &t1);
    mp_multiply(x, &t1, &t1);
    mp_add_fraction(&t1, 1, 2, &t1);
    mp_multiply(x, &t1, &t1);
    mp_add_integer(&t1, -1, &t1);
    mp_multiply(x, &t1, z);

    /* Solve using Newtons method */
    it0 = t = 5;
    while(1)
    {
        int ts2, ts3;

        /* t3 = (e^t3 - 1) */
        /* z = z - (t2 + t3 + (t2 * t3)) */
        mp_epowy(z, &t3);
        mp_add_integer(&t3, -1, &t3);
        mp_multiply(&t2, &t3, &t1);
        mp_add(&t3, &t1, &t3);
        mp_add(&t2, &t3, &t3);
        mp_subtract(z, &t3, z);
        if (t >= MP_T)
            break;

        /*  FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
         *  BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,
         *  WE CAN ALMOST DOUBLE T EACH TIME.
         */
        ts3 = t;
        t = MP_T;
        do {
            ts2 = t;
            t = (t + it0) / 2;
        } while (t > ts3);
        t = ts2;
    }

    /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
    if (t3.sign != 0 && t3.exponent << 1 > it0 - MP_T) {
        mperr("*** ERROR OCCURRED IN MP_LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
    }

    z->sign = -z->sign;
}