Fractalforums

Fractal Software => Programming => Topic started by: LionHeart on August 09, 2020, 11:11:18 AM

Title: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 09, 2020, 11:11:18 AM
I developed a C++ library for arbitrary precision arithmetic using MPFR. It seems to work very well until I try to use trig and log functions. These work well if I don't use multi-threading, but crash when I do. Small test programs show that the problems are not in the MPFR library, so I suspect I have a bug.

Does anyone know of any C++ library that allows mathematical operations on MPFR complex numbers?

If anyone wants to see my code, I am happy to share it. The full ManpWIN source code is available in downloads: https://fractalforums.org/index.php?action=downloads;sa=downfile&id=15 (https://fractalforums.org/index.php?action=downloads;sa=downfile&id=15)

Many thanks.

Title: Re: Is There a C++ Complex Library for MPFR?
Post by: Adam Majewski on August 09, 2020, 12:11:54 PM
http://www.multiprecision.org/mpc
http://arblib.org/
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: marcm200 on August 09, 2020, 12:18:59 PM
The kv library supports MPFR and transcendental functions - although I have only used it with double so far.

http://verifiedby.me/kv/index-e.html (http://verifiedby.me/kv/index-e.html)

Title: Re: Is There a C++ Complex Library for MPFR?
Post by: 3DickUlus on August 09, 2020, 12:53:31 PM
https://tspiteri.gitlab.io/gmp-mpfr-sys/mpc/Complex-Functions.html

built into GCC since version 4.5 requires GNU MPC version 0.8 or above, current (installed with my OS) version of gcc is 7.5.0 and mpc 3.1.0, so if you have the GCC compiler installed you already have multi precision complex types and functions.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 09, 2020, 01:08:25 PM
Thanks Adam Majewski, marcm200 and 3DickUlus.

I guess I didn't explain very well what I was asking. It's not the arithmetic functions, but the arithmetic operators.

These operators allow one to do arithmetic like:

int init_big_Formula_74()
    {
    if (!juliaflag)
zBig = qBig;
    zBig.x += param[0];
    zBig.y += param[1];
    return 0;
    }

int do_big_Formula_74()  // Flarium 106
    {
    double d;
    BigComplex z1, z3, zd, a;

    z1 = zBig;
    z3 = zBig.CCube();
    a = 1;
    zBig=((a-zBig-z3*zBig)/(zBig-(z3*4))-zBig)*qBig;  // sort of a good one
    zd = zBig-z1;
    d = zd.CSumSqr();
    return (d < MINSIZE);
    }

The problem to make a library for normal arithmetic operations on complex MPFR variables. Here is BigComplex class for the arithmetic operators that works well when I don't use multi-threading:

Code: [Select]
// BigComplex.h: interface for the Complex Bignum class.
//
//////////////////////////////////////////////////////////////////////

#pragma once

#include <math.h>
#include "BigDouble.h"
#include "Complex.h"

#define FALSE 0
#define TRUE 1
#define zerotol 1.e-50
#define BYTE unsigned char

