Translate C code to Xojo

Hello,
Any body who master a C language to translate math function (non linear equation root) to xojo ?

How many lines of code?
Can you post the code here?

Just post the code so we can have a look.

#include <math.h>
#define MAXIT 100

float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2,
float xacc)
{
void nrerror(char error_text[]);
int j;
float df,dx,dxold,f,fh,fl;
float temp,xh,xl,rts;

(*funcd)(x1,&fl,&df);
(*funcd)(x2,&fh,&df);
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
	nrerror("Root must be bracketed in rtsafe");
if (fl == 0.0) return x1;
if (fh == 0.0) return x2;
if (fl < 0.0) {
	xl=x1;
	xh=x2;
} else {
	xh=x1;
	xl=x2;
}
rts=0.5*(x1+x2);
dxold=fabs(x2-x1);
dx=dxold;
(*funcd)(rts,&f,&df);
for (j=1;j<=MAXIT;j++) {
	if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
		|| (fabs(2.0*f) > fabs(dxold*df))) {
		dxold=dx;
		dx=0.5*(xh-xl);
		rts=xl+dx;
		if (xl == rts) return rts;
	} else {
		dxold=dx;
		dx=f/df;
		temp=rts;
		rts -= dx;
		if (temp == rts) return rts;
	}
	if (fabs(dx) < xacc) return rts;
	(*funcd)(rts,&f,&df);
	if (f < 0.0)
		xl=rts;
	else
		xh=rts;
}
nrerror("Maximum number of iterations exceeded in rtsafe");
return 0.0;

}
#undef MAXIT

sorry , i use a blockquote but something wrong , in the begenning a difference ocurs (code format).
The script is Ok .

remove “;”
remove brackets { and }
put “end if” where they belong
end for where it should be
replace == with =
and it should be good to go almost.

It looks like it’s using the secant method:

I’ve already written Xojo code to do this. It’ll save having to wade through the horror of C code. Too late tonight, but I’ll dig it out tomorrow and post it.

#include <math.h>
#define MAXIT 100

float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2,
	float xacc)
{
	void nrerror(char error_text[]);
	int j;
	float df,dx,dxold,f,fh,fl;
	float temp,xh,xl,rts;

	(*funcd)(x1,&fl,&df);
	(*funcd)(x2,&fh,&df);
	if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
		nrerror("Root must be bracketed in rtsafe");
	if (fl == 0.0) return x1;
	if (fh == 0.0) return x2;
	if (fl < 0.0) {
		xl=x1;
		xh=x2;
	} else {
		xh=x1;
		xl=x2;
	}
	rts=0.5*(x1+x2);
	dxold=fabs(x2-x1);
	dx=dxold;
	(*funcd)(rts,&f,&df);
	for (j=1;j<=MAXIT;j++) {
		if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
			|| (fabs(2.0*f) > fabs(dxold*df))) {
			dxold=dx;
			dx=0.5*(xh-xl);
			rts=xl+dx;
			if (xl == rts) return rts;
		} else {
			dxold=dx;
			dx=f/df;
			temp=rts;
			rts -= dx;
			if (temp == rts) return rts;
		}
		if (fabs(dx) < xacc) return rts;
		(*funcd)(rts,&f,&df);
		if (f < 0.0)
			xl=rts;
		else
			xh=rts;
	}
	nrerror("Maximum number of iterations exceeded in rtsafe");
	return 0.0;
}
#undef MAXIT

I said earlier that it appears to use the secant method, but I now think it uses the bisection method, which is similar but less flexible than the secant method. In the bisection method, the two input values must bracket the root; otherwise it won’t converge. Here is my code for the secant method.

Function SecantRoot(xA As double, xB As double, errTolerance As double) as double
  'Secant method for finding root of function
  'Input parameters are two x values which
  'preferably bracket the root value, and if not, are at least near the root
  '
  'SecantTestFn is the name of the function whose root is to be found
  '
  const maxit = 100
  
  dim x0 As Double = xA
  dim x1 As Double = xB
  if SecantTestFn(x0)=0 then return x0
  if SecantTestFn(x1)=0 then return x1
  
  dim xn As Double
  for i as Integer = 1 to maxit
    xn = x1 - SecantTestFn(x1)*(x1-x0)/(SecantTestFn(x1)-SecantTestFn(x0))
    if SecantTestFn(xn)=0 or abs(1-x1/xn)<errTolerance then return xn
    x0=x1
    x1=xn
  next
  'Failed to find root within specified tolerance. So, return closest value
  return xn
End Function

Much of the confusing clutter in the posted C code is because one of the calling parameters is a pointer to the function to be minimized. My code has the function name hard coded. This could be changed by passing the function as a delegate.

1 Like

RTSafe is bisection + Newton-Raphson

See pdf page 390 (printed book 366): https://www.cec.uchile.cl/cinetica/pcordero/MC_libros/NumericalRecipesinC.pdf

Yes,

Yes the function posted rtsafe is from NumericalRecipesC.
the algo use the rtsafe to find the root, the author recommand rtsafe instead of NewtonRaphson
.

In the algo the author use a Lagrange Muliplier (optimization) wich is the smalest root of an a 6th deg equation.
I think rtsafe is used to insure the convergeance if a local extremum appears.

