Memory issue

[Craig Helfgott]

Okay, here is my exp subroutine, can you please tell me what I am doing
wrong?


mpf_class mpf_exp(mpf_class b){    // a = exp(b)
    mpf_class a;
    int numsquares, inv;
    mpf_class temp_exp;

    temp_exp = b;
    inv = 1;
    if (sgn(b) < 0) {
        temp_exp = -b;
        inv = -1;
    }

    numsquares = 0;
    while(cmp(temp_exp, 1) >= 0){
        mpf_div_2exp(temp_exp.get_mpf_t(), temp_exp.get_mpf_t(), 1);
        numsquares++;
    }

    unsigned long long int in_prec;
    in_prec = temp_exp.get_prec();

    // Here we use binary splitting.
    // e^{p/q} = \sum_{k=0}^{NTerms -1} (p/q)^k /k!
    //         = P(0, NTerms) / Q(0, NTerms)   if p^2 <= q <= NTerms.

    // Q(a,b) = b! q^{b-a-1} / a!
    // P(a,b) = b! q^{b-a-1} / a! + b! q^{b-a-2} p / (a+1)! + ... + b
p^{b-a-1}

    // m := (a+b)/2,  Q(a,b) = q * Q(a,m) * Q(m,b),   P(a,b) = q * P(a,m) *
Q(m,b) + p^{m-a} * P(m,b)
    // Q(a,a+1) = P(a,a+1) = a+1

    mpz_class numer, FinalP, FinalQ;
    mpf_class temp_numer;

    unsigned long long int lgNTerms, NTerms, lglgdenom, lgdenom;
    int i,j,k;

    lglgdenom = 0; lgdenom = 1;
    FinalP = 1;
    FinalQ = 1;

    // We want the number of terms to be 2^n, with (2^n)! q^{2^n} >
2^{PrecisionForExp}
    // So  n (2^n-1/2) - 1.45 2^n + 2^n lgdenom > PrecisionForExp
    i = lgdenom; lgNTerms = 0;
    while(i < PrecisionForExp) {
        lgNTerms++;
        i = ((lgNTerms -2 + lgdenom) << lgNTerms) - (lgNTerms+1)/2;
    }
    NTerms = 1 << lgNTerms;

    mpz_class P_values[lgNTerms+2], Q_values[lgNTerms+2];
    mpz_class tmp_P, tmp_Q, tmpvar1;

    while((lgdenom < in_prec) && (lgdenom < PrecisionForExp)) {
        if (lgdenom == 1) {
            mpf_mul_2exp(temp_exp.get_mpf_t(), temp_exp.get_mpf_t(), 1);
        } else {
            mpf_mul_2exp(temp_exp.get_mpf_t(), temp_exp.get_mpf_t(),
lgdenom/2);
        }
        temp_numer = floor(temp_exp);
        temp_exp = temp_exp - temp_numer;
        numer = temp_numer;
        if (sgn(temp_numer) > 0) {
            i = PrecisionForExp + 1;
            while(i > PrecisionForExp) {
                lgNTerms--;
                i = ((lgNTerms -2 + lgdenom) << lgNTerms) - (lgNTerms+1)/2;
            }
            lgNTerms++;
            NTerms = 1 << lgNTerms;
            // We run this as a depth-first traversal of a binary tree.
Each node (a,b)
            // has two children, (a,m) and (m,b) in that order.  We traverse
left child
            // then right child, then parent.
            for(i=1; i<= NTerms; i++) {
                tmp_P = i;
                tmp_Q = i;
                numer = temp_numer;
                k = i;
                j = 0;
    // m := (a+b)/2,  Q(a,b) = q * Q(a,m) * Q(m,b),   P(a,b) = q * P(a,m) *
Q(m,b) + p^{m-a} * P(m,b)
    // Q(a,a+1) = P(a,a+1) = a+1
                while((k%2) == 0) {
                    tmpvar1 = P_values[j] * tmp_Q;
                    mpz_mul_2exp(tmpvar1.get_mpz_t(), tmpvar1.get_mpz_t(),
lgdenom);
                    tmp_P = tmpvar1 + numer * tmp_P;
                    tmp_Q = Q_values[j] * tmp_Q;
                    mpz_mul_2exp(tmp_Q.get_mpz_t(), tmp_Q.get_mpz_t(),
lgdenom);
                    numer = numer * numer;
                    k /= 2;
                    j++;
                }
                P_values[j] = tmp_P;
                Q_values[j] = tmp_Q;
            }
            FinalP = FinalP * P_values[j];
            FinalQ = FinalQ * Q_values[j];
        }
        lglgdenom++;
        lgdenom *= 2;
    }

    mpq_class ratio;
    if (inv > 0) {
        ratio = FinalP;
        ratio = ratio / FinalQ;
    } else {
        ratio = FinalQ;
        ratio = ratio / FinalP;
    }
    a = ratio;

    if (numsquares > 0) {
        mpf_pow_ui(a.get_mpf_t(), a.get_mpf_t(), 1 << numsquares);
    }

    return a;
}



On Fri, Dec 17, 2010 at 3:39 AM, Torbjorn Granlund <tg at gmplib.org> wrote:

> Craig Helfgott <chelfgott at gmail.com> writes:
>
>   1) Do I need to run
>  mpf_clear if I am using mpf_class instead of mpf_t?
>
> You need to match new's and delete's in C++.
>
> Vague questions cannot reslt in anythiung but at best vague replies.  If
> you want god replies, state the facts carefully, with example code.
>
>  2)  Is there any reason
>  that g++ with the -lgmp and -lgmpxx flags would automatically parallelize
> a
>  loop at compile-time?
>
> See the linker documentation (man ld) for the function of the -l option.
>
> --
> Torbjörn
>