In that case - modulo an integer - I would use f:= x -> exp(-x) and then

  # shift to 0 ... newUpper and back again
  lower:= -1;
  F:= x -> f(modp(x-lower,2)+lower); # periodic extension by period=2

  myRange:= -6 .. 4; 
  plot(F, myRange, scaling = constrained);

Limitations are seen by (Digits=15)

  myRange:= 10^6 .. 10^6 + 10; 
  plot(F, myRange);
  #plot('F'(x), x=myRange); # use quotation marks for that!

Have you tried to use https://en.wikipedia.org/wiki/Fibonacci_number ? You would need to check the conventions for starting (those are not unique). They are partly based on closed formulae for Fib(n) and also give criteria to recognize naturals as Fibonacci numbers (have not checked them). And of course lots of references.

NB: For the first 2 your question is not "well posed", so just treat it as special case

You have a linear system in 3 variables and your task has only "trivial"
solutions, because the according matrix is invertible:

  with(LinearAlgebra):

  sys:= convert({E1, E2, E3}, rational);
  var:= [c, d, e];
  A,x:=GenerateMatrix(sys,var );

  Determinant(A); evalf(%);

I think that this does not make sense: you want to feed it to a C++, which is
to work in double precision, but your formulae are already too lengthy to feed
it (if they ever can be achieved).

For example: det_2 should not be computed by inverting J_2 and taking its det,
but taking det(J_2) and inverting. Which is the same, but much less expensive.

Already det(J_2) would fill ~ 5.000 pages and it is almost for sure, that this
neither can be feed to C++. Nor (after splitting) it will result in any value
which is reliable.

If you really want to go that way: to not use simplify/size (it will not give
you any advantage here, I tested it).

Perhaps some small advantage by codegen[makeproc] for the matrix and then
codegen[optimize] on that.

Now set up a proc, which computes the det of a 3x3 matrix. And feed it, having
numerical values. May be, some C++ libs already have that.

In this case you want to check the input(s) for fsolve.

For your Vector I get MM[1] = 660049668.737575*a[1](.1e-3) + ...
and you want to solve for a[1](.1e-3), ..., a[18](.1e-3)
What should that mean? You need variables.The last line has false curly brackets, the syntax is {seq(expr,k=1..18)}

It is easy to see, that the function has period = Pi and furthermore has a
point symmetry in Pi/2. So you can restrict to the interval 0 ... Pi/2 and
from that it is clear that there are at least 4 global extrema.

I have no explicit solution, but my guessing is: it is the first local extremum.

Here is brute (!) numerical approach.

f:=(cos(x)-cos((2*n+1)*x))/2/sin(x);
g:=eval(%, x= Pi/2*t); # so we search in t = 0 .. 1

nTst:=2000;
G:=eval(g, n=nTst);
2=numer(G);            # it is 2 at most 
fsolve(%, t=0, complex);
t0:=abs(Re(%));

diff(G,t)=0;
fsolve(%,t=t0);

#plot(G, t=0..nTst*t0);
plot(G, t=0..4*t0);

Edited 01.09.2013 to append a worksheet
http://www.mapleprimes.com/view.aspx?sf=151187/468374/MP_maximum_rational_.mws

Maple is a CAS - I do not understand why to expect a publication feature in top. It just makes it more complicated to be handle: some star instead a bar ... anyway: how about ` *` (holding a blank)?

Your sheet "crashes" for me: you only provide output ... but starring at it
I would aks: "solving for what variables?" Even if they would be linear the
result is parametric. And thus may be complicated.

Write the results into a variable (and only display parts which are needed) - nobody would read 1 Mio characters

?Int

Int( exp(x), x); value(%);

Find appended some coded suggestion. 

Using double precision I get ~ 40 msec for each of those integrals using NAG routines,
working on an usual office machine.

Not unusual - but one of the main problems is to get a reasonable (?!) result at all.


Decreasing required precision (either through epsilon = ... or digits = ...) in the
NAG integration routines does not change run times significantly. However that seems
to lead to errors only.

Shrinking the integration ranges to single precision also gives no advantages.

Expanding the range does seems not to avoid errors which may occure.