gooood , thank’s.

Ah, Numerical Recipes, of course! No wonder the code is unreadable.

rtsafe converted to Xojo. Not tested, and probably won’t work without some modifications, but this should give a rough idea. I have almost no experience with delegates. The last time I used them was many years ago. So hopefully someone will jump in here and correct what I’ve done.

function rtsafe (funcd as delegate, x1 as double, x2 as double) as double
  dim j as integer
  dim df,dx,dxold,f,fh,fl as double
  dim temp,xh,xl, rts as double
  funcd.invoke(x1,fl,df)
  funcd.invoke(x2,fh,df)
  
  if ((fl > 0.0 and fh > 0.0) or (fl < 0.0 and fh < 0.0)) then
    nrerror("Root must be bracketed in rtsafe")
    end if
  if (fl = 0.0) then return x1
  if (fh = 0.0) then return x2
  if (fl < 0.0) then
    xl = x1
    xh = x2
  else
    xh = x1
    xl = x2
  end if

  rts = 0.5 *(x1 + x2)
  dxold = abs(x2 - x1)
  dx = dxold
  funcd.invoke(rts,f,df)

  for j = 1 to MAXIT
    if ((((rts - xh) *df - f) *((rts - xl) *df - f) > 0.0) or (fabs(2.0 *f) > fabs(dxold *df))) then
      dxold = dx
       dx = 0.5 *(xh - xl)
       rts = xl + dx
       if (xl = rts) then return rts
     else
      dxold = dx
       dx = f / df
       temp = rts
       rts = rts - dx
       if (temp = rts) then return rts
     end if
    if (fabs(dx) < xacc) then return rts
     funcd.invoke(rts, f, df)
    if (f < 0.0) then
      xl = rts
     else
      xh = rts
     end if
  next
end function

Here is a more complete set of the Robert’s conversion.
I felt I needed to make clear the necessity of the ByRefs or it would fail miserably case the user misses it.


Public Sub funcdDelegate(value As Double, Byref result As double, Byref derivative As Double)


Public Function rtsafe(funcDeriv As funcdDelegate, x1 As Double, x2 As Double, maxIter As Integer = 100, accuracy As Double = 3E-8) as Double
  Var j As Integer
  Var df,dx,dxold,f,fh,fl As Double
  Var temp,xh,xl, rts As Double
  
  funcDeriv.Invoke(x1,fl,df)
  funcDeriv.Invoke(x2,fh,df)
  
  If ((fl > 0.0 And fh > 0.0) Or (fl < 0.0 And fh < 0.0)) Then
    Var ex As New RuntimeException
    ex.Message = "Root must be bracketed in rtsafe"
    Raise ex
    Return 0.0
  End If
  
  If (fl = 0.0) Then Return x1
  If (fh = 0.0) Then Return x2
  
  If (fl < 0.0) Then
    xl = x1
    xh = x2
  Else
    xh = x1
    xl = x2
  End If
  
  rts = 0.5 *(x1 + x2)
  dxold = abs(x2 - x1)
  dx = dxold
  
  funcDeriv.Invoke(rts,f,df)
  
  For j = 1 To maxIter
    
    If (((rts - xh) *df - f) *((rts - xl) *df - f) > 0.0) Or (abs(2.0 *f) > abs(dxold *df)) Then
      dxold = dx
      dx = 0.5 *(xh - xl)
      rts = xl + dx
      If (xl = rts) Then Return rts
    Else
      dxold = dx
      dx = f / df
      temp = rts
      rts = rts - dx
      If (temp = rts) Then Return rts
    End If
    
    If (abs(dx) < accuracy) Then Return rts
    
    funcDeriv.Invoke(rts, f, df)
    
    If (f < 0.0) Then
      xl = rts
    Else
      xh = rts
    End If
    
  Next
  
  Var ex As New RuntimeException
  ex.Message = "Maximum number of iterations exceeded in rtsafe"
  Raise ex
  
  Return 0.0
  
End Function


Public Sub funcd(value As Double, Byref result As double, Byref derivative As Double)
  // funcd example matching the type funcdDelegate
  // Do your math here receiving "value" and returning the result and the first derivative
End Sub


Public Sub main()
  
  // Example passing the delegated function funcd above
  Var root as double = rtsafe((New funcdDelegate(AddressOf funcd)), x1, x2)
  
End Sub

2 Likes

Rober Weaver,
the most work is done.
many thank’s.

Thanks. I was going to mention the ByRef requirement, but wasn’t sure of the syntax for the delegate definition. So, I thought it best to leave well enough alone, and let someone else do it properly.
Also, I see you caught the fabs() and changed it to abs(). I just noticed that I’d missed it.

1 Like

Hello Rick,
that what i was trying to do …
thank for you help.

1 Like

Hi Rick,
example to use rtsafe().

// Function F(x)
Fx(X as double) as double
return cos(x)-x * x * x //NewtonRaphson failed to converge to the root with
guess = 0

//First derivative
dFdx(X as double) as double
return -sin(x)-3 * x * x

With bisection [0,1] the convergeance done, root : 0.865479.
Could you show me by an example how to apply this example in rtsafe() ?
Not understand funcdDelegate() …