class BigComplex
    {
    public:
BigComplex(void) { }
BigComplex(const BigDouble & real, const BigDouble & imaginary)
    {
    x = real;
    y = imaginary;
    }

BigComplex(const BigComplex & Cmplx1)// Copy Constructor
    {
    mpfr_set(x.x, Cmplx1.x.x, MPFR_RNDN);
    mpfr_set(y.x, Cmplx1.y.x, MPFR_RNDN);
    }

BigComplex(const BigDouble & value)
    {
    x = value;
    y = 0;
    }

~BigComplex(void);
inline BigComplex * operator =(const BigComplex &); // Assignment Operator
inline BigComplex * operator =(const BigDouble &); // Assignment to a big double Operator
BigComplex * operator =(const Complex &); // Assignment to a complex Operator
BigComplex * operator =(const double &); // Assignment to a double Operator
BigComplex * operator+=(const BigComplex &);
BigComplex * operator+=(BigDouble&);
BigComplex * operator-=(const BigComplex &);
BigComplex * operator-=(BigDouble &);
BigComplex * operator*=(const BigComplex &);
BigComplex * operator*=(BigDouble &);
BigComplex * operator++(void); // prefix ++ operator
bool    operator==(BigComplex &);
BigComplex operator^(BigDouble &);
BigComplex operator^(BigComplex &);
BigComplex operator^(BYTE &);
BigComplex operator +(const BigComplex &); // Addition Operator
BigComplex operator +(const BigDouble &); // complex add by double Operator
BigComplex operator -(const BigComplex &); // Subtraction Operator
BigComplex operator -(const BigDouble &); // complex subtract by double Operator
BigComplex operator -(void); // unary minus
BigComplex operator *(const BigComplex &); // Multiplication Operator
BigComplex operator *(const BigDouble &); // complex multiply by double Operator
BigComplex operator /(const BigComplex &); // Division Operator
BigComplex operator /(const BigDouble &); // complex divide by double Operator
BigComplex CSqr(void); // square
BigComplex CCube(void); // cube
BigComplex CInvert(void); // invert
double    CSumSqr(void); // real squared + imaginary squared
BigComplex CPolynomial(int); // take a complex number to an integer power
BigComplex CSin(void); // sine of a complex number
BigComplex CCos(void); // cosine
BigComplex CTan(); // tangent
BigComplex CSinh(); // hyperbolic sine of a complex number
BigComplex CCosh(); // hyperbolic cosine
BigComplex CTanh(); // hyperbolic tangent
BigComplex CExp(void); // exponent
BigComplex CDouble(void); // double   r = 2*n
BigComplex CHalf(void); // half     r = n/2
BigComplex CLog(void); // log
BigDouble  CFabs(void); // abs
BigComplex BigComplexPower(BigComplex &); // a^b
BigDouble x, y;
    };

Here is the BigDouble class:

Code: [Select]
// BigDouble.h: interface for the Bignum class.
//
//////////////////////////////////////////////////////////////////////

#pragma once

#include <math.h>
#include "mpfr.h"

#define FALSE 0
#define TRUE 1
#define zerotol 1.e-50    // 1.e-14
#define SAFETYMARGIN 5 // Multiply number of digits by a safety margin

extern int decimals;

class BigDouble
    {
    public:
BigDouble(void);
~BigDouble(void);
BigDouble(const double & value)
    {
    bitcount = decimals * SAFETYMARGIN;
    mpfr_init2(x, bitcount);
    mpfr_set_d(x, value, MPFR_RNDN);
    }

BigDouble(const BigDouble & Big) // Copy Constructor
    {
    bitcount = decimals * SAFETYMARGIN;
    mpfr_init2(x, bitcount);
    mpfr_set(x, Big.x, MPFR_RNDN);
    }

BigDouble * operator =(const BigDouble &); // Assignment Operator
BigDouble * operator =(const double &); // Assignment to a double Operator
BigDouble * operator+=(const BigDouble &);
BigDouble * operator+=(double&);
BigDouble * operator-=(const BigDouble &);
BigDouble * operator-=(double &);
BigDouble * operator*=(const BigDouble &);
BigDouble * operator*=(double &);
BigDouble BigSqr(void); // square is faster than * because of symmetry
BigDouble BigSqrt(void); // square root
BigDouble BigInvert(void); // invert
BigDouble BigAbs(void); // abs()
BigDouble BigX2(void); // double   r = 2*n
BigDouble BigHalf(void); // half     r - n/2
BigDouble BigCos(); // cos
BigDouble BigSin(); // sin
BigDouble BigTan(); // tan
BigDouble BigSinh(); // sinh
BigDouble BigCosh(); // cosh
BigDouble BigTanh(); // tanh
BigDouble BigLog(); // log
BigDouble BigLog10(); // log10
long   BigDoubleToInt(void); // convert a big double to int
double   BigDoubleToDouble(void); // convert a big double to double

int   ChangePrecision(int n); // increase precision to n digits -1 means precision is out of range

BigDouble operator +(const BigDouble &); // Addition Operator
BigDouble operator +(const double &); // add to double Operator
BigDouble operator -(const BigDouble &); // Subtraction Operator
BigDouble operator -(const double &); // subtract from double Operator
BigDouble operator -(void); // unary minus
BigDouble operator *(const BigDouble &); // Multiplication Operator
BigDouble operator *(const double &); // Multiply by double Operator
BigDouble operator /(const BigDouble &); // Division Operator
BigDouble operator /(const double &); // divide by double Operator

bool   operator==(BigDouble &);
bool   operator<(BigDouble &);
bool   operator>(BigDouble &);
bool   operator==(double &);
bool   operator<(double &);
bool   operator>(double &);

mpfr_t x; // pointer to BigNum
mpfr_prec_t bitcount;
    };


