Moshier's Cephes is a classical numerical library, here a sketch for computing
sin(x), http://www.netlib.org/cephes/doubldoc.html#sin

I do not comment much on range reduction (but integer multiples of Pi will need
increased precision for large values), but sin( Pi/2 - xi ) = cos(xi) allows to
reduce to the interval 0 .. Pi/4, if caring for both trigonometric functions.

Now to handle the Taylor series of sin(x) one can read it as x + x^3*series(x^2)
and E := t -> (sin(t^(1/2))-t^(1/2))/t^(3/2) is that analytic (!) function, that
means that sin(x) = x + x^3*E(x^2). Do not divide by x, but substract it.

Similar for cos(x) = 1 - x^2*someSeries(x^2).

Using Digits = 18 it turns out, that degree = 6 is enough for an approximation,
which reduces computational costs and possible errors through higher degrees:

  M:=numapprox:-chebyshev(E(t), t=0..evalf(Pi/4), 10^(-Digits+0)):
  M:=convert(eval(M, T=orthopoly[T]),horner):
  M:=unapply(evalf[18](M),t);

  Sin:= x -> x + x^3*M(x^2);

For visualization use acer's code (I would add style=point).

For higher precision find out that -Sum((-t)^(-1+k)/(2*k+1)!,k = 1 .. infinity)
is a series presentation for E(t), which is -1/6*hypergeom([1],[2, 5/2],-1/4*t).

As alternating series with decreasing coefficients it allows to estimate an 
cut off error after n term.

Since one actually wants x^3*series(x^2) it turns out that n=8 is enough for
usual double precision, beyond one can use Lambert's W function to estimate
it - or just stops in a brute way, if no changes can be found. Or best rely
on Maple to do it.



