Debugger crashing?

Targets macOS
1

Can anybody please explain to me why my code crashes the debugger? I use XOJO 2023 R 1.1 with Bob Delaney’s ComplexMatrix Plugin on an iMac M1. Björn kindly ported it to be usable with Apple Silicon and offers it on his site einhugur.com as an open source project.

Var bet,re,im As Double
Var z,x,y As Integer
Var c1,c2 As Complex
Var tabel(300,0) As Double
Var faktor As Double = 0.00425
Var XrEcke As Double = -1.5
Var YiEcke As Double = -1.7

c1= New Complex(1.0,0.1)

For y=0 To 800
  For x=0 To 800
    re=x*faktor+XrEcke
    im=y*faktor+YiEcke
    c2 = New Complex(re,im)
    z=0
    Do
      c2=c2.Sinh.Pow(4)'<< crashes here when z=2
      c2=Abs(c2.Real)+Abs(c2.Imag)+c1
      bet=c2.Abs
      z=z+1
    Loop Until z>300 Or bet > 999
    tabel(x,y)=bet\6000
    If x=277 And y=0 Then Break
  Next
Next

If I run it and, after it broke into the debugger, step a few steps further, the debugger crashes. I think the culprit is the sinh. But why? The debugger shows no weird values until the moment of the crash. Maybe sinh has no proper exception handling.

Any hint appreciated

2

Perhaps put a Try/Catch round the sinh to pin down values etc? Just for testing.

3

Possible bug in the plugin. Maybe see if Björn can help?

4

Break the operations in more steps like

c3 = c2.Sinh
c2 = c3.Pow(4)

Observe the values and the new behavior.

2 Likes
5

Its a low level crash so debugger or try & catch has no chance.

From my tests then it was the Sinh that crashes.

I do not think we will see this one fixed unless someone with some sort of understanding of the actual math going on does it.

The Sinh in this case jumps to multiple routines internally.

I basically got no knowledge of complex math, so I can only fix obvious memory clobbering and obvious programming flaws in Bob’s old plugins.

Ideally the crash should be isolated better somehow though since it is very hard to work with it in such loop.

2 Likes
6

sinh basically comes down to a couple of exponents:

sinh θ = (e^θ − e^(−θ)) / 2

Numeracy, Maths and Statistics - Academic Skills Kit

Anything else it is doing is presumably to do with the matrix elements.

7

Relevant code:

This first one is where it starts that one only wraps his inner routine:

static REALobject SinhCM(REALobject instance)
{
	ClassData(ComplexMatrixClass, instance, ComplexMatrixData, data);

	if (data->nr < 1 || data->nc < 1 || data->nr != data->nc)
	{
		ThrowRowColumnMismatchException();
		REALLockObject(instance);
		return instance;
	}

	REALobject z = REALnewInstanceOfClass(&ComplexMatrixClass);
	ClassData(ComplexMatrixClass, z, ComplexMatrixData, zData);

	zData->nr = data->nr;
	zData->nc = data->nc;
	zData->array = (Complex**)malloc(zData->nr * sizeof(Complex*));
	for (int32_t i = 0; i < zData->nr; i++)
		zData->array[i] = (Complex*)malloc(zData->nc * sizeof(Complex));

	mySinhCM(zData->array, data->array, zData->nr);

	return z;

}/* SinhCM */

Then it goes to his inner routine:

void mySinhCM(Complex **z, Complex **x, int n)
{
	Complex		**temp;
	long			i, j;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		z[0][0].Re = sinh(x[0][0].Re)*cos(x[0][0].Im);
		z[0][0].Im = cosh(x[0][0].Re)*sin(x[0][0].Im);
		return;
	}
	
	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	EquateCM_CM(temp, x, n); // temp = x
	NegateCM(temp, n); // temp = -x
	myExpCM(z, temp, n); // z = exp(-x)
	EquateCM_CM(temp, z, n); // temp = exp(-x)
	
	myExpCM(z, x, n); // z = exp(x)
	SubCM_CM(z, z, temp, n); // z = exp(x) - exp(-x)
	DivCM_D(z, 2, n); // z = sinh(x)
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
		
	for(i=0;i<n;i++)
		free(temp[i]);
	
	free(temp);
	
}/* mySinhCM */