I would like to know if there is a library of these complex operations on MPFR. Does this make sense?

Thanks.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: marcm200 on August 09, 2020, 01:50:27 PM
   ~BigComplex(void);
   inline   BigComplex * operator =(const BigComplex &);   // Assignment Operator
Probably not related to your question, but I'm curious: Could ypu explain these a bit to me?

Is there a reason the destructor is not declared virtual?

You return a pointer of BigComplex in the assignment operators, I suppose for speed reasons - but I can't see a constructor taking a pointer. Is this implicitly converted to a reference in your compiler in code like a=b=c; ? Is there a reason not returning a reference? (I'm always looking for technical options to speed-optimized my code).

Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 09, 2020, 10:20:11 PM
Hi marcm200,

Thanks for trying to help me sort this mess out ;)

Part of my problem is that I am not a good C++ programmer :)  There is no reason not to make them virtual.

I think this is where my difficulty lies. I don't really have a constructor / destructor for the BigComplex class as these functions are performed by the BigDouble variables x and y. I would rather have these variables declared as mpfr_t than BigDouble. Then I would have to construct / destruct these variables.

Code: [Select]
// BigComplex.h: interface for the Complex Bignum class.
//
//////////////////////////////////////////////////////////////////////

#pragma once

#include <math.h>
#include "BigDouble.h"
#include "Complex.h"

#define FALSE 0
#define TRUE 1
#define zerotol 1.e-50
#define BYTE unsigned char

class BigComplex
    {
    public:
BigComplex(void);
~BigComplex(void);
BigComplex(const double & value)
    {
    bitcount = decimals * SAFETYMARGIN;
    mpfr_init2(x, bitcount);
    mpfr_set_d(x, value, MPFR_RNDN);
    mpfr_init2(y, bitcount);
    mpfr_set_d(y, value, MPFR_RNDN);
    }

BigComplex(const BigDouble & Big) // Copy Constructor
    {
    bitcount = decimals * SAFETYMARGIN;
    mpfr_init2(x, bitcount);
    mpfr_set(x, Big.x, MPFR_RNDN);
    mpfr_init2(y, bitcount);
    mpfr_set(y, Big.x, MPFR_RNDN);
    }


BigComplex operator =(const BigComplex &); // Assignment Operator
BigComplex operator =(const BigDouble &); // Assignment to a double Operator
BigComplex operator+=(const BigComplex &);
BigComplex operator+=(BigDouble&);
BigComplex operator-=(const BigComplex &);
BigComplex operator-=(BigDouble &);
BigComplex operator*=(const BigComplex &);
BigComplex operator*=(BigDouble &);
BigComplex operator++(void);
bool    operator==(BigComplex &);
BigComplex operator^(BigDouble &);
BigComplex operator^(BigComplex &);
BigComplex operator^(BYTE &);
BigComplex operator +(const BigComplex &); // Addition Operator
BigComplex operator +(const BigDouble &); // complex add by double Operator
BigComplex operator +(const double &); // complex add by double Operator
BigComplex operator -(const BigComplex &); // Subtraction Operator
BigComplex operator -(const BigDouble &); // complex subtract by double Operator
BigComplex operator -(const double &); // complex subtract by double Operator
BigComplex operator -(void); // unary minus
BigComplex operator *(const BigComplex &); // Multiplication Operator
BigComplex operator *(const BigDouble &); // complex multiply by double Operator
BigComplex operator *(const double &); // complex multiply by double Operator
BigComplex operator /(const BigComplex &); // Division Operator
BigComplex operator /(const BigDouble &); // complex divide by double Operator
BigComplex operator /(const double &); // complex divide by double Operator
BigComplex CSqr(void); // square
BigComplex CCube(void); // cube
BigComplex CInvert(void); // invert
double    CSumSqr(void); // real squared + imaginary squared
BigComplex CPolynomial(int); // take a complex number to an integer power
BigComplex CSin(void); // sine of a complex number
BigComplex CCos(void); // cosine
BigComplex CTan(); // tangent
BigComplex CSinh(); // hyperbolic sine of a complex number
BigComplex CCosh(); // hyperbolic cosine
BigComplex CTanh(); // hyperbolic tangent
BigComplex CExp(void); // exponent
BigComplex CDouble(void); // double   r = 2*n
BigComplex CHalf(void); // half     r = n/2
BigComplex CLog(void); // log
BigDouble  CFabs(void); // abs
BigDouble  CReal(void); // return real component
BigDouble  CImag(void); // return imaginary component
BigComplex BigComplexPower(BigComplex &); // a^b
mpfr_t x, y; // pointer to BigNum
mpfr_prec_t bitcount;
    };