Edited (after acer's comment): the code for the plot, g = as f, 
but using Sin instead of S

 xRange:= 0..Pi/4;
 yRange:= -17..-15.5;
 plot([g,F],xRange,numpoints=201,
   view=[xRange,yRange],color=[pink,blue], 
   style=point, symbol=[circle, cross],
   myPlotDefault):
 plot(-16,xRange, color=grey):
 plots:-display(%%,%);

The relative error is ~ 1 DBL_EPSILON, may be better to use chebpade or minmax.

combine(GAL, radical) assuming x::real;

Almost always I use 15 Digits and often use rationals, and in most
cases I plot the situation:

restart;
Digits:=15;
-2.711505682*x-7.65*ln(3-x)-3/8*x^3+8.422482772;
E:=convert(%, rational); #plot(E, x=-6 .. 3);
plot(E, x=-1 .. 2);

This suggestes, there may be 3 real solutions.

RootOf(E,x); [allvalues(%)]; evalf(%);

tells me: there are 2 Reals and 1 complex at x~1/4

eps:=1e-6;
plot(E, x=1/4 - eps .. 1/4 + eps);

shows me, that x=1/4 is not a zero, it is ~ 2e-9.

But you must be aware, that your input is only 10 decimals.
So it may be that is just by the input being not exact enough:

op(E); eval(%, x=1/4); evalf(%); evalf[10](%); 

      -0.6778764205, -7.738746974, -0.005859375000, 8.422482772

In that constellation it is not unlikely that just a rounding error
may be the reason, that one does not see a zero, the result is
by adding the terms.

Edited: be aware, that your factor 7.65 at the log is only 3 decimals.
7.649999997 resp 7.65000000213486 should give a zero in x=1/4
for Digits = 10 resp Digits = 15. Try to check that first.

Remembered ther is a command in the LinearAlgebra package,
?LinearAlgebra[GenerateMatrix] may be what you are looking for

For me it seems difficult to see what you intend to achieve (especially from the
C source only), and why.

As suggested in your other thread I first would try to get the Maple script to be
ok, before calling it from something external.

And if I use some external, I carefully care for input and return types, your
function returns an integer, so I would write define_external that way (it is
too long ago, that I remember that wrapper things).

If I would need all the values of the various derivatives I would use (after the
prec-computations) just 1 function with 3 parameters u,v,w and an Array. That
is: why do you need the RTable?

And as already said: I would try to find a way, that the construction already works
in one of the systems (=Maple), before trying to interface with another (=C).


And for finding reasons for crashes I would try to reduce to the most simple case.

There are several things I do not understand in your C source (though it is some
time ago): MATRIX_OFFSET_C_RECT - what is it, have you checked indexing errors?
ALGEB DirectFastPointer - have you check correct type in calling? And so on ...

To locate errors you have to break it down, give the most simple situation with
just one function and most simple calls. Just to be sure, that there are not
multiple problems and that it worked before by lucky incidence.

Just my 2 cents.

But more or less: I do not understand what you want to do.

1. I would prefer to see the code. Not the cryptic thing in your post.

2. Some times ago I did such stuff (replicating Maple commands in C).
It is not so much fun. But it only makes sense (for me), if one has done
the work in Maple first (= have a good code working in Maple) and after
that go that painful route.

3. It may be worth to put all the Maple code into some proc, which
then is being called from C (if ever needed)

4. Maintaining such is a pain.

Why and for what do you you need it?

The appended sheet gives a solution through Maple without more theory:

Re-write as summing +-tan(d°), d = 1,3,5 ...89. Group as tan(d°) + tan(90°-d°),
convert to exp and simplify. Then Maple luckily finds the exact sum, it is 45.

MP_hirnyk_SumTan_ele.mws

The exact answer is 45. And here is a sketch, using (minimal) polynomials for
the sum of those algebraic numbers.

First note that Product( x - a(j), j=1 .. n) has coefficent -Sum(a(j), j=1..n)
in degree n-1 (and the product in degree 0), for example look at

  Product( x - a(j), j=1 .. 4); value(%): expand(%): collect(%, x);
  coeff(%,x, degree(%,x) - 1);

Now there is some luck with that task: the minimal polynomials for cos(...)
are all in x^2, if one takes the same arguments as for the tan(...).

That means, that one can derive the minimal polynomials for the tan(...) in a
simple way (without resolvents and doing the sin case), see Robert Israel's
contribution in http://www.mathkb.com/Uwe/Forum.aspx/math/77031/A-Question-about-Sums-of-Tangents

Addtionally it turns out that for these 45 terms only 6 different polynomials
actually occure they are of degree 1, 2, 4, 6, 8, 24. These degrees add up to
45 and we have 45 different terms as summands.

That means: the summands are exactly the zeros of those 6 polynomials.

Thus adding the coefficients in topdegree - 1 gives the sum (up to sign).

See attached sketchy file

MP_hirnyk_SumTan.mws

You really want 5 = 3*k = 3*1 = 3 ? Why ?

Note that your coefficients cover a wide range, from 1e-16 to 1e+18

Using 15 Digits and 'Isolate' a check for the result shwos, that they are large
and with 75 Digits one gets 1e-28 should be zero.

For 'Bivariate' it works, if one accepts that results are not checked:

BivariatePolynomial(sys , [lam1, lam2], verify=false);

This however returns undefines solution.

Guess, that are the reasons.

May be you want to dig for condition numbers (I hate that, but here
it seems to be needed)

PS: reallly, you should post things in a form ready for copy & paste.

You may look up ?polar (and the link polar coordinates ) and ?evalc

A 2-argument arctan is commonly used and your exp
would have an imaginary I = sqrt(-1)

Where the help pages should be more clear about that

1/myExpr = R*D^2 / ( R*D^2 + s*L + s^2*R*L*C );
isolate(%, myExpr);
expand(%);

                                           2
                                   s L    s  L C
                      myExpr = 1 + ---- + ------
                                      2      2
                                   R D      D

For the adjusted question (powering the Sum) I think it is

Limit( (1+2/n)^n, n=infinity) = exp(1)^2.

For that I used the asymptotic series in the already mentioned Temme paper
in a quite brute way to have hypergeom([1, n],[n+1],2) ~ (-1-2/(n+1)) and
then the Sum reduces to 1 + 2/n.

Ok, some reasoning is needed, why one can do the limit in 2 steps.

Since cos is invertible over that interval one can try:

Int(cos(a*x)*cos(x)^(a-2), x = 0 .. (1/2)*Pi);
Change(%, x=arccos(t),t);

                                  0

It seems that M15 does not find an anti-derivative here?

Numerical tests indicate that -1 is the limit. And an idea to prove it is as follows:

The Sum evaluates in terms of LerchPhi, to be convert to 2F1 and it can be seen,
that one has to show that hypergeom([1, n],[n+1],2) ~ -1.

The usual transformations for 2F1 do not apply due to the parameter constellation.

But having a look into papers of N Temme to get ideas beyond Maple abilities I took
"Large parameter cases of the Gauss hypergeometric function".

Now hypergeom([1, n],[n+1],z) ~ hypergeom([1, n],[n],z)  = (-1+z) for |z| < 1
by taking the summation version of 2F1-functions.

I have not followed his reasoning why this holds for the analytic extensions (and
into the branch cut).

Hopefully somebody finds a more simple way :-(

Edited: this is for the original question (without powering).