Added metaMPFR in the tools directory of MPFR.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@7274 280ebfd0-de03-0410-8827-d642c229c3f4

3 files changed, 1467 insertions, 0 deletions

diff --git a/tools/metaMPFR/metaMPFR_common.mpl b/tools/metaMPFR/metaMPFR_common.mpl
new file mode 100644
index 000000000..7917be78b
--- /dev/null
+++ b/tools/metaMPFR/metaMPFR_common.mpl
@@ -0,0 +1,521 @@
+#
+# These procedures come from the source code of GFun.
+#
+getname:=proc(yofz::function(name), y, z)
+  y:=op(0,yofz);
+  if type(y,'procedure') then error `not an unassigned name`,y fi;
+  z:=op(yofz)
+end proc:
+
+
+#
+# returns the smallest i such that u(n+i) appears in a recurrence
+#
+minindex := proc(rec,u,n)
+  min(op(map(op,indets(rec,'specfunc'('linear'(n),u)))))-n
+end proc:
+
+
+#
+# returns the largest i such that u(n+i) appears in a recurrence
+#
+maxindex := proc(rec,u,n)
+  max(op(map(op,indets(rec,'specfunc'('linear'(n),u)))))-n
+end proc:
+
+
+#
+# A recurrence of the form a(n+d) = p(n)/q(n) a(n) is represented through a record:
+# OneTermRecurrence : record(order, numerator, denominator)
+#
+`type/OneTermRecurrence` := 'record(order, numerator, denominator)':
+
+
+#
+#checkOneTermRecurrence
+# Input: a recurrence rec (either with or without initial conditions).
+#          If it has initial conditions, they are ignored.
+#        a(n): the name of the sequence and the name of the variable.
+#
+# Output:
+# This procedure checks that rec is a recurrence of the form a(n+d) = p(n)/q(n) a(n)
+# If the check succeeds, it returns the corresponding record. If it fails, an error is
+# returned.
+#
+checkOneTermRecurrence := proc(rec, aofn)::OneTermRecurrence;
+  local r, d, a, n, term1, term2, res;
+
+  getname(aofn, a, n):
+  if type(rec, 'set') then
+    r:=select(has, rec, n);
+    if nops(r)>1
+      then error `invalid recurrence`, rec
+    fi:
+    if nops(r)=0
+      then error "%1 does not appear in the recurrence", n
+    fi:
+    r := op(r):
+  else r:=rec:
+  fi:
+  if type(r,'`=`')
+    then r:=op(1,r)-op(2,r)
+  fi:
+  if indets(r,'specfunc'('anything',a)) <> indets(r,'specfunc'('linear'(n),a))
+    then error "the recurrence contains elements that are not linear in %1", n
+  fi:
+  if nops(r) <> 2
+    then error "the recurrence contains %1 terms (expected 2)", nops(r)
+  fi:
+  r := subs(n=n-minindex(r, a, n), r):
+  d := maxindex(r, a, n):
+
+  term1 := select(has, r, a(n)):
+  term2 := select(has, r, a(n+d)):
+
+  res := factor( -(term1/a(n)) / (term2/a(n+d)) ):
+
+  Record( 'order'=d, 'numerator' = numer(res), 'denominator' = denom(res) )
+end proc:
+
+
+#
+# my_factors factorizes p the same way as factors(p) would do except that the constant part is computed
+# differently. We assume here that p has integer coefficients, and we want to factorize it over polynomials
+# with integer coefficients. my_factors ensures that the factors have integer coefficients.
+#
+my_factors := proc(p)
+  local L, c, fact, i, my_c, my_fact, q:
+  L := factors(p):
+  c := L[1]: fact := L[2]:
+  my_c := c: my_fact := []:
+  for i from 1 to nops(fact) do
+    q := denom(fact[i][1]):
+    my_fact := [ op(my_fact), [ fact[i][1]*q, fact[i][2] ] ]:
+    my_c := my_c / (q^fact[i][2]):
+  od:
+  [ my_c, my_fact]:
+end proc:
+
+
+#
+# This procedure decomposes a one-term recurrence with the following form:
+# a(n+d) = c * s1(n)/s1(n+d) * s2(n+d)/s2(n) * p(n)/q(n) * a(n)
+#
+# Known issue: this procedure assumes that the only variables involved are n and x with their usual meaning.
+#
+decomposeOneTermRecurrence := proc(formalRec::OneTermRecurrence, res_cste, res_s1, res_s2, res_p, res_q)
+  local p, q, cste, s1, s2, d, L, i, tmp, exponent, r, polyring;
+  p := formalRec:-numerator:
+  q := formalRec:-denominator:
+  d := formalRec:-order:
+  s1 := 1:
+  L := op(2,my_factors(p)): # L contains the non trivial factors of p
+  for i from 1 to nops(L) do
+    tmp := L[i][1]: exponent := L[i][2]:
+    r := gcd(tmp^exponent, subs(n=n-d, q)):
+    p := quo(p,r,n): q := quo(q, subs(n=n+d, r),n): s1 := s1 * r:
+  od:
+
+  s2 := 1:
+  L := op(2,my_factors(p)): # L contains the *remaining* non trivial factors of p
+  for i from 1 to nops(L) do
+    tmp := L[i][1]: exponent := L[i][2]:
+    r := gcd(tmp^exponent, subs(n=n+d, q)):
+    p := quo(p, r, n): q := quo(q, subs(n=n-d, r), n): s2 := s2 * r:
+  od:
+
+  # Finally we look for the constant part (with respect to n) of p/q
+  cste := op(1, my_factors(p))/op(1, my_factors(q)):
+  p := p/op(1, my_factors(p)): q := q/op(1, my_factors(q)):
+  polyring := RegularChains[PolynomialRing]([n,x]):
+  L := op(2, my_factors(p)):
+  for i from 1 to nops(L) do
+    if RegularChains[MainVariable](L[i][1], polyring) = x
+      then cste := cste * L[i][1]^L[i][2]: p := quo(p,L[i][1]^L[i][2],x):
+    fi:
+  od:
+  L := op(2, my_factors(q)):
+  for i from 1 to nops(L) do
+    if RegularChains[MainVariable](L[i][1], polyring) = x
+      then cste := cste / L[i][1]^L[i][2]: q := quo(q,L[i][1]^L[i][2],x):
+    fi:
+  od:
+
+  res_cste := cste;
+  res_s1 := s1;
+  res_s2 := s2;
+  res_p := simplify(p);
+  res_q := simplify(q);
+end proc:
+
+
+#
+#coeffrecToTermsrec
+# Input: a linear recurrence rec (either with or without initial conditions).
+#           a(n): the name of the sequence and the name of the variable.
+#           x: a value or symbolic name
+#
+# Output:
+# The recurrence satisfied by a(n)*x^n. Note that this recurrence is also denoted by a(n).
+# If initial conditions were provided, corresponding initial conditions are computed.
+#
+coeffrecToTermsrec := proc(rec, aofn, x)
+  local a,n,L,r,cond,d,i,tmp,c,res;
+  getname(aofn, a, n):
+  if type(rec, 'set') then
+    L := selectremove(has, rec, n):
+    r := L[1]:
+    if nops(r)>1
+      then error `invalid recurrence`, rec
+    fi:
+    if nops(r)=0
+      then error "%1 does not appear in the recurrence", n
+    fi:
+    r := op(r):
+    cond := L[2]:
+  else r := rec:
+  fi:
+  d := maxindex(r, a, n):
+  L := indets(r,'specfunc'('linear'(n),a)):
+  if indets(r,'specfunc'('anything',a)) <> L
+    then error "the recurrence contains elements that are not linear in %1", n
+  fi:
+  L := map(op, L):
+  for i from 1 to nops(L) do
+    r := subs(a(op(i,L))=a(op(i,L))*x^(d-op(i,L)+n), r):
+  od:
+  if cond<>'cond' then
+    c := {}:
+    for i from 1 to nops(cond) do
+      tmp := op(i, cond): # tmp should have the form 'a(k) = cste'
+      if not type(tmp,'`=`') then error "Invalid initial condition: %1", tmp: fi:
+      L := selectremove(has, {op(tmp)}, a):
+      if (nops(L[1]) <> 1) or (nops(L[2])<>1)
+        then error "Invalid initial condition: %1", tmp:
+      fi:
+      tmp := op(1, L[1]): # tmp has the form 'a(k)'
+      c := {op(c), tmp = op(1, L[2])*x^op(tmp)}:
+    od:
+    res := {r, op(c)}:
+  else res := r:
+  fi:
+  res:
+end proc:
+
+
+#
+# This procedure removes the conditions of the form a(k)=0 from the initial conditions of rec
+# It returns a list L = [L1, L2, ...] where Li = [k, expr] representing the condition a(k)=expr.
+# Moreover, it asserts that the Li are ordered by increasing k.
+#
+removeTrivialConditions := proc(rec, aofn)
+  local a,n,i,L,tmp,c,cond,k:
+  getname(aofn, a, n):
+  if not type(rec, 'set') then
+    error "%1 is not a recurrence with initial conditions", rec
+  else
+    L := selectremove(has, rec, n):
+    cond := L[2]:
+    if nops(cond)=0
+      then error "%1 does not contain initial conditions", rec
+    fi:
+  fi:  
+  c := []:
+  for i from 1 to nops(cond) do
+    tmp := op(i, cond): # tmp should have the form 'a(k) = cste'
+    if not type(tmp,'`=`') then error "Invalid initial condition: %1", tmp: fi:
+    L := selectremove(has, {op(tmp)}, a):
+    if (nops(L[1]) <> 1) or (nops(L[2])<>1)
+      then error "Invalid initial condition: %1", tmp:
+    fi:
+    if op(1, L[2])<>0 then c := [op(c), [op(op(1, L[1])), op(1, L[2])]]: fi:
+  od:
+  # We check that the conditions are ordered by increasing k.
+  if (nops(c)=0) then return c: fi:
+  k := c[1][1]:
+  for i from 2 to nops(c) do
+    if (c[i][1]<=k)
+    then error "Unexpected error in removeTrivialConditions: the conditions are not correctly ordered (%1)\n", c
+    else k := c[i][1]
+    fi:
+  od:
+  c:
+end proc:
+
+
+#
+# findFixpointOfDifferences: takes a set L of integer and returns the smallest set S
+# containing L and such that for each i, S[i]-S[i-1] \in S
+findFixpointOfDifferences := proc(L)
+  local res, i:
+  res := L:
+  for i from 2 to nops(L) do
+    res := { op(res), L[i]-L[i-1] }:
+  od:
+  if (res=L) then return res else return findFixpointOfDifferences(res) fi:
+end proc:
+
+
+#
+# error_counter functions allows one to follow the accumulation of errors in each variable.
+#   an error_counter is a list of the form [[var1, c1], [var2, c2], ... ]
+#   where the vari are variable names and the ci indicate how many approximation errors
+#   are accumulated in vari.
+#
+
+#
+# This procedure initializes the counter associated with variable var to 1 (and creates it if needed.)
+# It returns an up-to-date error_counter.
+init_error_counter := proc (var, error_counter)
+  local i, res:
+  res  := error_counter:
+  for i from 1 to nops(res) do
+    if (res[i][1]=var) 
+    then res[i][2] := 1:
+         return res:
+    fi
+  od:
+  res := [op(res), [var, 1]]:
+end:
+
+
+#
+# This procedure adds a given number to the counter associated with variable var.
+# It returns an up-to-date error_counter.
+add_to_error_counter := proc (var, n, error_counter)
+  local i, res:
+  res  := error_counter:
+  for i from 1 to nops(res) do
+    if (res[i][1]=var) 
+    then res[i][2] := res[i][2]+n:
+         return res:
+    fi
+  od:
+  res := [op(res), [var, n]]:
+end proc:
+
+#
+# This procedure sets the value of the counter associated with variable var.
+# It returns an up-to-date error_counter.
+set_error_counter := proc(var, n, error_counter)
+  local i,err:
+  err  := error_counter:
+  for i from 1 to nops(err) do
+    if (err[i][1]=var) 
+    then err[i][2] := n:
+         return err:
+    fi
+  od:
+  err := [op(err), [var, n]]:  
+end proc:
+
+#
+# This procedure initializes the counter associated to the multiplication of var2 and var3,
+# putting the result in variable var1. 
+# It returns an up-to-date error_counter.
+error_counter_of_a_multiplication := proc (var1, var2, var3, error_counter)
+  local i, res, c2, c3:
+  c2 := 0: c3 := 0:
+  for i from 1 to nops(error_counter) do
+    if (error_counter[i][1]=var2) then c2 := error_counter[i][2] fi:
+    if (error_counter[i][1]=var3) then c3 := error_counter[i][2] fi:
+    if (error_counter[i][1]=var1)
+    then
+      res := [ op(error_counter[1..i-1]), op(error_counter[i+1..nops(error_counter)]) ]
+    fi:
+  od:
+  if (res = 'res') then res := error_counter fi:
+  res := [op(res), [var1, c2+c3+1]]:
+end:
+
+#
+# Copies the error counter of var2 into var1
+error_counter_on_copy := proc(var1, var2, error_counter)
+  local i, err, c2:
+  c2 := 0:
+  for i from 1 to nops(error_counter) do
+    if (error_counter[i][1] = var2) then c2 := error_counter[i][2] fi:
+    if (error_counter[i][1] = var1)
+    then
+      err := [ op(error_counter[1..i-1]), op(error_counter[i+1..nops(error_counter)]) ]
+    fi:
+  od:
+  if (err = 'err') then err := error_counter fi:
+  if (c2 <> 0) then err := [op(res), [var1, c2]] fi:
+end proc:
+
+
+#
+# Returns the value of the error counter associated to a variable
+find_in_error_counter := proc(var, error_counter) 
+  local i:
+  for i from 1 to nops(error_counter) do
+    if (error_counter[i][1] = var) then return error_counter[i][2] fi:
+  od:
+  return 0:
+end proc:
+
+#
+# generate_multiply_rational(fd, var1, var2, r, error_counter, indent) generates code for performing
+# var1 = var2*r in MPFR
+#   fd is the file descriptor in which the code shall be produced.
+#   var1 and var2 are strings representing variable names. r is a Maple rational number.
+#   error_counter is an error_counter (as described above).
+#   indent is an optional argument. It is a string used to correctly indent the code. It is prefixed to any
+#   generated line. Hence, if indent="  ", the generated code will be indented by 2 spaces. 
+# An up-to-date error_counter is returned.
+generate_multiply_rational := proc(fd, var1, var2, r, error_counter, indent:="")
+  local p,q,err:
+  err := error_counter:
+  if (whattype(r)<>'fraction') and (whattype(r)<>'integer')
+  then error "generate_multiply_rational used with non rational number %1", r: fi:
+  if (abs(r)=1)
+  then
+    if (var1=var2)
+    then
+      if (r<>1) then fprintf(fd, "%sMPFR_CHANGE_SIGN (%s);\n", indent, var1) fi:
+      return err:
+    else 
+      if (r=1)
+        then fprintf(fd, "%smpfr_set (%s, %s, MPFR_RNDN);\n", indent, var1, var2):
+        else fprintf(fd, "%smpfr_neg (%s, %s, MPFR_RNDN);\n", indent, var1, var2):
+      fi:
+      return error_counter_on_copy(var1, var2, err):
+    fi
+  fi:
+  # Now, r is a rational number different from 1.
+  p := numer(r): q := denom(r):
+  if (abs(p)<>1)
+  then
+    fprintf(fd, "%smpfr_mul_si (%s, %s, %d, MPFR_RNDN);\n", indent, var1, var2, p):
+    err := error_counter_of_a_multiplication(var1, var2, "", err):
+    if(q<>1)
+    then
+      fprintf(fd, "%smpfr_div_si (%s, %s, %d, MPFR_RNDN);\n", indent, var1, var1, q):
+      err := error_counter_of_a_multiplication(var1, var1, "", err):
+    fi:
+  else
+    fprintf(fd, "%smpfr_div_si (%s, %s, %d, MPFR_RNDN);\n", indent, var1, var2, p*q):
+    err := error_counter_of_a_multiplication(var1, var2, "", err):
+  fi:
+  return err:
+end proc:
+
+
+#
+# generate_multiply_poly is the same as generate_multiply_rational but when r is a rational fraction.
+# The fraction r must have the form p/q where p and q are polynomials with integer coefficients.
+# Moreover, the gcd of the coefficients of p must be 1. Idem for q.
+# The procedure returned a list [m, d, err] where m is the set of indices k such that
+# a mpfr_mul_sik function is needed and idem for d with mpfr_div_sik.
+# err is an up-to-date error counter.
+generate_multiply_poly := proc(fd, var1, var2, r, error_counter, indent:="")
+  local p,q,Lp,Lq,n,i,j,var, required_mulsi, required_divsi, err:
+  err := error_counter:
+  required_mulsi := {}:
+  required_divsi := {}:
+  p := numer(r): q := denom(r):
+  Lp := my_factors(p): Lq := my_factors(q):
+  if (Lp[1] <> 1)
+    then error "generate_multiply_poly: an integer can be factored out of %1", p:
+  fi:
+  if (Lq[1] <> 1)
+    then error "generate_multiply_poly: an integer can be factored out of %1", q:
+  fi:
+  Lp := Lp[2]: Lq := Lq[2]:
+  var := var2:
+  if (nops(Lp) <> 0)
+  then
+    n := 0:
+    for i from 1 to nops(Lp) do n := n + Lp[i][2] od:
+    if (n=1)
+    then
+      fprintf(fd, "%smpfr_mul_si (%s, %s", indent, var1, var):
+    else
+      required_mulsi := { op(required_mulsi), n }:
+      fprintf(fd, "%smpfr_mul_si%d (%s, %s", indent, n, var1, var):
+    fi:
+    for i from 1 to nops(Lp) do
+      for j from 1 to Lp[i][2] do
+        fprintf(fd, ", %a", Lp[i][1]):
+      od:
+    od:
+    fprintf(fd, ", MPFR_RNDN);\n"):
+    err := set_error_counter(var1, n+find_in_error_counter(var, err) , err):
+    var := var1:
+  fi:
+  if (nops(Lq) <> 0)
+  then
+    n := 0:
+    for i from 1 to nops(Lq) do n := n + Lq[i][2] od:
+    if (n=1)
+    then
+      fprintf(fd, "%smpfr_div_si (%s, %s", indent, var1, var):
+    else
+      required_divsi := { op(required_divsi), n }:
+      fprintf(fd, "%smpfr_div_si%d (%s, %s", indent, n, var1, var)
+    fi:
+    for i from 1 to nops(Lq) do
+      for j from 1 to Lq[i][2] do
+        fprintf(fd, ", %a", Lq[i][1])
+      od:
+    od:
+    fprintf(fd, ", MPFR_RNDN);\n"):
+    err := set_error_counter(var1, n+find_in_error_counter(var, err) , err):
+    var := var1:
+  fi:
+  if (var1 <> var) then
+    fprintf(fd, "%smpfr_set (%s, %s, MPFR_RNDN);\n", indent, var1, var):
+    err := set_error_counter(var1, find_in_error_counter(var, err) , err):
+  fi:
+  return [required_mulsi, required_divsi, err]:
+end proc:
+
+
+#
+# This function generates the code of a procedure mpfr_mul_uin or mpfr_div_uin
+#
+generate_muldivsin := proc(op, n)
+  local i, var:
+  if ((op <> "mul") and (op <> "div"))
+    then error "Invalid argument to generate_muldivuin (%1). Must be 'mul' or 'div'", op
+  fi:
+  if (whattype(n) <> 'integer')
+  then error "Invalid argument to generate_muldivuin (%1). Must be an integer.", n
+  fi:
+
+  if (op="mul") then var := "MUL" else var := "DIV" fi:
+
+  printf("__MPFR_DECLSPEC void mpfr_div_si%d _MPFR_PROTO((mpfr_ptr, mpfr_srcptr,\n", n):
+  for i from n to 2 by -2 do
+    printf("                                               long int, long int,\n"):
+  od:
+  if (i=1)
+  then
+    printf("                                               long int, mpfr_rnd_t));\n"):
+  else
+    printf("                                               mpfr_rnd_t));\n")
+  fi:
+
+  printf("\n\n\n"):
+  printf("void\n"):
+  printf("mpfr_%s_si%d (mpfr_ptr y, mpfr_srcptr x,\n", op, n):
+  for i from n to 2 by -2 do
+    printf("              long int v%d, long int v%d,\n", n-i+1, n-i+2):
+  od:
+  if (i=1)
+  then
+    printf("              long int v%d, mpfr_rnd_t mode)\n", n):
+  else
+    printf("              mpfr_rnd_t mode)\n")
+  fi:
+  printf("{\n"):
+  printf("  long int acc = v1;\n"):
+  printf("  mpfr_set (y, x, mode);\n"):
+   for i from 2 to n do
+    printf("  MPFR_ACC_OR_%s (v%d);\n", var, i):
+  od:
+  printf("  mpfr_%s_si (y, y, acc, mode);\n", op):
+  printf("}\n"):
+  return:
+end proc:
diff --git a/tools/metaMPFR/metaMPFR_straightforwardAlgo.mpl b/tools/metaMPFR/metaMPFR_straightforwardAlgo.mpl
new file mode 100644
index 000000000..3a86f87a7
--- /dev/null
+++ b/tools/metaMPFR/metaMPFR_straightforwardAlgo.mpl
@@ -0,0 +1,703 @@
+read("metaMPFR_common.mpl"):
+
+FUNCTION_SERIES := 0:
+FUNCTION_SERIES_RATIONAL := %+1:
+CONSTANT_SERIES := %+1:
+
+###############################################################################
+######################### We can now generate the code ########################
+###############################################################################
+# This procedure generate code for a straightforward evaluation of the series
+# corresponding to a recurrence.
+#   rec is a recurrence with inital conditions
+#   type can take three values, depending on which series should be evaluated:
+#     FUNCTION_SERIES -> produces code for evaluating sum( a(i)*x^i )
+#     FUNCTION_SERIES_RATIONAL -> produces code for evaluating sum( a(i)*(p/q)^i )
+#     CONSTANT_SERIES -> produces code for evaluating sum( a(i) )
+#   name is the name the should be given to the produced procedure.
+#   fofx is an optional parameter. It is the function being implemented.
+#     -> if provided, this argument will be used to (heuristically) find the limits
+#        of the function at +/-oo and find its asymptotical behavior.
+#   *IMPORTANT NOTE*: it must be a function of the variable 'x'.
+#                     Moreover, neither _f nor x shall be assigned at the time of
+#                     calling the function
+#    
+generateStraightforwardAlgo := proc(rec, aofn, type, name, filename, fofx := _f(x))
+  local a, b, n, fd, init_cond, nc, d, exponents, formalRec, c, s1, s2, p, q, hardconstant, f0, i0, ri0, i, j, var, var1, var2, var3, guard_bits, required_mulsi, required_divsi, temp, error_counter, error_in_loop, maxofti, additional_error:
+
+  getname(aofn, a, n):
+
+  required_mulsi := {}:
+  required_divsi := {}:
+  fd := fopen(filename, WRITE):
+
+  # We check that we have a recurrence of the form a(n+d)=r(n)*a(n)
+  # and we extract its canonical form:
+  #    a(n+d) = c * s1(n)/s1(n+d) * s2(n+d)/s2(n) * p(n)/q(n) * a(n)
+  #
+  formalRec := checkOneTermRecurrence(rec, a(n)):
+  decomposeOneTermRecurrence(formalRec, c, s1, s2, p, q):
+  d := formalRec:-order:
+
+  # We keep only non-trivial inital conditions
+  init_cond := removeTrivialConditions(rec, a(n));
+  nc := nops(init_cond):
+  exponents := {1, d}:
+  for i from 1 to nc do
+    exponents := { op(exponents), init_cond[i][1] }
+  od:
+  exponents := findFixpointOfDifferences(exponents):
+  exponents := exponents minus {0}:
+
+  error_counter := []:
+
+  fprintf(fd, "/* Evaluation by a straightforward algorithm */\n"):
+  fprintf(fd, "/* Code automatically generated by metaMPFR. */\n"):
+
+  fprintf(fd, "static int\n"):
+  if (type=FUNCTION_SERIES)
+    then fprintf(fd, "mpfr_%a (mpfr_ptr res, mpfr_srcptr x, mpfr_rnd_t rnd)\n", name):
+  elif (type=FUNCTION_SERIES_RATIONAL)
+    then fprintf(fd, "mpfr_%a (mpfr_ptr res, int u, int v, mpfr_rnd_t rnd)\n", name):
+  elif (type=CONSTANT_SERIES)
+    then fprintf(fd, "mpfr_%a (mpfr_ptr res, mpfr_rnd_t rnd)\n", name):
+  fi:
+  fprintf(fd, "{\n"):
+
+
+  ######################################################
+  #################### Declarations ####################
+  ######################################################
+
+  fprintf(fd, "  MPFR_ZIV_DECL (loop);\n"):
+  fprintf(fd, "  MPFR_SAVE_EXPO_DECL (expo);\n"):
+  fprintf(fd, "  mpfr_prec_t wprec;             /* working precision */\n"):
+  fprintf(fd, "  mpfr_prec_t prec;              /* target precision */\n"):
+  fprintf(fd, "  mpfr_prec_t err;               /* used to estimate the evaluation error */\n"):
+  fprintf(fd, "  mpfr_prec_t correctBits;       /* estimates the number of correct bits*/\n"):
+  fprintf(fd, "  unsigned long int k;\n"):
+  fprintf(fd, "  unsigned long int conditionNumber;        /* condition number of the series */\n"):
+  fprintf(fd, "  unsigned assumedExp;           /* used as a lowerbound of -EXP(f(x)) */\n"):
+  fprintf(fd, "  int r;                         /* returned ternary value */\n"):
+  fprintf(fd, "  mpfr_t s;                      /* used to store the partial sum */\n"):
+  
+  if (whattype(c) = 'fraction') or (whattype(c) = 'integer')
+    then hardconstant := 0
+    else hardconstant := 1
+  fi:
+  if (type=FUNCTION_SERIES) then hardconstant := 1 fi:
+
+  if (hardconstant=1)
+  then
+    if (type=CONSTANT_SERIES)
+      then fprintf(fd, "  mpfr_t x%d;                      /* used to store %a */\n", d, c):
+    elif (type=FUNCTION_SERIES_RATIONAL)
+      then fprintf(fd, "  mpfr_t x%d;                      /* used to store %a */\n", d, c*(u/v)^d):
+    elif (type=FUNCTION_SERIES)
+      then fprintf(fd, "  mpfr_t x%d;                      /* used to store %a */\n", d, c*x^d):
+    fi:
+  fi:
+
+  if (type=FUNCTION_SERIES)
+  then
+    fprintf(fd, "  mpfr_t tmp;\n"):
+    if (nops(exponents)-1 >= 1) then fprintf(fd, "  mpfr_t ") fi:
+    for i from 1 to nops(exponents)-1 do
+      fprintf(fd, "x%d", exponents[i]):
+      if (i<nops(exponents)-1) then fprintf(fd, ", ") else fprintf(fd, ";                     /* used to store x^i */\n") fi
+    od:
+  fi:
+
+  fprintf(fd, "  mpfr_t ");
+  for i from 1 to nc do
+    fprintf(fd, "tip%d", init_cond[i][1]):
+    if (i<nc) then fprintf(fd, ", ") else fprintf(fd, "; /* used to store successive values of t_i */\n") fi:
+  od:
+  fprintf(fd, "  int "):
+  for i from 1 to nc do
+    fprintf(fd, "test%d", init_cond[i][1]):
+    if (i<nc) then fprintf(fd, ", ") else fprintf(fd, ";\n") fi:
+  od:
+
+  fprintf(fd, "\n"):
+  fprintf(fd, "  /* Logging */\n"):
+  if (type=FUNCTION_SERIES)
+    then fprintf(fd, "  MPFR_LOG_FUNC ( (\"x[%%#R]=%%R rnd=%%d\", x, x, rnd), (\"res[%%#R]=%%R\", res, res) );\n\n")
+  elif (type=FUNCTION_SERIES_RATIONAL)
+    then fprintf(fd, "  MPFR_LOG_FUNC ( (\"x=u/v with u=%%d and v=%%d, rnd=%%d\", u, v, rnd), (\"res[%%#R]=%%R\", res, res) );\n\n")
+  else fprintf(fd, "  MPFR_LOG_FUNC ( (\"rnd=%%d\", rnd), (\"res[%%#R]=%%R\", res, res) );\n\n")
+  fi:
+
+
+  ######################################################
+  #################### Special cases ###################
+  ######################################################
+
+  if ( (type=FUNCTION_SERIES) or (type=FUNCTION_SERIES_RATIONAL) )
+  then
+    fprintf(fd, "  /* Special cases */\n"):
+ 
+    if (type=FUNCTION_SERIES)
+      then fprintf(fd, "  if (MPFR_UNLIKELY (MPFR_IS_NAN (x)))\n")
+      else fprintf(fd, "  if (MPFR_UNLIKELY (v==0))\n")
+    fi:
+    fprintf(fd, "    {\n"):
+    fprintf(fd, "      MPFR_SET_NAN (res);\n"):
+    fprintf(fd, "      MPFR_RET_NAN;\n"):
+    fprintf(fd, "    }\n"):
+
+    if (init_cond[1][1] > 0) then f0 := 0 else f0 := init_cond[1][2] fi:
+    if (type=FUNCTION_SERIES)
+      then fprintf(fd, "  if (MPFR_UNLIKELY (MPFR_IS_ZERO (x)))\n")
+      else fprintf(fd, "  if (MPFR_UNLIKELY (u==0))\n")
+    fi:
+    if (whattype(f0) = 'integer')
+    then
+      fprintf(fd, "    {\n"):
+      fprintf(fd, "      return mpfr_set_si (res, %a, rnd);\n", f0):
+      fprintf(fd, "    }\n")
+    else
+      printf("You need to provide a function mpfr_%a0(mpfr_t, mpfr_rnd_t) that evaluates %a with correct rounding\n", name, f0):
+      fprintf(fd, "    {\n"):
+      fprintf(fd, "      return mpfr_%a0 (res, rnd);\n", name):
+      fprintf(fd, "    }\n")
+    fi:
+
+    if (type=FUNCTION_SERIES)
+    then
+      fprintf(fd, "  if (MPFR_UNLIKELY (MPFR_IS_INF (x)))\n"):
+      fprintf(fd, "    {\n"):
+      for i from -1 to 1 by 2 do  # Trick to handle both -oo and +oo
+        if (i<0)
+          then fprintf(fd, "      if (MPFR_IS_NEG (x))\n"):
+          else fprintf(fd, "      else\n")
+        fi:
+        fprintf(fd, "        {\n"):
+        if (i<0) then f0 := "m" else f0 := "p" fi:
+        if (fofx <> _f(x))
+        then
+          f0 := limit(fofx, x=i*infinity):
+          if (whattype(f0)='integer')
+          then fprintf(fd, "          return mpfr_set_si (res, %a, rnd);\n", f0):
+          elif (f0 = infinity) or (f0 = -infinity)
+          then
+            fprintf(fd, "          MPFR_SET_INF(res);\n"):
+            if (f0>0)
+            then fprintf(fd, "          MPFR_SET_POS(res);\n")
+            else fprintf(fd, "          MPFR_SET_NEG(res);\n")
+            fi:
+            fprintf(fd, "          MPFR_RET(0);\n"):
+          else
+            if (i<0) then f0 := "m" else f0 := "p" fi:
+          fi: 
+        fi:
+        if ((f0 = "p") or (f0 = "m"))
+        then
+          printf("You need to provide a function mpfr_%a%sinf(mpfr_t, mpfr_rnd_t) that evaluates lim(f(x), x=", name, f0):
+          if (f0 = "m") then printf("-") fi:
+          printf("inf) with correct rounding.\n"):
+          fprintf(fd, "          return mpfr_%a%sinf (res, rnd);\n", name, f0)
+        fi:
+        fprintf(fd, "        }\n")
+      od:
+      fprintf(fd, "    }\n")
+    fi:
+  fi:
+
+
+  ######################################################
+  ################## Precomputations ###################
+  ######################################################
+
+  fprintf(fd, "\n"):
+  fprintf(fd, "  /* Save current exponents range */\n"):
+  fprintf(fd, "  MPFR_SAVE_EXPO_MARK (expo);\n\n"):
+
+  if ( (type=FUNCTION_SERIES) or (type=FUNCTION_SERIES_RATIONAL) )
+  then  fprintf(fd, "  /* FIXME: special case for large values of |x| ? */\n\n")
+  fi:
+
+  # Note : prec is the value such that we will try to compute an approximation
+  # with relative error smaller than 2^(1-prec).
+  # Several things may happen:
+  #   1) We do not achieve the intended error: this is because we badly estimated the exponent of the result
+  #   2) We achieve the error, but it is not sufficient to decide correct rounding (Ziv's bad case)
+  # Against 1), we can try to make our estimate of the exponent better with any heuristic.
+  # Against 2), we can consider more guard bits. 11 guard bits seem a good value for the beginning
+  #   (statistically, we expect to fail in less than 0.1 % of the cases)
+  # wprec is the precision used during the computation, in order to ensure the final relative error 2^(1-prec)
+  #
+  fprintf(fd, "  /* We begin with 11 guard bits */\n"):
+  fprintf(fd, "  prec = MPFR_PREC (res) + 11;\n"):
+  fprintf(fd, "  MPFR_ZIV_INIT (loop, prec);\n"):
+
+  # TODO: the value here can be chosen completely heursitically. We could do something much better
+  # when fofx is known by using asympt(fofx, x, 1). A clean implementation appears to be complex though.
+  # We must catch errors if the development does not exist (e.g. AiryAi(-x));
+  # We must find a separation after which the asymptotic behavior is valid (e.g. x>1)
+  printf("The code contains a variable assumedExp arbitrarily set to 10. You can put any value heuristically chosen. The closer it is to -log_2(|f(x)|), the better it is.\n"):
+  fprintf(fd, "  assumedExp = 10; /* TIP: You can put any value heuristically chosen. The closer it is to -log_2(|f(x)|), the better it is */\n"):
+
+  # TODO: find a way of putting a rigorous value here.
+  # This value *must* be rigorous: the safety of the implementation relies on it.
+  # Precisely, we need to have sum(|a(i)*x^i|) <= 2^conditionNumber
+  fprintf(fd, "  conditionNumber = xxx; /* FIXME: set a value such that sum(|a(i)*x^i|) <= 2^conditionNumber */\n"):
+  printf("The code contains a variable conditionNumber that you must manually set to a suitable value, in order to ensure that sum_{i=0}^{infinity} |a(i)*x^i| <= 2^conditionNumber\n"):
+  fprintf(fd, "  wprec = prec + ERRORANALYSISPREC + conditionNumber + assumedExp;\n"):
+
+
+  ######################################################
+  ################## Initialisations ###################
+  ######################################################
+
+  if (hardconstant=1)
+  then
+    fprintf(fd, "  mpfr_init (x%d);\n", d):
+  fi:
+  if (type = FUNCTION_SERIES)
+  then
+    fprintf(fd, "  mpfr_init (tmp);\n"):
+    for i from 1 to nops(exponents)-1 do
+      fprintf(fd, "  mpfr_init (x%d);\n", exponents[i])
+    od:
+  fi:
+
+  for i from 1 to nc do
+    fprintf(fd, "  mpfr_init (tip%d);\n", init_cond[i][1])
+  od:
+  fprintf(fd, "  mpfr_init (s);\n\n"):
+
+
+  ######################################################
+  ########## Ziv' loop: setting the precision ##########
+  ######################################################
+
+  fprintf(fd, "  /* ZIV' loop */\n"):
+  fprintf(fd, "  for (;;)\n"):
+  fprintf(fd, "    {\n"):
+
+  fprintf(fd, "      MPFR_LOG_MSG ((\"Working precision: %%d\\n\", wprec, 0));\n\n"):
+  if (hardconstant=1)
+  then
+    fprintf(fd, "      mpfr_set_prec (x%d, wprec);\n", d):
+  fi:
+  if (type = FUNCTION_SERIES)
+  then
+    fprintf(fd, "      mpfr_set_prec (tmp, wprec);\n"):
+    fprintf(fd, "      if(mpfr_get_prec(x) > wprec)\n"):
+    fprintf(fd, "        mpfr_set_prec (x1, wprec);\n"):
+    fprintf(fd, "      else\n"):
+    fprintf(fd, "        mpfr_set_prec(x1, mpfr_get_prec(x));\n"):
+    for i from 2 to nops(exponents)-1 do
+      fprintf(fd, "      mpfr_set_prec (x%d, wprec);\n", exponents[i])
+    od:
+  fi:
+
+  for i from 1 to nc do
+    fprintf(fd, "      mpfr_set_prec (tip%d, wprec);\n", init_cond[i][1])
+  od:
+  fprintf(fd, "      mpfr_set_prec (s, wprec);\n\n"):
+
+
+  ######################################################
+  ############ Ziv' loop: initial conditions ###########
+  ######################################################
+
+  fprintf(fd, "      mpfr_set_ui (s, 0, MPFR_RNDN);\n"):
+
+  if (type = FUNCTION_SERIES)
+  then
+    fprintf(fd, "      mpfr_set (x1, x, MPFR_RNDN);\n"):
+    error_counter := init_error_counter("x1", error_counter):
+    for i from 2 to nops(exponents) do
+        fprintf(fd, "      mpfr_mul (x%d, x%d, x%d, MPFR_RNDN);\n", exponents[i], exponents[i-1], exponents[i]-exponents[i-1]):
+	var1 := sprintf("x%d", exponents[i]):
+        var2 := sprintf("x%d", exponents[i-1]):
+        var3 := sprintf("x%d", exponents[i]-exponents[i-1]):
+	error_counter := error_counter_of_a_multiplication(var1, var2, var3, error_counter):
+    od:
+  fi:
+
+  for i from 1 to nc do
+    i0 := init_cond[i][1]:
+    ri0 := init_cond[i][2]:  # We implement t_{i0} <- ri0
+    if (whattype(ri0)='integer') or (whattype(ri0)='fraction')
+    then
+      if (type = FUNCTION_SERIES) and (i0 <> 0)
+        then 
+          var :=  sprintf("      mpfr_mul_si (tip%d, x%d, ", i0, i0):
+          var1 := sprintf("tip%d", i0):
+          var2 := sprintf("x%d", i0):
+          error_counter := error_counter_of_a_multiplication(var1, var2, "", error_counter):
+        else 
+          var := sprintf("      mpfr_set_si (tip%d, ", i0):
+          var1 := sprintf("tip%d", i0):
+          error_counter := init_error_counter(var1, error_counter):
+      fi:
+      fprintf(fd, "%s%d, MPFR_RNDN);\n", var, numer(ri0)):
+      if (whattype(ri0)='fraction')
+      then
+        fprintf(fd, "      mpfr_div_si (tip%d, tip%d, %d, MPFR_RNDN);\n", i0, i0, denom(ri0)):
+        var1 := sprintf("tip%d", i0):
+        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
+      fi:
+    else
+      printf("You need to provide a function mpfr_%a%d (mpfr_t, mpfr_rnd_t) that evaluates %a with faithful rounding.\n", name, i0, ri0): 
+      fprintf(fd, "      mpfr_%a%d (tip%d, MPFR_RNDN);\n", name, i0, i0):
+      var1 := sprintf("tip%d", i0):
+      error_counter := init_error_counter(var1, error_counter):
+      if (type = FUNCTION_SERIES) and (i0 <> 0)
+      then
+        fprintf(fd, "      mpfr_mul (tip%d, tip%d, x%d, MPFR_RNDN);\n", i0, i0, i0):
+        var1 := sprintf("tip%d", i0):
+        var2 := sprintf("x%d", i0):
+        error_counter := error_counter_of_a_multiplication(var1, var1, var2, error_counter):
+      fi:
+    fi:
+
+    if (type = FUNCTION_SERIES_RATIONAL) and (i0 <> 0)
+    then
+      var1 := sprintf("tip%d", i0):
+      if (i0 = 1)
+      then fprintf(fd, "      mpfr_mul_si (tip%d, tip%d, ", i0, i0):
+      else
+        required_mulsi := { op(required_mulsi), i0 }:
+        fprintf(fd, "      mpfr_mul_si%d (tip%d, tip%d, ", i0, i0, i0):
+      fi:
+      for j from 1 to i0 do
+        fprintf(fd, "u, "):
+        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
+      od:
+      fprintf(fd, "MPFR_RNDN);\n"):
+      if (i0 = 1)
+      then fprintf(fd, "      mpfr_div_si (tip%d, tip%d, ", i0, i0):
+      else
+        required_divsi := { op(required_divsi), i0 }:
+        fprintf(fd, "      mpfr_div_si%d (tip%d, tip%d, ", i0, i0, i0):
+      fi:
+      for j from 1 to i0 do
+        fprintf(fd, "v, "):
+        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
+      od:
+      fprintf(fd, "MPFR_RNDN);\n"):
+    fi:
+
+    fprintf(fd, "      mpfr_add (s, s, tip%d, MPFR_RNDN);\n\n", i0):
+  od:
+
+  if (whattype(c) = 'integer') or (whattype(c) = 'fraction')
+  then
+    if (type = FUNCTION_SERIES) and (c = -1)
+    then
+      fprintf(fd, "      MPFR_CHANGE_SIGN (x%d);\n", d):
+    elif (type = FUNCTION_SERIES) and (c <> 1)
+    then
+      fprintf(fd, "      mpfr_mul_si (x%d, x%d, %d, MPFR_RNDN);\n", d, d, numer(c)):
+      var1 := sprintf("x%d", d):
+      error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
+      if (whattype(c) = 'fraction')
+      then
+        fprintf(fd, "      mpfr_div_si (x%d, x%d, %d, MPFR_RNDN);\n", d, d, denom(c)):
+        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
+      fi:
+    fi:
+  else
+    printf("You need to provide a function mpfr_%a_cste (mpfr_t, mpfr_rnd_t) that evaluates %a with faithful rounding.\n", name, c):
+    if (type = CONSTANT_SERIES)
+    then
+      fprintf(fd, "      mpfr_%a_cste (x%d, MPFR_RNDN);\n", name, d):
+      var1 := sprintf("x%d", d):
+      error_counter := init_error_counter(var1, error_counter):
+    elif (type = FUNCTION_SERIES) then
+      fprintf(fd, "      mpfr_%a_cste (tmp, MPFR_RNDN);\n", name):
+      error_counter := init_error_counter("tmp", error_counter):
+      fprintf(fd, "      mpfr_mul (x%d, tmp, x%d, MPFR_RNDN);\n", d, d):
+      var1 := sprintf("x%d", d):
+      error_counter := error_counter_of_a_multiplication(var1, "tmp", var1, error_counter):
+    elif (type = FUNCTION_SERIES_RATIONAL) then
+      fprintf(fd, "      mpfr_%a_cste (x%d, MPFR_RNDN);\n", name, d):
+      var1 := sprintf("x%d", d):
+      error_counter := init_error_counter(var1, error_counter):
+      if (d = 1)
+      then fprintf(fd, "      mpfr_mul_si (x%d, x%d, ", d, d):
+      else
+        required_mulsi := { op(required_mulsi), d }:
+        fprintf(fd, "      mpfr_mul_si%d (x%d, x%d, ", d, d, d):
+      fi:
+      for j from 1 to d do
+        fprintf(fd, "u, "):
+        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
+      od:
+      fprintf(fd, "MPFR_RNDN);\n"):
+      var1 := sprintf("x%d", d):
+      if (d = 1)
+      then fprintf(fd, "      mpfr_div_si (x%d, x%d, ", d, d):
+      else
+        required_divsi := { op(required_divsi), d }:
+        fprintf(fd, "      mpfr_div_si%d (x%d, x%d, ", d, d, d):
+      fi:
+      for j from 1 to d do
+        fprintf(fd, "v, "):
+        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
+      od:
+      fprintf(fd, "MPFR_RNDN);\n"): 
+    fi:
+  fi:
+
+  ######################################################
+  ######### Ziv' loop: evaluation of the series ########
+  ######################################################
+
+  fprintf(fd, "\n"):
+  fprintf(fd, "      /* Evaluation of the series */\n"):
+  fprintf(fd, "      k = %d;\n", d):
+  fprintf(fd, "      for (;;)\n"):
+  fprintf(fd, "        {\n"):
+  if (init_cond[1][1] <> 0) then fprintf(fd, "          k += %d;\n", init_cond[1][1]) fi:
+
+  for i from 1 to nc do
+    error_in_loop := 0:
+    i0 := init_cond[i][1]:
+    if (hardconstant = 1)
+    then 
+      fprintf(fd, "          mpfr_mul (tip%d, tip%d, x%d, MPFR_RNDN);\n", i0, i0, d):
+      var1 := sprintf("x%d", d):
+      error_in_loop := error_in_loop + 1 + find_in_error_counter(var1, error_counter):
+    else
+      var := sprintf("tip%d", i0):
+      temp := generate_multiply_rational(fd, var, var, c, [[var, error_in_loop]], "          "):
+      error_in_loop := find_in_error_counter(var, temp):
+      if (type = FUNCTION_SERIES_RATIONAL)
+      then
+        if (d=1)
+        then
+          fprintf(fd, "          mpfr_mul_si (tip%d, tip%d, u, MPFR_RNDN);\n", i0, i0):
+          fprintf(fd, "          mpfr_div_si (tip%d, tip%d, v, MPFR_RNDN);\n", i0, i0):
+	  error_in_loop := error_in_loop + 2:
+        else
+          required_mulsi := { op(required_mulsi), d }:
+          fprintf(fd, "          mpfr_mul_si%d (tip%d, tip%d", d, i0, i0):
+          for j from 1 to d do fprintf(fd, ", u") od:
+          fprintf(fd, ", MPFR_RNDN);\n"):
+	  error_in_loop := error_in_loop + d:
+
+          required_divsi := { op(required_divsi), d }:
+          fprintf(fd, "          mpfr_div_si%d (tip%d, tip%d", d, i0, i0):
+          for j from 1 to d do fprintf(fd, ", v") od:
+          fprintf(fd, ", MPFR_RNDN);\n" ):
+	  error_in_loop := error_in_loop + d:
+        fi
+      fi
+    fi:
+    var := sprintf("tip%d", i0):
+
+    temp := generate_multiply_poly(fd, var, var, subs(n=k-d, p/q), [[var, error_in_loop]], "          "):
+    required_mulsi := { op(required_mulsi), op(temp[1]) }:
+    required_divsi := { op(required_divsi), op(temp[2]) }:
+    error_in_loop := find_in_error_counter(var, temp[3]):
+    temp := generate_multiply_poly(fd, "tmp", var, subs(n=k, s2/s1), [[var, error_in_loop]], "          "):
+    required_mulsi := { op(required_mulsi), op(temp[1]) }:
+    required_divsi := { op(required_divsi), op(temp[2]) }:
+
+    fprintf(fd, "          mpfr_add (s, s, tmp, MPFR_RNDN);\n"):
+    
+    if (i<nc) then fprintf(fd, "\n          k += %d;\n", init_cond[i+1][1]-i0)
+    else  fprintf(fd, "\n          k += %d;\n", d-i0)
+    fi:
+  od:
+
+
+  ######################################################
+  ################### Error analysis ###################
+  ######################################################
+
+  maxofti := 0: # store the maximum of the error counters of the initial conditions
+  for i from 1 to nc do
+    var := sprintf("tip%d", init_cond[i][1]):
+    if find_in_error_counter(var, error_counter) > maxofti
+      then maxofti := find_in_error_counter(var, error_counter):
+    fi:
+  od:
+  additional_error := find_in_error_counter("tmp", temp[3]) - error_in_loop:
+
+
+  ######################################################
+  #### Ziv' loop: stopping criterion for the series ####
+  ######################################################
+
+  # The first neglected term is tk, so the remainder is made by
+  # tk + t(k+d) + t(k+2d).... and the corresponding series
+  # beginning with t(k+1), t(k+2), etc. up to t(k+d-1).
+  #
+  # We have t(k0+d) = c*s1(k0)/s1(k0+d) * s2(k0+d)/s2(k0)* p(k0)/q(k0) * x^d t(k0)
+  # (where x=u/v or x=1 in cases of rational series or constant series)
+  # So it suffices that:
+  # forall k0>=k, |c * s1(k0)/s1(k0+d) * s2(k0+d)/s2(k0) * p(k0)/q(k0) * x^d| <= 1/2     (1)
+  #
+  # If this is true, we can bound |tk + t(k+d) + t(k+2d) + ...| by 2tk
+  # so the total remainder is bounded by 2*nc*tk.
+  #
+  # global_test depends on k and we must satisfy:
+  #  "if (global_test) then (1) holds".
+
+  fprintf(fd, "          global_test = xxx; /* FIXME: set the value in order to ensure that, whenever global_test is true, we have: forall k'>=k, |r(k')*x^d| <= 1/2, where r is the fraction such that a(n+d)=r(n)a(n)*/\n"):
+  printf("The code contains a variable global_test that you must manually set to a suitable value, in order to ensure that when global_test is true, the following holds:\n"):
+  printf("        forall k'>=k, |r(k')*x^d| <= 1/2, where r is the fraction such that a(n+d)=r(n)a(n)\n"):
+  guard_bits := 1+1+ceil(log[2](nc)):
+  for i from 1 to nc do
+    i0 := init_cond[i][1]:
+    fprintf(fd, "          test%d = ( (!MPFR_IS_ZERO(s))\n", i0):
+    fprintf(fd, "                    && ( MPFR_IS_ZERO(tip%d)\n", i0):
+    fprintf(fd, "                         || (MPFR_EXP(tip%d) + (mp_exp_t)prec + %d <= MPFR_EXP(s))\n", i0, guard_bits):
+    fprintf(fd, "                         )\n"):
+    fprintf(fd, "                    );\n"):
+  od:
+  fprintf(fd, "          if (");
+  for i from 1 to nc do
+    fprintf(fd, "test%d && ", init_cond[i][1]):
+  od:
+  fprintf(fd, "global_test)\n"):
+  fprintf(fd, "            break;\n"):
+  fprintf(fd, "        }\n\n"):
+
+
+  ######################################################
+  ############## Ziv' loop: testing final ##############
+  ######################################################
+
+  fprintf(fd, "      MPFR_LOG_MSG ((\"Truncation rank: %%lu\\n\", k));\n\n"):
+  fprintf(fd, "      err = ERRORANALYSISK + conditionNumber - MPFR_GET_EXP (s);\n\n"):
+  fprintf(fd, "      /* err is the number of bits lost due to the evaluation error */\n"):
+  fprintf(fd, "      /* wprec-(prec+1): number of bits lost due to the approximation error */\n"):
+  fprintf(fd, "      MPFR_LOG_MSG ((\"Roundoff error: %%Pu\\n\", err));\n"):
+  fprintf(fd, "      MPFR_LOG_MSG ((\"Approxim error: %%Pu\\n\", wprec-prec-1));\n\n"):
+  fprintf(fd, "      if (wprec < err+1)\n"):
+  fprintf(fd, "        correct_bits=0;\n"):
+  fprintf(fd, "      else\n"):
+  fprintf(fd, "        {\n"):
+  fprintf(fd, "          if (wprec < err+prec+1)\n"):
+  fprintf(fd, "            correct_bits =  wprec - err - 1;\n"):
+  fprintf(fd, "          else\n"):
+  fprintf(fd, "            correct_bits = prec;\n"):
+  fprintf(fd, "        }\n\n"):
+  fprintf(fd, "      if (MPFR_LIKELY (MPFR_CAN_ROUND (s, correct_bits, MPFR_PREC (y), rnd)))\n"):
+  fprintf(fd, "        break;\n\n"):
+
+  fprintf(fd, "      if (correct_bits == 0)\n"):
+  fprintf(fd, "        {\n"):
+  fprintf(fd, "          assumed_exponent *= 2;\n"):
+  fprintf(fd, "          MPFR_LOG_MSG ((\"Not a single bit correct (assumed_exponent=%%lu)\\n\",\n"):
+  fprintf(fd, "                         assumed_exponent));\n"):
+  fprintf(fd, "          wprec = prec + ERRORANALYSISK + conditionNumber + assumed_exponent;\n"):
+  fprintf(fd, "        }\n"):
+  fprintf(fd, "      else\n"):
+  fprintf(fd, "        {\n"):
+  fprintf(fd, "          if (correct_bits < prec)\n"):
+  fprintf(fd, "            { /* The precision was badly chosen */\n"):
+  fprintf(fd, "              MPFR_LOG_MSG ((\"Bad assumption on the exponent of %s(x)\", 0));\n", name):
+  fprintf(fd, "              MPFR_LOG_MSG ((\" (E=%%ld)\\n\", (long) MPFR_GET_EXP (s)));\n"):
+  fprintf(fd, "              wprec = prec + err + 1;\n"):
+  fprintf(fd, "            }\n"):
+  fprintf(fd, "          else\n"):
+  fprintf(fd, "            { /* We are really in a bad case of the TMD */\n"):
+  fprintf(fd, "              MPFR_ZIV_NEXT (loop, prec);\n\n"):
+
+  fprintf(fd, "              /* We update wprec */\n"):
+  fprintf(fd, "              /* We assume that K will not be multiplied by more than 4 */\n"):
+  fprintf(fd, "              wprec = prec + ERRORANALYSIS4K + conditionNumber\n"):
+  fprintf(fd, "                - MPFR_GET_EXP (s);\n"):
+  fprintf(fd, "            }\n"):
+  fprintf(fd, "        }\n\n"):
+
+  fprintf(fd, "    } /* End of ZIV loop */\n\n"):
+  fprintf(fd, "  MPFR_ZIV_FREE (loop);\n\n"):
+  fprintf(fd, "  r = mpfr_set (y, s, rnd);\n\n"):
+
+
+  ######################################################
+  ################ Clearing everything #################
+  ######################################################
+
+  fprintf(fd, "  mpfr_clear (s);\n"):
+  if (hardconstant=1)
+  then
+    fprintf(fd, "  mpfr_clear (x%d);\n", d):
+  fi:
+  if (type = FUNCTION_SERIES)
+  then
+    fprintf(fd, "  mpfr_clear (tmp);\n"):
+    for i from 1 to nops(exponents)-1 do
+      fprintf(fd, "  mpfr_clear (x%d);\n", exponents[i]):
+    od:
+  fi:
+
+  for i from 1 to nc do
+    fprintf(fd, "  mpfr_clear (tip%d);\n", init_cond[i][1]):
+  od:
+
+  fprintf(fd, "\n"):
+  fprintf(fd, "  MPFR_SAVE_EXPO_FREE (expo);\n"):
+  fprintf(fd, "  return mpfr_check_range (y, r, rnd);\n"):
+  fprintf(fd, "}\n"):
+  
+  fclose(fd):
+
+  for i from 1 to nops(required_mulsi) do
+    printf("You need to provide a mpfr_mul_si%d function.\n", required_mulsi[i]):
+    printf("  -> This can be achieved by a call to generate_muldivsin(\"mul\", %d):\n",  required_mulsi[i]):
+  od:
+  for i from 1 to nops(required_divsi) do
+    printf("You need to provide a mpfr_div_si%d function.\n", required_divsi[i]):
+    printf("  -> This can be achieved by a call to generate_muldivsin(\"div\", %d):\n",  required_divsi[i]):
+  od:
+
+
+  ######################################################
+  ################## Error analysis ####################
+  ######################################################
+  
+  printf("\n\n"):
+  printf("Before the loop, we have "):
+  for i from 1 to nc do
+    var := sprintf("tip%d", init_cond[i][1]):
+    printf("%s {%d}", var, find_in_error_counter(var, error_counter)):
+    if (i <> nc) then printf(", ") else printf("\n") fi:
+  od:
+  printf("Each step of the loop adds another {%d}\n", error_in_loop):
+  if (additional_error <> 0)
+  then printf("Moreover, the multiplication by %a adds another {%d} to each term before it is summed.\n", subs(n=k, s2/s1), additional_error)
+  fi:
+  printf("Finally, we have s = sum_(i=0)^(k-1) ( ti{%d + %dk} )\n", maxofti + additional_error + 1 - error_in_loop, error_in_loop):
+  printf("We bound it by {(k+%d)*2^(%d)}\n", ceil( (maxofti + additional_error + 1 - error_in_loop)/error_in_loop), ceil(log[2](error_in_loop))):
+
+  a := ceil( (maxofti + additional_error + 1 - error_in_loop)/error_in_loop):
+  b := ceil(log[2](error_in_loop)):
+
+  if (a > 0)
+    then var := sprintf("MPFR_INT_CEIL_LOG2 (prec + %d)", a)
+    elif (a=0) then  var := sprintf("MPFR_INT_CEIL_LOG2 (prec)")
+    else sprintf("MPFR_INT_CEIL_LOG2 (prec - %d)", -a)
+  fi:
+  if (b > 0) then var := sprintf("%s + %d", var, b+2) fi:
+  var := sprintf("sed -n -i 's/ERRORANALYSISPREC/%s/g;p' %s", var, filename):
+  system(var):
+
+  if (a > 0)
+    then var := sprintf("MPFR_INT_CEIL_LOG2 (k + %d)", a)
+    elif (a=0) then  var := sprintf("MPFR_INT_CEIL_LOG2 (k)")
+    else sprintf("MPFR_INT_CEIL_LOG2 (k - %d)", -a)
+  fi:
+  if (b > 0) then var := sprintf("%s + %d", var, b+2) fi:
+  var := sprintf("sed -n -i 's/ERRORANALYSISK/%s/g;p' %s", var, filename):
+  system(var):
+
+  if (a > 0)
+    then var := sprintf("MPFR_INT_CEIL_LOG2 (k + %d)", a)
+    elif (a=0) then  var := sprintf("MPFR_INT_CEIL_LOG2 (k)")
+    else sprintf("MPFR_INT_CEIL_LOG2 (k - %d)", -a)
+  fi:
+  var := sprintf("%s + %d", var, b+4):
+  var := sprintf("sed -n -i 's/ERRORANALYSIS4K/%s/g;p' %s", var, filename):
+  system(var):
+
+end proc:
diff --git a/tools/metaMPFR/metaMPFR_tests.mpl b/tools/metaMPFR/metaMPFR_tests.mpl
new file mode 100644
index 000000000..0a350804b
--- /dev/null
+++ b/tools/metaMPFR/metaMPFR_tests.mpl
@@ -0,0 +1,243 @@
+read ("metaMPFR_straightforwardAlgo.mpl"):
+
+f := AiryAi(x):
+deq := holexprtodiffeq(f, y(x)):
+rec := diffeqtorec(deq, y(x), a(n)):
+name_of_function := op(0,f):
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+f := erf(x):
+deq := holexprtodiffeq(f, y(x)):
+rec := diffeqtorec(deq, y(x), a(n)):
+name_of_function := op(0,f):
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = -(6*n+1)*(6*n+2)*(6*n+3)*(6*n+4)*(6*n+5)*(6*n+6)*a(n)/( (n+1)^3*(3*n+1)*(3*n+2)*(3*n+3)*12288000 ), a(0)=1 }:
+name_of_function := alpha:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+
+rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1 }:
+name_of_function := test0a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi }:
+name_of_function := test1a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1 }:
+name_of_function := test2a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi }:
+name_of_function := test3a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=2 }:
+name_of_function := test4a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi, a(1)=0}:
+name_of_function := test5a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1, a(1)=Pi }:
+name_of_function := test6a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(1)=0}:
+name_of_function := test7a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=0, a(2)=0 }:
+name_of_function := test8a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=0, a(1)=Pi, a(2)=2 }:
+name_of_function := test9a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=0, a(1)=0, a(2)=1 }:
+name_of_function := test10a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+7) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(4)=1, a(6)=2  }:
+name_of_function := test11a:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):
+
+
+rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1 }:
+name_of_function := test0b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi }:
+name_of_function := test1b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1 }:
+name_of_function := test2b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi }:
+name_of_function := test3b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=2 }:
+name_of_function := test4b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi, a(1)=0}:
+name_of_function := test5b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1, a(1)=Pi }:
+name_of_function := test6b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(1)=0}:
+name_of_function := test7b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=0, a(2)=0 }:
+name_of_function := test8b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=0, a(1)=Pi, a(2)=2 }:
+name_of_function := test9b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=0, a(1)=0, a(2)=1 }:
+name_of_function := test10b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+rec := { a(n+7) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(4)=1, a(6)=2  }:
+name_of_function := test11b:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):
+
+
+rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1 }:
+name_of_function := test0c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi }:
+name_of_function := test1c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1 }:
+name_of_function := test2c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi }:
+name_of_function := test3c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=2 }:
+name_of_function := test4c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi, a(1)=0}:
+name_of_function := test5c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1, a(1)=Pi }:
+name_of_function := test6c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(1)=0}:
+name_of_function := test7c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=0, a(2)=0 }:
+name_of_function := test8c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=0, a(1)=Pi, a(2)=2 }:
+name_of_function := test9c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+3) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=0, a(1)=0, a(2)=1 }:
+name_of_function := test10c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+
+rec := { a(n+7) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(4)=1, a(6)=2  }:
+name_of_function := test11c:
+name_of_file := sprintf("%a.c", name_of_function):
+printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
+generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):
+