And as you can see then his inner routine calls several others…

static inline void EquateCM_CM(Complex **z, Complex **x, int n)
{
	long	i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			z[i][j] = x[i][j];

}/* EquateCM_CM */

static inline void NegateCM(Complex **z, int n)
{
	long	i, j;
		
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			NegateCC(z[i][j]);

}/* NegateCM */

oid myExpCM(Complex **z, Complex **x, int n)
{
	long			i, j, k;
	Complex			**temp, **currentTerm, **zCurrent;
	double			temp1;
	double			maxEl;
	int             m;
	Complex			**xt;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		temp1 = exp(x[0][0].Re);
		z[0][0].Re = temp1*cos(x[0][0].Im);
		z[0][0].Im = temp1*sin(x[0][0].Im);
		return;
	}
	
	// find maximum element of x in absolute value
	maxEl = 0;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			if(AbsC(x[i][j])>maxEl)
				maxEl = AbsC(x[i][j]);
	
	if(maxEl==0)
	{
		// z = unit matrix
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
				if(i==j)
				{
					z[i][j].Re = 1.;
					z[i][j].Im = 0;
				}
				else
				{
					z[i][j].Re = 0;
					z[i][j].Im = 0;
				}
		return;
	}
	
	xt = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		xt[i] = (Complex *)malloc(n*sizeof(Complex));
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			xt[i][j].Re = x[i][j].Re;
			xt[i][j].Im = x[i][j].Im;
		}
	
	m = 1;
	if(maxEl<(double)0x7FFFFFFF)
	{
		m = (int)maxEl;
		m++;
		if(m>1)
		{
			// divide matrix xt by m
			for(i=0;i<n;i++)
				for(j=0;j<n;j++)
				{
					xt[i][j].Re = xt[i][j].Re/m;
					xt[i][j].Im = xt[i][j].Im/m;
				}
		}
	}

	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	currentTerm = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		currentTerm[i] = (Complex *)malloc(n*sizeof(Complex));
		
	zCurrent = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		zCurrent[i] = (Complex *)malloc(n*sizeof(Complex));
	
	
	
	// make z, currentTerm, and zCurrent Unit matrices
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			if(i==j)
			{
				z[i][j].Re = currentTerm[i][j].Re = zCurrent[i][j].Re = 1;
				z[i][j].Im = currentTerm[i][j].Im = zCurrent[i][j].Im = 0;
			}
			else
			{
				z[i][j].Re = currentTerm[i][j].Re = zCurrent[i][j].Re = 0;
				z[i][j].Im = currentTerm[i][j].Im = zCurrent[i][j].Im = 0;
			}
	
	
	for(k=1; ;k++)
	{
		MulCM_CM(temp, xt, currentTerm, n);
		EquateCM_CM(currentTerm, temp, n);
		DivCM_D(currentTerm, k, n);
		
		AddCM_CM(z, z, currentTerm, n);
		
		if(isEqualCM_CM(z, zCurrent, n))
			break;
		
		EquateCM_CM(zCurrent, z, n);
	}
	
	if(m>1)
	{
		myPowerCM(z, zCurrent, m, n);
	}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	for(i=0;i<n;i++)
	{
		free(xt[i]);
		free(temp[i]);
		free(currentTerm[i]);
		free(zCurrent[i]);
	}
	
	free(xt);
	free(temp);
	free(currentTerm);
	free(zCurrent);
	
}/* myExpCM */

static inline void SubCM_CM(Complex **z, Complex **x, Complex **y, int n)
{
	long	i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			SubCC(z[i][j], x[i][j], y[i][j]);
	
}/* SubCM_CM */

static inline void SubCC(Complex& z, const Complex& x, const Complex& y)
{
	z.Re = x.Re - y.Re;
	z.Im = x.Im - y.Im;

}/* SubCC */



static inline void DivCM_D(Complex **z, int k, int n)
{
	long	i, j;
	
	if(k==0)
		return;
		
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			DivCC_D(z[i][j], (double)k);

}/* DivCM_D */

static inline void DivCC_D(Complex& z, double x)
{
	z.Re = z.Re/x;
	z.Im = z.Im/x;

}/* DivCC_D */