However this fails when I try to address x or y. I hope all this makes sense.

Any hits are greatly appreciated.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: quaz0r on August 10, 2020, 01:52:56 AM
lionheart, unfortunately for you i think you may have reinvented the wheel, which is almost always the case when you dont look hard enough for a library and instead assume it doesnt exist.  check out mpfr c++

http://www.holoborodko.com/pavel/mpfr/
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 10, 2020, 02:02:55 AM
Hi quaz0r,

Thanks for replying.

This does not include complex numbers. There is no problem in my BigDouble library. The issues occur only when I use the BigComplex library.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: marcm200 on August 10, 2020, 09:00:49 AM
@LionHeart: The times I used multithreading (via OpenMP) I continuously made two mistakes: Assume a variable is local to a thread when I declared it outside and the thread messed with it, giving strange and unreproducible results in the parallel case while working perfectly single-threaded

And crashing almost always occured when this falsely assumed local variable did some memory allocation or deallocation. As your BigComplex constructor calls mpfr's functions (and I assume those allocate memory dnyamically) - are all those BigComplex declarations local to the threads? (Maybe your multithreading approach takes already care of this, then please ignore my answer).
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 10, 2020, 09:09:53 AM
Hi marcm200,

That's a good point. The code is complex (pun intended) and there may be a number of places where I might have issues. I'll start checking.

Thanks for the tip.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: sjhalayka on August 10, 2020, 06:29:39 PM
You know that there is the std::complex<T> type in C++, right?

https://github.com/sjhalayka/basic_fractals/blob/master/2d.cpp
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 10, 2020, 10:44:29 PM
Hi sjhalayka,

Yes thanks. That doesn't help with arbitrary precision however.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: sjhalayka on August 11, 2020, 10:37:12 PM
My bad.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 11, 2020, 10:45:27 PM
No worries, sjhalayka.

You took the time to answer  ;)
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: marcm200 on August 12, 2020, 09:26:41 AM
If you need a control to check your results (I used it with my custInt and custComplex struct for the subdivision algorithm:

Code: [Select]
typedef double kviabase; // change to MPFR arbitrary precision floating point type
typedef kv::interval<kviabase> kvia;
typedef kv::complex<kvia> kvcplxia;

// sample
kvcplxia cpi=kvcplxia(kv::constants<kvia>::pi(),kvia(0));
kvcplxia a=exp(2*cpi*kvcplxia(kvia(0),kvia(1))*cremer_number);

Afaik, this also should (I haven't tried) work with iabase being set to an MPFR datatype. But it is more fun to implement one's own class.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 12, 2020, 09:34:27 AM
Hi marcm200,

I'm sorry, I don't understand what you are doing here. Please clarify.

Many thanks.
Title: Re: Is There a C++ Complex Library for MPFR?
Post by: marcm200 on August 12, 2020, 09:53:03 AM
@LionHeart: The code posted is an example computing an interval enclosed value e^(2*i*pi*angle), currently with double precision as endpoints of the intervals (kviabase). If you set this to an MPFR type, you could compute a higher precision interval enclosure.

Computing the same formula with your BigDouble and BigComplex class at specific values or the angle, pi etc, you could then check if your result is near the interval (due to non-controlled rounding it might not be directly in).

I used this to compare my interval struct (double, complex and for small matrices)'s results to the ones computed with the kv library to see if my calculation path was correct. It was also a good means once I got stuck with my implementation to use the kv numbertype and continue writing the rest of the code (generating a subdivision image) - and then turn back to the struct code.

Title: Re: Is There a C++ Complex Library for MPFR?
Post by: LionHeart on August 12, 2020, 12:53:16 PM
Thanks for explaining, marcm200.

I understand now. Using typedef is a clever way to do it.

Thanks.