Bjorn I won’t get the chance to try out your latest complex matrix version on bitbucket for a couple of weeks due to some trips coming up, but will update you when I do.
Thanks again,
Jonathan.
It is a bit of work to syntax update. The plugin should be fine though.
I have tested the fp plugin on my M1 mac with excellent results. I have tested about half of the math functions and they agree well with Mathematica.
However, I believe that the fpgamma(complex argument) function is not working properly and is returning incorrect results. If anyone else could check this, I would be grateful.
I imagine its the same with the old plugin ?
As in not new issue.
Thanks
Björn
(Note it seems to have 2 variants of fpGamma)
fpGamma(x As BigFloat) As BigFloat
fpGamma(x As BigComplex) As BigComplex
As far as I can tell, fpGamma(x As BigFloat) As BigFloat works while fpGamma(x As BigComplex) As BigComplex does not. And I suspect the old plug-in has the same issue.
1 Like
I did have an issue with very small valued numbers in matrix operations - had to pre-scale up the numbers and that solved the issue. Sorry for lack of detail, it was a while ago, but wondering if some of your math issues might be similar small-number ones.
I use a random number generator that scans from -5 to 5 and includes very small values. Only with Gamma() have I found any issues, but I have not tested the Matrix operations at all, only the basic math operations. My testing continues and I will make sure to test very small values and report back.
I have no knowledge to fix complex math. But I can dig up the relevant code path and post the code if someone Math savy is able to see something obvious in it.
I would be happy to take a look at the Gamma() code and see if I can find the errors. Gamma also fails for some negative real numbers too.
Declaration starts here:
{ (REALproc) fpGamma_BC, REALnoImplementation, "fpGamma(x As BigComplex) As BigComplex", REALconsoleSafe },
Which calls this one:
static REALobject fpGamma_BC(REALobject x)
{
REALobject z;
BigComplexData *zData, *xData;
xData = (BigComplexData *) REALGetClassData(x, &BigComplexClass);
if(xData==nil)
return x;
z = REALnewInstanceWithClass(globalBigComplexRef);
zData = (BigComplexData *) REALGetClassData(z, &BigComplexClass);
Gamma(zData->c, xData->c);
return z;
}/* fpGamma_BC */
Still no math logic there, but has redundant nil check which can never happen. And not efficient way of creating the result class. (Those are just because I have not refactored this part yet). Near bottom it calls the actual math function, taking in the data of the result class instance as first parameter and data of the current class as 2nd parameter.
That takes us here:
void Gamma(fpComplex& z, const fpComplex& x)
{
if(!x.im.i.n)
{
z.re = Gamma(x.re);
z.im = 0;
return;
}
if(x.re.i.n>=0)
{
GammaCalc(z, x);
return;
}
// now use reflection formula
fpComplex xt, zt;
xt = 1. - x;
GammaCalc(zt, xt);
z = Pi()/(zt*sin(Pi()*xt));
}/* Gamma */
The Gamma in the first iff there takes us here:
fp Gamma(const fp& x)
{
fp z;
Gamma(z, x);
return z;
}/* Gamma */
// uses Spouge formula for Gamma(x+1)
static void GammaCalc(fpComplex& z, const fpComplex& x)
{
fpComplex xt;
INT32 a, P, p;
p = getBlockPrec();
setBlockPrec(p+1);
xt = x - 1.; // so we get Gamma(x)
P = getDecPrec();
a = (P-log(P)/log(10))*659/526 - .5; // change log to Lg2
a = a+1;
//cout << "a = " << a << endl; // a = 37 for decPrec = 32
z = GammaSum(xt, a);
z = z*pow(xt+a, (xt+.5))*exp(-(xt+a));
setBlockPrec(p);
}/* GammaCalc */
static fpComplex GammaSum(const fpComplex& x, INT32 a)
{
fpComplex z;
GammaSum(z, x, a);
return z;
}/* GammaSum, */
// Spouge formula sum from Wikipedia
// Sqrt(2*Pi) + Sum(k=1 to N) c(k)/(x+k)
// where N = a - 1 and
// c(0) = sqrt(2*Pi)
// c(k) = (-1)^(k-1) (a-k)^(k-.5) e^(a-k) / (k-1)! k = 1. 2. ... , a-1
static void GammaSum(fpComplex& z, const fpComplex& x, INT32 a)
{
INT32 k, p;
fp ck, den, temp, myE;
p = getBlockPrec();
setBlockPrec(1.4*p);
exp(myE, 1.);
// c(1)
temp = a - 1;
den = exp(-temp);
z = sqrt(2*Pi());
z = z + sqrt(temp)/((x+1.)*den);
for(k=2;k<a;k++)
{
den = den*(k-1)*myE;
temp = a - k;
ck = pow(temp, k-.5)/den;
if(k==2*(k/2))
ck = -ck;
z = z + ck/(x+k);
}
setBlockPrec(p);
}/* GammaSum */
Sort of hard to see here. You can also pull the project from my Git page. And if you open in Visual Studio on Windows for example and go in fpPlugin.cpp, and search for
static REALobject fpGamma_BC(REALobject x)
Then thats where it begins…
You can then drill down as it calls subroutines by clicking them and pressing F12.
Thanks so much Bjorn! I will look into this. I wrote my own routing for Gamma so I should be able to find the issues.
1 Like
The code that you posted seems solid, but I have discovered that the pow() function is causing a serious problem. Indeed, fppow(complex, complex) invariably gives an incorrect result. I have written my own routine for gamma() based on the fpPlugin and it works well, However, the critical section which requires fppow(a,b) has to be replaced by
fpexp( b*fplog(a) ) if either a or b are complex (bigcomplex).
I am happy to send anyone interested a copy of my gamma routine based on fpPlugin. Just contact me via email.
On another matter, was it not the case the Xojo was going to allow us to make plugins using Xojo? Whatever happened to that idea?
Plugins in Xojo only presumably would mean “Packaging code” So classes and control like you can create now in Xojo except you could package it in compiled component. (Instead of open source Xojo class or Encrypted Xojo class).
Other limitations you would not go past, Xojo would still be Xojo, slower than C++ code, threading limitations etc.
Though there certainly would be good use cases for such like some sort of controls. That one would be fast to make in Xojo code. But complex algorithms or something requiring speed or would be very poor case for it.
1 Like
I am not sure based on your post what I should be changing. hmmm…
would then…
fpexp( b*fplog(a)
fix the pow function also ?
Bjorn,
Perhaps we should start with the pow function. It fails for all the complex tests, so anyone working with this plugin and complex numbers will likely benefit. If you could post the code I will work on it. That might take care of gamma() too!
This would be where the pow goes
Note that his z in first parameter seems is always the return value
// z = x^y
// we have a problem when xRe<0 and xIm going through zero
// this causes a discontinuous change in theta when we take the log
// so we should test xRe and adjust as necessary
// but this just moves the discontinuity elsewhere,
// so I do it this way to be consistent with my Extended Plugin
// 4/14/2010: I changed my mind.
void pow(fpComplex& z, const fpComplex& x, const fpComplex& y)
{
fpComplex temp;
fp myPi;
if(!y)
{
z = 1.; // so 0^0 = 1
return;
}
if(!x)
{
z = 0.; // actually undefined if y.re<0
return;
}
if(x==1)
{
z = 1.;
return;
}
myPi = Pi();
log(temp, x);
/*
if(x.re.i.n<0)
{
if(x.im.i.n<0)
temp.im = temp.im + 2*myPi;
}
*/
mul(temp, y, temp);
exp(z, temp);
}/* pow */
I have no idea what Bob is trying to do here, but it clearly does not work. I am not very comfortable in C, but I recommend that we simply replace this section with
z = exp( y*log(x) );
I’m not familiar with logs of complex arguments, but a regular log of real values is undefined for many argument values. It looks as if Bob’s pow function is checking for arguments that would lead to singularities, and only performing the equivalent z = exp( y*log(x) ) calculation if the singular cases don’t occur. So, it would be wise to ensure that without those checks the singular cases will still return correct values, and won’t throw an exception.
Well, looks like he did it, but because the parameters are type fpComplex, he made it in steps using a fpComplex temp variable:
temp = log(x)
temp = y * temp
z = exp(temp)
In the excluded part he looks for some special case when he finds negative values in components (one needs to read fpComplex to know) of both real and imaginary axis and adds 2𝜋 to the imaginary part of log(x) when he finds it.