A math engine with support for complex numbers, calculus, and fractions

Question

I've been creating a math engine in C++. I have decided now is a good time to ask for feedback from the community. Right now the main "good" things it has is complex numbers, calculus and a fraction class. I am planning to add set theory sometime in the future, although I don't know how useful it would be to programmers. Since my project is for brushing up my C++ skills, but I'm making it as if I was making it as a general-purpose library.

For nsqrt and sqrt; the parameter precision represent the number of steps of a Newton-Raphson approximation to be performed

Is there anything I should improve on, change, add, etc.?

#include <iostream>
using namespace std;
namespace math{
const double mini = double (numeric_limits<double>::min());
const double e = 2.718281828459045;
const double pi = 3.141592653589793238462643383;
int def_prec = 16;
}

template <class t = double>
t exp(t base, int exponent){
if(exponent<0){
return 1/exp(base,-exponent);
} else if(exponent>0){
double output = 1;
for(int i=1;i<=exponent;i++){output *= base;}
return output;
} else {return 1;}
    
}

double nsqrt(double x, int precision = math::def_prec){
    double sum = 1;
    for (int i=1;i<=precision;i++){
        sum -= ((sum*sum)-x)/(2*sum);
    }
    return sum;
}

double esqrt(double num, double precision = math::def_prec){
    double h = num/precision;
    double y = 1;
    for (double x=1;x!=num;x+=h){
        y += h*(1/(2*y));
    }
    return y;
}

double sqrt(double num, int precision = math::def_prec){
    return nsqrt(num,precision);
}

double root(double num, int root, int precision = math::def_prec){
 double sum = 1;
  for (int i=1;i<=precision;i++){
      sum -= ((exp(sum,root))-num)/(root*exp(sum,root-1));
  }
  return sum;
}


class complex{
private:
    double r,i;
public:
    complex(double r_in,double i_in): r(r_in), i(i_in) {};
    double get_re() {return r;}
    double get_im(){return i;}
    double magnitude(int precision = math::def_prec){ return sqrt((r*r)+(i*i),precision);}
    complex operator + (complex a){
        r += a.get_re();
        i += a.get_im();
        return (*this);
    }
    complex operator + (double a){r += a;return (*this);}
    complex operator - (complex a){
        r -= a.get_re();
        i -= a.get_im();
        return (*this);
}
    complex operator - (double a){r -= a;return (*this);}
    complex operator * (complex a){
        r = (r*a.get_re())-(i*a.get_im());
        i = (r*a.get_im())+(i*a.get_re());
        return (*this);
    }
    complex operator * (double a){r *=a;i *= a;return (*this);}
    complex operator / (complex a){
           r = (r/a.get_re())-(i/a.get_im());
           i = (r/a.get_im())+(i/a.get_re());
           return (*this);
       }
    complex operator / (double a){r /=a;i /= a;return (*this);}
    void operator = (double a[2]){r=a[0];i=a[1];}
    void operator = (complex a){r=a.get_re();i=a.get_im();}
    bool operator == (complex a){return (r==a.get_re() && i==a.get_im());}
    bool operator == (double a[2]){return (r==a[0] && i==a[1]);}
    bool operator != (complex a){return (r!=a.get_re() || i!=a.get_im());}
    bool operator != (double a[2]){return (r!=a[0] || i!=a[1]);}
    void operator += (complex a){
        r += a.get_re();
        i += a.get_im();
    }
    bool operator > (complex a){return (r==a.get_re())? (i>a.get_im()) : (r>a.get_re());}
    bool operator < (complex a){return (r==a.get_re())? (i<a.get_im()) : (r<a.get_re());}
    bool operator <= (complex a){return (r==a.get_re())? (i<=a.get_im()) : (r<=a.get_re());}
    bool operator >= (complex a){return (r==a.get_re())? (i>=a.get_im()) : (r>=a.get_re());}
    void operator += (double a){r += a;}
    void operator -= (complex a){
        r -= a.get_re();
        i -= a.get_im();
    }
    void operator -= (double a){r -= a;}
    void operator *= (complex a){
        r = (r*a.get_re())-(i*a.get_im());
        i = (r*a.get_im())+(i*a.get_re());
    }
    void operator *= (double a){r *= a;i *= a;}
    void operator /= (complex a){
        r = (r/a.get_re())-(i/a.get_im());
        i = (r/a.get_im())+(i/a.get_re());
    }
    void operator /= (double a){r /= a;i /= a;}
    void operator ++ (){r++;}
    void operator -- (){r--;}
    
