# Memory issue

Craig Helfgott chelfgott at gmail.com
Sat Dec 18 01:59:12 CET 2010

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
>