I played a bit in trying to make the code a bit more efficient (avoiding the symbolics
I used within it), but that would give only about 10 % in speed or so (but makes it
hard to change).

So only brushing it up for Compiler:-Compile and using that might have advantages.


Looking at some examples (in the enclosed sheet) it seems, that the problem is quite
fragile for inputs. Thus it may be, that a proper scaling or reformulation should be
done (for which one usually needs to know a 'bit' of the application theory and the
possible ranges of all the parameters and variables).

NumIntExample_work.mws
NumIntExample_work.pdf 

If you plot one of the integrals then you see, that it very peaked near zero (try some
interval +- 1e-3 and thus one can guess, that there is no need to integrate over the
full axis.

Indeed the integrand is of form exp(rational of degree (4,2))/sqrt(polynom of degree 4),
where the rational is ~ 2^15*z^2 + ...

Plotting that A:=rational around zero you see, that it like a quadric, but quickly exceeds
values, for which exp makes sense. With 15 Digits that is ~ +- 700.

Now fsolve(A = - 700, z) gives ~ 5*1e-3, beyond exp(A) will be zero and in theory will
dominate the denominator (one can confirm by checking your case of parameters, I think)

I have not look deeper, but that should be the way to look at the numerical problem
and then there should be no problem to use NAG in full precision, there is only a very
small range near zero where the task lives.

(and I have not look for the problem of increasing time)

Though I would not use isolate here:

(r^2*cos(theta)^2+r^2*sin(theta)^2)^2 = a^2*(r^2*cos(theta)^2-r^2*sin(theta)^2);
simplify(%);

                     4    2  2              2
                    r  = a  r  (2 cos(theta)  - 1)

isolate(%, r^2); 
                                     4
                      2             r
                     r  = ----------------------
                           2              2
                          a  (2 cos(theta)  - 1)

combine(%);
                                     4
                          2         r
                         r  = ---------------
                               2
                              a  cos(2 theta)

I was pondering (and that formulation of 'brute integration' certainly is the wrong way).

By symmetry one can reduce to the square over 0 .. Pi (and takes 4 times of that). I get

Int(2/(-x^2+1)^(1/2)*EllipticK(2*I*x*2^(1/2))-
  6*EllipticE(2*I*x*2^(1/2))/(-x^2+1)^(1/2)/(8*x^2+1),x = 0 .. 1);
         1
        /                     1/2                       1/2
       |   2 EllipticK(2 I x 2   )   6 EllipticE(2 I x 2   )
       |   ----------------------- - ----------------------- dx
       |           2     1/2            2     1/2     2
      /         (-x  + 1)            (-x  + 1)    (8 x  + 1)
        0

and after switching to hypergeomtrics and "lower b" for the 2nd term one has to show

1 = Int(1/(-x^2+1)^(1/2)*hypergeom([1/2, 1/2],[1],-8*x^2)-
        3/(-x^2+1)^(1/2)*hypergeom([1/2, 3/2],[1],-8*x^2),x = 0 .. 1);

         1
        /                                 2
       |   hypergeom([1/2, 1/2], [1], -8 x )
  1 =  |   ---------------------------------
       |                2     1/2
      /              (-x  + 1)
        0

                                            2
           3 hypergeom([1/2, 3/2], [1], -8 x )
         - ----------------------------------- dx
                         2     1/2
                      (-x  + 1)

Maple solves that integral - but in useless terms of MeijerG (only numerical).
BTW: at least the first term can be 'solved' using Prudnikov Vol 3 p. 153, 2.4.16.1

I have not done such for longer time. And all that stuff was / is (?) not documented
or explained nicely, it needs some trail-and-error.

Anyway: already in actual programs one would do some checks on data input before
processing - in your case you would check for numerical values.

Then one needs something like MapleToFloat - your first example is just an integer.

And there are commands like executing procedure and not a statement. May be a way
is to put all the critical things into some proc and just call that.


Overall it is somewhat painful work and a start to learn it are the examples which
come with Maple (have not looked, whether they are still delivered).


NB: I see that you have 64 bit ... just be aware that all fits together.