    ostream& operator<<(ostream& os)
    {
        os << r << "+" << i << "i";
        return os;
    }
    
    void operator = (int a){r=a;}
};

complex exp(complex base, int exponent){
if(exponent<0){
complex one(1,0);
return one/exp(base,-exponent);
} else if(exponent>0){
complex output(1,0);
for(int i=1;i<=exponent;i++){output *= base;}
return output;
} else {return base/base;}
    
}


template <class t = int>
t gcd(t X, t Y) {
    t pre = 0;
    if (X%Y == 0) {return (X>Y)? Y : X;}
    pre = X%Y;
    X = Y;
    Y = pre;
    while (X%Y != 0) {
        pre = X%Y;
        X = Y;
        Y = pre;
    }
    return pre;
}

int round(double x){return (x-int(x)>0.5)? int(x)+1 : int(x);}

class fraction{
private:
    int n,d;
public:
    fraction(int n_in,int d_in): n(n_in),d(d_in){};
    void simplify(){n/=gcd(n,d);d/=gcd(n,d);}
    double get_double(){return n/d;}
    int get_num(){return n;}
    int get_den(){return d;}
    fraction operator + (fraction a){
        n = (n*a.get_den())+(a.get_num()*d);
        d *= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator + (int a){n += a*d;return (*this);}
    fraction operator - (fraction a){
        n = (n*a.get_den())-(a.get_num()*d);
        d *= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator - (int a){n -= a*d;return (*this);}
    fraction operator * (fraction a){
        n *= a.get_num();
        d *= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator * (int a){n *= a;return (*this);}
    fraction operator / (fraction a){
        n /= a.get_num();
        d /= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator / (int a){n /= a;return (*this);}
    void operator += (fraction a){
        n = (n*a.get_den())+(a.get_num()*d);
        d *= a.get_den();
        simplify();
    }
    void operator += (int a){n += a*d;}
    void operator -= (fraction a){
        n = (n*a.get_den())-(a.get_num()*d);
        d *= a.get_den();
        simplify();
    }
    void operator -= (int a){n -= a*d;}
    void operator *= (fraction a){
        n *= a.get_num();
        d *= a.get_den();
        simplify();
    }
    void operator *= (int a){n*=a;}
    void operator /= (fraction a){
        n /= a.get_num();
        d /= a.get_den();
        simplify();
    }
    void operator /= (int a){n /= a;}
    void operator = (fraction a){n = a.get_num();d=a.get_den();}
    void operator = (int a){n=a;d=1;}
    
    bool operator < (fraction a){return (n*a.get_den()<a.get_num()*d);}
    bool operator > (fraction a){return (n*a.get_den()>a.get_num()*d);}
    bool operator >= (fraction a){return (n*a.get_den()>=a.get_num()*d);}
    bool operator <= (fraction a){return (n*a.get_den()<=a.get_num()*d);}
    bool operator == (fraction a){return (n*a.get_den()==a.get_num()*d);}
    bool operator != (fraction a){return (n*a.get_den()!=a.get_num()*d);}
    
    bool operator < (int a){return (n<a*d);}
    bool operator > (int a){return (n>a*d);}
    bool operator >= (int a){return (n>=a*d);}
    bool operator <= (int a){return (n<=a*d);}
    bool operator == (int a){return (n==a*d);}
    bool operator != (int a){return (n!=a*d);}
    operator double () {return n/d;}
};

fraction tofrac(double x, int precision = 1000000000){
    long gcd_ = gcd(round(x-int(x) * precision), precision);
    fraction output(round(x-int(x)*precision)/gcd_,precision/gcd_);
    return output;
}

template <class t = double>
t limit(double (*f)(t),t to) {
    return (f(to+math::mini)+f(to-math::mini))/2;
}

template <class t = double>
t derivative(double (*f)(t),t x){
    return ((f(x+math::mini)-f(x))/math::mini+(f(x-math::mini)-f(x))/(-math::mini))/2;
}

template <class t = double>
t sigma(double (*f)(t),int from, int to){
    t out = 0;
    for (int n=from;n<=to;n++){out += f(n);}
    return out;
}

template <class t = double>
t pi(double (*f)(t),int from, int to){
    t out = 1;
    for (int n=from;n<=to;n++){out *= f(n);}
    return out;
}

double factorial1(int x) {
    double out = 1;
    for(int n=x;n>0;n--){out*=n;}
    return out;
}

template <class t = double>
t sin(t x, int precision = math::def_prec){
  t output = 0;
  for(int n=0;n<=precision;n++){
    output += exp(-1,n)*(exp(x,(2*n)+1)/factorial1((2*n)+1));
  }
    return output;
}

template <class t = double>
t cos(t x, int precision = math::def_prec){
    t output = 0;
    for(int n=0;n<=precision;n++){
      output += exp(-1,n)*(exp(x,2*n)/factorial1(2*n));
    }
      return output;
}

template <class t = double>
t tan(t x, int precision = math::def_prec){return sin(x,precision)/cos(x,precision);}

template <class t = double>
t ln(t x, int precision = math::def_prec){
    t output = 0;
    for(int n=1;n<=precision;n++){
        output += (exp(x-2,n)*exp(-1,n-1))/n;
    }
      return output;
}

complex ln(complex x, int precision = math::def_prec){
    complex output(0,0);
    for(int n=1;n<=precision;n++){
        output += (exp(x-1,n)*exp(-1,n-1))/n;
    }
      return output;
}

template <class t = double>
t log(t base, t x, int precision = math::def_prec){return ln(x,precision)/ln(base,precision);}

template <class t = double>
t derivative_fast(t (*f)(t),t x){return (f(x+1)-f(x-1))/2;}

double pi(int precision = math::def_prec){
    double output = 0;
    for(int i=0;i<=precision+(precision%2);i++){
        output += exp(-1.0,i)/((2*i)+1);
    }
    return output*4;
}

//this pi algorithm isn't my work
double spigotpi(long digits){
    
    long LEN = (digits / 4 + 1) * 14;  //nec. array length
    
    long array[LEN];                   //array of 4-digit-decimals
    long b;                        //nominator prev. base
    long c = LEN;                  //index
    long d;                        //accumulator and carry
    long e = 0;                    //save previous 4 digits
    long f = 10000;                //new base, 4 decimal digits
    long g;                        //denom previous base
    long h = 0;                    //init switch
    
        for (; (b = c -= 14) > 0;) {    //outer loop: 4 digits per loop
            for (; --b > 0;) {      //inner loop: radix conv
                d *= b;            //acc *= nom prev base
                if (h == 0)
                    d += 2000 * f;   //first outer loop
                else
                    d += array[b] * f;   //non-first outer loop
                g = b + b - 1;           //denom prev base
                array[b] = d % g;
                d /= g;            //save carry
            }
            h =  (("%d"), e + d / f);//print prev 4 digits
            std::cout << h;
            d = e = d % f;         //save currently 4 digits
                                   //assure a small enough d
        }
        if (getchar())
        {
            return 0;
        }
        return 0;

    
}
double pi2(long precision){
    double output = 1;
    for (long k = precision;k>0;k--){
        output *= 1.0+(k/((2.0*k)+1.0));
        cout << output << endl;
    }
    return output;
}

double e(int precision = math::def_prec){
    double output = 0;
    for(int i=0;i<=precision;i++){
        output += 1/factorial1(i);
    }
    return output;
}

double abs(complex x, int precision){return x.magnitude(precision);}

template <class t = double>
double abs(t x){return (x<0)? -x : x;}


complex sqrt(complex x, int precision = math::def_prec){
    if (x.get_im()!=0){
    complex output((1/sqrt(2,precision))*sqrt(x.magnitude()+x.get_re(),precision),((x.get_im()/abs(x.get_im()))/sqrt(2,precision))*sqrt(x.magnitude()-x.get_re(),precision));
    return output;
    } else {
        complex output(sqrt(x.get_re()),0);
        return output;
    }
}

template <class t = double>
t quadratic(t a,t b, t c,t& otherroot, int precision = math::def_prec){
    otherroot = (-b+sqrt((b*b)-4*a*c,precision))/(2*a);
    return (-b-sqrt((b*b)-4*a*c,precision))/(2*a);
}

complex e_to_xi(double x, int precision = math::def_prec){
   complex output(cos(x,precision),sin(x,precision));
    return output;
}

template <class t>
t e_to_x(t x, int precision = math::def_prec){
    t output = 0;
    for(int i=0;i<=precision;i++){
        output += exp(x,i)/factorial1(i);
    }
    return output;
}

complex e_to_x(complex x, int precision = math::def_prec){return e_to_xi(x.get_im(),precision)*e_to_x(x.get_re(),precision);}

complex exp(complex base, complex exponent, int precision = math::def_prec){return e_to_x(ln(base)*exponent);}

template <class t = double>
t exp(t base, double exponent, int precision = math::def_prec){
    fraction asfrac(0,0);
    asfrac = tofrac(exponent);
    return exp(root(base,asfrac.get_den(),precision),asfrac.get_num());
}

template <class t = double>
t zeta(t s, int precision = math::def_prec){
    t output(0,0);
    for(int n=1.0;n<=precision;n++){output += 1/exp(n,s);}
    return output;
}

template <class t = double>
t eta(t s, int precision = math::def_prec){
    t output = 0;
    for(int n=1;n<=precision;n++){output += exp(-1,(n-1)%2)/exp(n,s);}
    return output;
}


template <class c=double>
c rsum(double (*f)(c), int n,c x[n],c t[n-1]){
    c output = 0;
    for(int i=0;i<=n-1;i++){
        output += f(t[i])*(x[i+1]-x[i]);
    }
    return output;
}

template <class t=double>
t inf(t a, t b){return (a<b)? a:b;}

template <class t=double>
t sup(t a, t b){return (a>b)? a:b;}

template <class c=double>
c ldsum(double (*f)(c), int n,c x[n]){
    c output = 0;
    for(int i=0;i<=n-1;i++){
        output += f(inf(x[i],x[i+1]))*(x[i+1]-x[i]);
    }
    return output;
}

template <class c=double>
c udsum(double (*f)(c), int n,c x[n]){
    c output = 0;
    for(int i=0;i<=n-1;i++){
        output += f(sup(x[i],x[i+1]))*(x[i+1]-x[i]);
    }
    return output;
}
template <class c=double>
c rint(c a, c b, c (*f)(c), int precision=math::def_prec){
    c partition[precision];
    for(int i=0;i<=precision;i++){
        partition[i] = a+(((b-a)/precision)*(i+1));
    }
    return rsum(f, precision, partition, partition); //use the lower bounds as tags
}

template <class c=double>
c udint(c a, c b, c (*f)(c), int precision = math::def_prec){
    c P[precision];
}

Answer 1

using namespace std

Given that you're using template and declare your class, you're working in a header-only library. However, you should never use using namespace in header. Not only is it bad practice, but in this case you might conflict with std::complex.

Surprising operator semantics

Also, almost all operators don't follow the usual guidelines and are therefore either misleading or even plain wrong. Let's have a look at operator-():

    complex operator - (complex a){
        r -= a.get_re();
        i -= a.get_im();
        return (*this);
}

Let's stop. This operator changes the value of the current item. That's completely unexpected. If I have

complex a = foo();
complex b = bar();
complex a_cp = a;
complex b_cp = b;
complex c = a - b;

then I expect a == a_cp and b == b_cp. However, this isn't true. This also prevents us from using operator-() on constant values.

Code duplication

Next, we see a lot of duplication. Let's stay at that operator and look at its compound assignment variant:

    complex operator - (double a){r -= a;return (*this);}

    void operator -= (complex a){
        r -= a.get_re();
        i -= a.get_im();
    }
    void operator -= (double a){r -= a;}

We can immediately see that there is duplication between - and -=. It's almost the same code. Usually all compound assignment operations return T& instead of void or T, so let's fix that first:

    complex& operator -= (complex a){
        r -= a.get_re();
        i -= a.get_im();
        return *this;
    }

Scott Meyers: prefer non-member, non-friend functions

Now operator- can be rewritten using operator-=, which is also an item in Scott Meyers' book:

complex operator-(const complex &a, const complex &b) {
    complex tmp = a;
    tmp -= b;
    return tmp;
}

Note that the operator above is free-standing. It's not a member function. This allows us to use it with types that can be converted into complex.

Provide conversion constructors

The code above is missing a double or even double[2] variant. However, that's not a problem if we have a fitting constructor at hand, e.g.

class complex {
...
public:
     complex(double r)    : r(r), i(0) {}
     complex(double c[2]) : r(c[0]), i(c[1]) {}

     complex(const complex&) = default;
     complex(complex&&)      = default;
     complex& operator=(const complex&) = default;
     complex& operator=(complex&&) = default;
}

While we're at it, we should also add the usual copy and move constructors resp. assignments.

Use const& where possible

Almost all functions in complex copy their argument. That's not necessary for them.

Use const for functions that may not change the class internals

The get_re() and get_im() functions don't change the complex among others, so mark them as const:

class complex {
...
    double get_re() const {
        return r;
    }
    double get_im() const {
        return i;
    }
...

Use a code formatter

The code formatting and style changes all too much. Stick to a readable style and enforce it, either manually (which takes time) or with a tool like clang-format. Your IDE/editor might have a built-in feature.

Answer 2

Note: This isn't comprehensive, I'm adding to the previous reviews as well.

Numeric Constants

Take a look at numeric constants in C++. There's no need for you to define them yourself.

Dangerous Iteration

    double h = num/precision;
    double y = 1;
    for (double x=1;x!=num;x+=h){
        y += h*(1/(2*y));
    }

It's not guaranteed that x will ever reach the value of num. Here, you could use the iteration style for (double x=1; x<=num; x+=h). However, this seems like the nth approximation to the sqare root, so why not use an integer counter? After all, you're not using x at all inside the loop!

Unusual Iteration

  for (int i=1;i<=precision;i++){

The idiomatic way to write that is to start with zero and then to compare for equality:

  for (int i=0; i!=precision; ++i) {

Unit Tests

This isn't strictly part of reviewing your code. Still, for libraries like this, using Test-Driven Development (TDD) is probably the standard. I believe that some of the faults in your code would have become obvious.

Documentation

What are nsqrt(), esqrt() and sqrt()? Why add another sqrt() overload at all? Why not use the one provided by the standard library? These three are just examples of things you need to put into a suitable context so that others can use your code.

In addition, putting down a few of your thoughts will make the code easier to maintain for your "now +6months" future self. Also, sometimes you will find that an initial thought is actually a stupid idea when you explain it to someone else. The same thing happens if you write it down as comment (which should be for someone else as well).

Answer 3

As a complement to the other answers, here are a few brief comments on the math part.

Disclaimer: I'm no numerical analyst, so there are probably more substantial efficiencies than the ones I suggest. I would say that efficiencies are important, since usually math code will be the backbone of a numerical application.

First a few general principles:

• As a general principle, you want to try and minimize the number of operations you make, to minimize errors.

• Errors do exist. A compiled C program stores its numbers in binary (or hex if you prefer). This has some unexpected consequences. A typical example is that the number 1/10 (one over ten) is periodic in base 2. The consequence of this is that it only stored as an approximation, and that any computation involving this number will carry errors over. This happens left and right with any usual floating point implementation. Very often it does not generate problems; but in a library you want to minimize computations, as mentioned above.

• Any function implemented in your library is likely to be used repeatedly by a program. This makes efficiency fairly important.

Now to a (very) non-exhaustive list of concrete manifestations of the above:

• You exp function simply multiplies the base, exponent times. This can be streamlined a little (timewise) by keeping products, and this can be achieved by using the binary representation of the exponent. For instance, if you want to do exp(base, 20), you notice that 20 is 2^4 + 2^2. So you can do (base^2)^4 + (base^2)^2. This makes you calculate three products and one sum instead of 20 products.

• You definitely want to do complex multiplication in polar form, where the products turn to sums.

• The limit and the derivative are calculated as (f(x+h)+f(x-h))/2 and ((f(x+h)-f(x))/h+(f(x-h)-f(x))/(-h))/2. This is bad numerically. In numerical calculuations, you want to avoid very small denominators, as this augments your errors greatly.

• In the sine and the cosine, you are re-calculating the factorial for each term of the series, instead of calculating it as you go: in the term n+1, the divisor is d_{n+1} = d_n * (n+1). Same with the power of X and the sign.

• The tangent function would be better served by using its own Taylor series than by calculating both the sine's and cosine's and dividing. In any case, you need to use reduction to the first quadrant and even to the first half of the first quadrant.

Answer 4

In addition to the other answer,

1. Put all your code inside a namespace. Why are only constants inside a math namespace? Also, math might be too common a name to use as namespace.

A good idea would something like this:

namespace my_math_library
{
    namespace constants
    {
       // put constants here
    }

    // rest of the stuff here
} 

---

2. You can declare all constants, functions and classes as constexpr.

---

3. template<class t = double>. You don't need to actually provide a default template argument, since the template argument will just get deduced from the function argument.

---

4. Why are some of your functions templated and others are not? (nsqrt, esqrt)

---

5. for(double x = 1; x != num; x += h). Comparing two floating point numbers using == or != is a recipe for disaster. Compare it against a tolerance (epsilon) e.g. std::numeric_limits<double>::epsilon()

---

6. x-int(x) * precision. Be careful of operator precedence. int(x) * precision is evaluated first. For any value of x >= ~4, you are invoking undefined behavior since the value doesn't fit inside an int.

---

7. t derivate(double (*f)(t), t x). Instead of passing a function pointer, pass a templated argument:

template<typename Func, typename T>
T derivate(Func&& func, T x)

That way, you can also pass in lambdas and functors.

---

8. rsum(double (*f)(c), int n,c x[n],c t[n-1] is not valid. Also, use consistent notation when using templates.

---

template <class c=double>
c udint(c a, c b, c (*f)(c), int precision = math::def_prec){
    c P[precision];
}

Variable length arrays are illegal in C++. Also, function isn't doing anything.

Answer 5

Coming from a math background and with no C++ experience, here's what I notice:

You jump between if (condition) { body inline } and

if (condition) {
    body by itself
}

Stick with the second.

Your function and parameter names need to be more descriptive. Failing that, you should comment why you're doing what you're doing.

precision should be iterations. Even better would be to leave it as precision (or tolerance), and calculate enough iterations until you are within precision of the answer (or at least your current answer changes by at most precision).

Generally, exp means e^x. It would be better to name yours power.

What is the difference between nsqrt and esqrt? I see that nsqrt is the Newton-Rahpsom method. I don't understand esqrt. Others have remarked about the problem of x!=num. I'll point out that (1) you should have started with x=0, and (2) this is a convoluted way to have precision iterations, so you would have wanted for (int i=1;i<=precision;i++) and dispensed with h.

I'm pretty sure that you can remove the if portion of gcd and it will still work.

Shouldn't class frac have some sort of method to return a string representation?

tofrac should make use of simplify. Others have noticed that it's broken, eg, anything from 0 to .5 will return 0. I think you wanted numerator=round(x*precision); return fraction(numerator, precision).simplify();. Now precision doesn't mean iterations, it means "my guess for a denominator". But this approach will not work well for denominators that don't divide precision. It would be better to return the fraction with the smallest denominator that is within precision of x. This can be done in O(denominator) time.

limit should document its limitations. It will work ok for a removable discontinuity. It won't work for other types of limits.

sigma should be named sum or summation. pi should definitely be named prod or product.

Several of your Taylor polynomial loops use exp(-1,n) to calculate a sign. It will be faster to just have sign *= -1 in the loop, although that won't be too big of a deal if you're only running a loop 16 times (the default).

You have a constant pi and a function pi that returns an approximation of the constant (in addition to the pi that should be prod). Why do you have a function to approximate a constant you already have stored? Similarly, you have a function e to approximate your constant e.

What is spigotpi supposed to do? Print pi? Or something else? A math engine shouldn't be printing things; that should be the calling program's decision.

What is pi2? Again with printing.

complex sqrt should also return a complex number if x.get_re()<0 (actually, does this even get called in that case? It should). You need an intermediate value somewhere in there; I can't follow what you're doing. And I'm suspicious that you have sqrt(2).

inf and sup should be min and max. inf and sup would only make sense if you can have an infinite array in the input, and I don't know of a programming language that could handle that. (Or maybe it could return a lower bound within tolerance of the true inf, but you'd need a lot more code to pull that off.)

In rsum, you're assuming that x[i]<=t[i]<=x[i+1]. You should document that, and possibly assert it as well. If rsum stands for Riemann sum, then ldsum must be lower Reimann sum and udsum must be upper Riemann sum. What's the d stand for?

Does rint mean "Riemann sum with equally spaced partitions"? You should document that. In your comment, I would say "left hand endpoint" instead of "lower bound" (I would assume "lower bound" to be ldsum).

Off topic:

Series with denominators that are O(n) converge horribly slowly (in ln and pi). You should find better algorithms. zeta and eta will also converge slowly, but there's not much you can do with those.