A bit of a total spaghetti as you can see.

I do not know at what stage it crashes. (it would be easier to realize that if knowing how to generate the crash without the loop in Xojo)

8

Problem actually seems to be here I think

void mySin(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp 		x, z, currentTerm, xSqr;
	long		n;
	
	z.Re = x.Re = currentTerm.Re = xRe;
	z.Im = x.Im = currentTerm.Im = xIm;
	
	Mul(xSqr, x, x); // here in case z=x
	
	zRe = z.Re;
	zIm = z.Im;
	
	for(n=2; ;n+=2)
	{
		Mul(currentTerm, xSqr, currentTerm);
		currentTerm.Re = -currentTerm.Re/(n*(n+1));
		currentTerm.Im = -currentTerm.Im/(n*(n+1));
		
		z.Re = z.Re + currentTerm.Re;
		z.Im = z.Im + currentTerm.Im;
		
		if(z.Re==zRe && z.Im==zIm)
			break;
		
		zRe = z.Re;
		zIm = z.Im;
	}
	
}/* mySin */

That for loop of his becomes endless

Where:
z.Re = -nan
zRe = -nan
z.lm=-nan
zlm=nan

I do not know what would be right there so that the math ends right, but in this case it never makes it out of the for loop.

I guess maybe nan can be interpreted as same as -nan ?

9

And better way here to reproduce is and test:


var r as Complex = new Complex(991.0359637310328935, 0.1)

r = r.Sinh

And when testing same in C++ standard lib:

#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <ranges>
 
int main()
{
    using namespace std::complex_literals;
    std::cout << std::fixed << std::setprecision(1);
 
    std::complex<double> z1(991.0359637310328935, 0.1);
   
    std::cout << "sinh" << z1 << " = " << std::sinh(z1) << std::endl;
    
    std::cout << sinh(z1).real();
}

Then that says correct result would be:
sinh(991.0,0.1) = (inf,inf)

10

Thank you Björn for giving us insights into the function of Sinh. But I can’t change it. So I somehow solved it with the workaround of writing my own code for what I intended Sinh to do. That finally works.

11

Can you share the solution for Bjön, the others.

This help too.

Marking the case Resolved also help.

12

It is no solution to the Sinh problem; it’s my own code to circumvent the utilisation of the Sinh method of the plugin. Here is the updated code:

Var bet,re,im,w,si,e As Double
Var z,x,y As Integer
Var c1,c2 As Complex
Var tabel(800,800) As Integer
Var faktor As Double = 0.00425
Var XrEcke As Double = -1.5
Var YiEcke As Double = -1.7

e=2.71828182845904523536029
c1= New Complex(1.0,0.1)

For y=0 To 800
  For x=0 To 800
    re=x*faktor+XrEcke
    im=y*faktor+YiEcke
    c2 = New Complex(re,im)
    z=0
    Do
      w=ATan(im/re)
      si=0.5*(Pow(e,w)-Pow(e,-w))
      si=Pow(si,4)
      c2=c2+si
      re=Abs(c2.Real)
      im=Abs(c2.Imag)
      c2=New Complex(re,im)
      bet=c2.Abs
      z=z+1
    Loop Until z=2
    tabel(x,y)=Round(bet)*7
  Next
Next
Break
1 Like
13

Instantiating objects is not a fast process and you’re doing it 1,283,202 times. Would you not be better using:

c2.Real = Abs( c2.Real )
c2.Imag = Abs( c2.Imag )

That way you just keep reusing the same object.

Equally, you could move the “c2 = New Complex(1,1)” outside of the loops, and then use

Unless that constructor is doing something else I’m not aware of.

14

Thank you for your useful hints @Ian_Kennedy

15

unfortunately, that doesn’t work because:

julia.hyperbolic, line 10
This method cannot accept an assigned value (it lacks an Assigns parameter).
c2.Imag=y*faktor+YiEcke
1 Like
16

Probably Immutable. Can’t be direct assigned, operations create a new result object.

17

That’s a pity, it would theoretically be much faster. Having looked at one Xojo plug-in that supports Complex it doesn’t do anything in the Constructor other than assign the real and imaginary parts.