complex math

[srvaldez]

inspired by José Roca http://www.planetsquires.com/protect/fo ... 8#msg31078
more function can of course be added.
complex.bi

Code: Select all

'
'    FILE: Complex.cpp
'  AUTHOR: Rob Tillaart
' VERSION: 0.1.09
' PURPOSE: library for Complex math for Arduino
'     URL: http:'arduino.cc/playground/Main/ComplexMath
'
' Released to the public domain
'
' translated to FreeBASIC by srvaldez 2017 sep 18

Type complex
	As Double re, im
	
	Declare Constructor()
	Declare Constructor(Byval re As Double)
	Declare Constructor(Byval re As Double, Byval im As Double)

	Declare Operator Let ( Byval rhs As Long )
	Declare Operator Let ( Byval rhs As Longint )
	Declare Operator Let ( Byval rhs As Integer )
	Declare Operator Let ( Byval rhs As Double )
	Declare Operator Let ( Byref rhs As complex )
	Declare Operator Let ( Byref rhs As String )
	Declare Operator Cast ( ) As String
End Type

Const pi =3.141592653589793
Const pi2=6.283185307179586
''' uncomment if you wish
'''Dim Shared As Const Complex i = complex(0,1)

Constructor complex()
	this.re=Cdbl(0)
	this.im=Cdbl(0)
End Constructor

Constructor complex(Byval re As Double)
	this.re=re
	this.im=Cdbl(0)
End Constructor

Constructor complex(Byval re As Double, Byval im As Double)
	this.re=re
	this.im=im
End Constructor

Operator complex.let( Byval rhs As Long )
	this.re=Cdbl(rhs)
	this.im=Cdbl(0)
End Operator

Operator complex.let( Byval rhs As Longint )
	this.re=Cdbl(rhs)
	this.im=Cdbl(0)
End Operator

Operator complex.let( Byval rhs As Integer )
	this.re=Cdbl(rhs)
	this.im=Cdbl(0)
End Operator

Operator complex.let( Byval rhs As Double )
	this.re=rhs
	this.im=Cdbl(0)
End Operator

Operator complex.let( Byref rhs As complex )
	this.re=rhs.re
	this.im=rhs.im
End Operator

Operator complex.let (Byref s As String)
	Dim As String b, a=Ucase(Ltrim(Rtrim(s)))
	Dim As Long c, i, k, l=Len(a), Sgn1=1, sgn2=1
	If Right(a,1)<>"I" Then
		Print "not a valid complex number"
		this.re=0 : this.im=0
		Exit Operator
	End If
	if a="I" or a="+I" then
		this.re=0
		this.im=1
		exit operator
	elseif a="-I" then
		this.re=0
		this.im=-1
		exit operator
	end if
	i=l-2
	While i>=0 Andalso a[i]=32
		i-=1
	Wend
	b=chr(a[i])
	If b<>"*" Then
		if b="+" or b="-" then
			if b="-" then this.im=-1 else this.im=1
			i+=1
		else
			Print "not a valid complex number"
			this.re=0 : this.im=0
			Exit Operator
		end if
	End If
	a=Rtrim(Left(a,i))
	If Left(a,1)="-" Then
		Sgn1=-1
		a=Mid(a,2)
	End If
	l=Len(a)
	c=Instr(a,"E")
	k=Instr(c+1,a,Any "+-")
	If k=c+1 Then
		k=Instr(k+1,a,Any "+-")
	End If
	If Chr(a[k-1])="-" Then
		sgn2=-1
	End If
	this.re=Sgn1*Val(Left(a,k-1))
	if this.im=0 then this.im=sgn2*Val(Mid(a,k+1))
End Operator

Operator complex.cast ( ) As String
	Dim As String c=" "
	If (this.re)<0 Then c = ""
	c += Str(this.re)
	if this.im=1 then
		c+=" + i"
	elseif this.im=-1 then
		c+=" - i"
	elseif this.im<0 Then
		c+=" - "+Str(Abs(this.im))+"*i"
	Elseif this.im>0 Then
		c+=" + "+Str(this.im)+"*i"
	End If
	Operator = c
End Operator

#Macro hypot(x,y)
	(Sqr(((x)*(x))+((y)*(y))))
#Endmacro

#Define modulus hypot

#Macro phase(x,y)
	(Atan2(y,x))
#Endmacro

Function reciprocal(Byval x As complex) As complex
    Dim As Double f = 1.0/(x.re*x.re + x.im*x.im)
    Dim As Double r = x.re*f
    Dim As Double i = -x.im*f
    Return Complex(r,i)
End Function

Function cabs(Byval x As complex) As complex
    Return Complex(Sqr(x.re*x.re+x.im*x.im),0)
End Function

Operator + (Byval lhs As complex, rhs As complex) As complex
    Return Complex(lhs.re + rhs.re, lhs.im + rhs.im)
End Operator

Operator + (Byval lhs As complex, rhs As Double) As complex
    Return Complex(lhs.re + rhs, lhs.im)
End Operator

Operator + (Byval lhs As Double, rhs As complex) As complex
    Return Complex(lhs + rhs.re, rhs.im)
End Operator

Operator - (rhs As complex) As complex 'negate
    Return Complex(-rhs.re, -rhs.im)
End Operator

Operator - (Byval lhs As complex, rhs As complex) As complex
    Return Complex(lhs.re - rhs.re, lhs.im - rhs.im)
End Operator

Operator - (Byval lhs As complex, rhs As Double) As complex
    Return Complex(lhs.re - rhs, lhs.im)
End Operator

Operator - (Byval lhs As Double, rhs As complex) As complex
    Return Complex(lhs - rhs.re, -rhs.im)
End Operator

Operator * (Byval lhs As complex, rhs As complex) As complex
    Dim As Double re = lhs.re * rhs.re - lhs.im * rhs.im
    Dim As Double im = lhs.re * rhs.im + lhs.im * rhs.re
    Return Complex(re, im)
End Operator

Operator * (Byval lhs As complex, rhs As Double) As complex
    Dim As Double re = lhs.re * rhs
    Dim As Double im = lhs.im * rhs
    Return Complex(re, im)
End Operator

Operator * (Byval lhs As Double, rhs As complex) As complex
    Dim As Double re = lhs * rhs.re
    Dim As Double im = lhs * rhs.im
    Return Complex(re, im)
End Operator

Operator / (Byval lhs As complex, rhs As complex) As complex
    Dim As Double f = 1.0/(rhs.re*rhs.re + rhs.im*rhs.im)
    Dim As Double re = (lhs.re * rhs.re + lhs.im * rhs.im) * f
    Dim As Double im = (lhs.im * rhs.re - lhs.re * rhs.im) * f
    Return Complex(re, im)
End Operator

Operator / (Byval lhs As complex, rhs As Double) As complex
    Dim As Double f = 1.0/(rhs*rhs)
    Dim As Double re = (lhs.re * rhs) * f
    Dim As Double im = (lhs.im * rhs) * f
    Return Complex(re, im)
End Operator

Operator / (Byval lhs As double, rhs As complex) As complex
    Dim As Double f = 1.0/(rhs.re*rhs.re + rhs.im*rhs.im)
    Dim As Double re = (lhs * rhs.re) * f
    Dim As Double im = (-lhs * rhs.im) * f
    Return Complex(re, im)
End Operator

Function csquare(Byval x As complex) As complex
    Dim As Double re = x.re * x.re - x.im * x.im
    Dim As Double im = 2 * x.re * x.im
    Return Complex(re, im)
End Function

Function csqr(Byval x As complex) As complex
    Dim As Double m = hypot(x.re,x.im)
    Dim As Double re = Sqr(0.5 * (m+x.re))
    Dim As Double im = Sqr(0.5 * (m-x.re))
    If x.im < 0 Then im = -im
    Return Complex(re, im)
End Function

Function cexp(Byval x As complex) As complex
    Dim As Double e = Exp(x.re)
    Return Complex(e * Cos(x.im), e * Sin(x.im))
End Function

Function clog(Byval x As complex) As complex
    Dim As Double m = modulus(x.re, x.im)
    Dim As Double p = phase(x.re, x.im)
    If p > pi Then p -= pi2
    Return Complex(Log(m), p)
End Function

Function cpow(Byval x As complex, Byval c As Const Complex) As complex
    Dim As Complex t = clog(x)
    t = t * c
    Return cexp(t)
End Function

Function clogn(Byval x As complex, Byval c As Const Complex) As complex
    Return clog(x)/clog(c)
End Function

Function clog10(Byval x As complex) As complex
    Return clogn(x, 10)
End Function

Function sinh(Byval x As Double) As Double
	Dim As Double y=Exp(x)
	Return (y-1.0/y)*0.5
End Function

Function cosh(Byval x As Double) As Double
	Dim As Double y=Exp(x)
	Return (y+1.0/y)*0.5
End Function

Function tanh(Byval x As Double) As Double
	Dim As Double y=Exp(x)
	Return (y-1.0/y)/(y+1.0/y)
End Function

Function csin(Byval x As complex) As complex
    Dim As Double s = Sin(x.re)
    Dim As Double c = Sqr(1.0-s*s)
    Return Complex(s * cosh(x.im), c * sinh(x.im))
End Function

Function ccos(Byval x As complex) As complex
    Dim As Double s = Sin(x.re)
    Dim As Double c = Sqr(1.0-s*s)
    Return Complex(c * cosh(x.im), -s * sinh(x.im))
End Function

Function ctan(Byval x As complex) As complex
    '/* faster but 350 bytes longer!!
    Dim As Double s = Sin(x.re)
    Dim As Double c = Cos(x.re)
    Dim As Double sh = sinh(x.im)
    Dim As Double ch = cosh(x.im)
    '' return Complex(s*ch, c*sh) / Complex(c*ch, -s*sh)
    Dim As Double r0 = s*ch
    Dim As Double i0 = c*sh
    Dim As Double cre = c*ch
    Dim As Double cim = -s*sh
    Dim As Double f = 1.0/(cre*cre + cim*cim)
    Dim As Double r = r0 * cre + i0 * cim
    Dim As Double i = r0 * cim - i0 * cre
    Return Complex(r * f, -i * f)
    
    '''return c_sin() / c_cos();
End Function

Dim Shared As Const Complex one = 1

Function gonioHelper1(Byval x As complex, Byval mode As Const Byte) As complex
    Dim As Complex c = csqr(one - csquare(x))
    If mode = 0 Then
        c = c + x * Complex(0,-1)
    Else
        c = x + c * Complex(0,-1)
    End If
    c = clog(c) * Complex(0,1)
    Return c
End Function

Function casin(Byval x As complex) As complex
    Return gonioHelper1(x, 0)
End Function

Function cacos(Byval x As complex) As complex
    Return gonioHelper1(x, 1)
End Function

Function catan(Byval x As complex) As complex
    Return (Complex(0,-1) * clog(Complex(x.re, x.im - 1)/Complex(-x.re, -x.im - 1))) * 0.5
End Function

Function ccsc(Byval x As complex) As complex
    Return one / csin(x)
End Function

Function csec(Byval x As complex) As complex
    Return one / ccos(x)
End Function

Function ccot(Byval x As complex) As complex
    Return one / ctan(x)
End Function

Function cacsc(Byval x As complex) As complex
    Return one / casin(x)
End Function

Function casec(Byval x As complex) As complex
    Return one / cacos(x)
End Function

Function cacot(Byval x As complex) As complex
    Return one / catan(x)
End Function

'HYPERBOLICUS I
function csinh(Byval x As complex) as complex
    dim as complex y=cexp(x)
    return (y-reciprocal(y))*0.5
end function

function ccosh(Byval x As complex) as complex
    dim as complex y=cexp(x)
    return (y+reciprocal(y))*0.5
end function

function ctanh(Byval x As complex) as complex
	return csinh(x) / ccosh(x)
end function

function gonioHelper2(Byval x As complex, byval mode as const byte) as complex
    dim as Complex c = csquare(x)
    if mode = 0 then
        c = c+1
    else
        c = c-1
    end if
    c = clog(x + csqr(c))
    return c
end function

function casinh(Byval x As complex) as complex
	return gonioHelper2(x,0)
end function

function cacosh(Byval x As complex) as complex
	return clog(x+csqr(x-one)*csqr(x+one))
	'return gonioHelper2(x,1)
end function

function catanh(Byval x As complex) as complex
    dim as Complex c = clog(x + one)
    c = c - clog(-(x - one))
    return c * 0.5
end function

'HYPERBOLICUS II
function ccsch(Byval x As complex) as complex
	return one / csinh(x)
end function

function csech(Byval x As complex) as complex
	return one / ccosh(x)
end function

function ccoth(Byval x As complex) as complex
	return one / ctanh(x)
end function

function cacsch(Byval x As complex) as complex
	return casinh(one/x)
end function

function casech(Byval x As complex) as complex
	return cacosh(one/x)
end function

function cacoth(Byval x As complex) as complex
	return catanh(one/x)
end function

[edit] changed the string casting functions to allow i without a factor, so "1+i" is accepted.
test-complex.bas

Code: Select all

#include "complex.bi"
Dim As complex x=complex(2,5),y=complex(3,7),z
Dim As Double r=1

z=x-y
Print "z = ";x;" - ";y;" = ";z
z=ctan(z)
z=catan(z)
Print "catan(ctan(";z;")) = ";z
Print "assignment from string"
z="-1e+10 +   1e-100   *    i   "
Print "z = ";Chr(34)+"-1e+10 +   1e-100   *    i   "+Chr(34);" = ";z
z="+100 - 0.0000001 * i"
Print "z = ";Chr(34)+"+100 - 0.0000001 * i"+Chr(34);" = ";z
Print "passing a double value to a complex function"
Print "dim r as double = 1 : print catan(r) -> ";catan(r)

[srvaldez]

I have a question for the experts, the complex inverse hyperbolic cosine in the original code was

Code: Select all

function cacosh(Byval x As complex) as complex
	return gonioHelper2(x,1)
end function

if z=-1-2*i
then cacosh(ccosh(z)) ->-1-2*i
however, Mathematica, Maple, Maxima, pari/gp all give 1+2*i, so what's the right answer?
I changed the code so that it conforms with the above mentioned software, but it makes no sense to me.

[jdebord]

With FBMath I get (1 + 2i), too.

[srvaldez]

hello jdebord
could this be a case where both answers a correct?
[edit] the answer is yes, it's a multi-valued function http://scipp.ucsc.edu/~haber/webpage/arc.pdf
but the choice is made in favor of the positive answer.

[dodicat]

From first principles.

Code: Select all

Type complex
    As Double re,im
    declare operator cast() as string
End Type

operator complex.cast() as string
return str(re) +" , " +str(im)+" i"
end operator


Operator *(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)'n
End Operator

Operator +(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re+n2.re,n1.im+n2.im)'n
End Operator

Operator -(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re-n2.re,n1.im-n2.im)
End Operator

Operator /(n1 As complex,n2 As complex) As complex
Dim As Double d = n2.re*n2.re+n2.im*n2.im
Return Type<complex>((n1.re*n2.re+n1.im*n2.im)/d,(n1.im*n2.re - n1.re*n2.im)/d)
End Operator

Function CExp( Z As Complex) As Complex
    Return Type<Complex> (Exp(z.re) * Cos(Z.im), Exp(z.re) * Sin(Z.im))
End Function 

Function Clog(z As complex) As complex
    Dim As Double r=Sqr(z.re*z.re+z.im*z.im)
    Dim As Double theta=Atan2(z.im,z.re)
    Return  Type<complex>(Log(r),theta)
End Function

Operator ^(a As complex,b As complex) As complex
If a.re = 0 And a.im = 0 Then
    If b.re = 0 And b.im = 0 Then
        Return Type<Complex> (1, 0) 
    Else
        Return Type<Complex> (0, 0)
    End If  
End If
Return CExp(B*CLog(A))
End Operator

Function cosh(Byval x As Double) As Double
    cosh = (Exp(x) + Exp(-x)) / 2
End Function

Function sinh(Byval x As Double) As Double
    sinh = (Exp(x) - Exp(-x)) / 2 
End Function
 
Function Ccos(z As complex) As complex
    Return Type<complex>(Cos(z.re)*cosh(z.im),-Sin(z.re)*sinh(z.im))
End Function

function Ccosh(z as complex) as complex
    dim as complex i=type(0,1)
    return ccos(i*z)
end function


FUNCTION Carccos(byval z AS complex) AS complex
    dim as complex tmp=clog(z+ (z*z-type(1,0))^ type(.5,0))
    return type(0,-1)*tmp
END FUNCTION


function Carccosh(z as complex) as complex
    return (((z-type(1,0))^type(.5,0)) / ((type(1,0)-z)^type(.5,0)))*Carccos(z) 
    end function

dim as complex z=type(1,-2)

dim as complex A=ccosh(z)

print z
print A
print cArccosh(a)

sleep



 

Which is a bit silly because jdeborg has a library for this already.

[dodicat]

You can try values here :
http://en.cppreference.com/w/cpp/numeric/complex/acosh
Gcc (G++) seems a crazy environment.
It seems each build has no backward compatibility.
to run this uses c++17 and some concepts??
And jbeborg is jdebord (sorry, typo)

[srvaldez]

thank you dodicat for your example and link, my reason for posting this is because I am a fan of public domain code, though I usually like to give credit where credit is due.

[jdebord]

@srvaldez : Thank you for the Haber reference. Very useful !

@dodicat : I checked the C++ examples with FBMath. For the first two I get the same result : (0.000000,1.047198). Probably FB does not make a difference between 0.000000 and -0.000000

For the two other examples I get the same results than C++

By the way, while making an example program for FBMath I noticed that, in order to retrieve the original complex number, you have to apply the "arc" function first, e. g. ccosh(cacosh(z)) = z

[dodicat]

Perhaps, if it known that a complex function has two answers, the type complex could be adjusted to serve.
Using (from the C++ site)
log(z + sqr(z+1)*sqr(z-1) for the complex arccosh function, and a second answer being the same formula with z.im part being flipped in sign.

The function could return this custom tweak.
edit : any number of answers (e.g. complex roots)

Code: Select all

Type complex
    As Double re,im
    As boolean MoreSolutions           'for more solutions this=true
    As Double re2(Any),im2(Any)        'the other solutions
    Declare Constructor(As Double=0,As Double=0)
    Declare Operator Cast() As String
End Type

Constructor complex(real As Double,imag As Double)
re=real
im=imag
End Constructor

#include "crt/stdio.bi"

Function CRound(Byval x As Double,Byval precision As Integer=30) As String
    If precision>30 Then precision=30
    Dim As zstring * 40 z:Var s="%." &str(Abs(precision)) &"f"
    sprintf(z,s,x)
    Return Rtrim(Rtrim(z,"0"),".")
End Function 

Operator complex.cast() As String
Dim As long places=15,tb=30
Dim As String sr,si
If re>=0 Then sr="+" Else sr=""
If im>=0 Then si="+" Else si=""
If MoreSolutions=false Then
    Print sr+Cround(re,places);Tab(tb);si+Cround(im,places); " i"
Else
    Print sr+str(re);Tab(tb);si+Cround(im,places); " i"
    For n As Long=1 To Ubound(re2)
        If re2(n)>=0 Then sr="+" Else sr=""
        If im2(n)>=0 Then si="+" Else si=""
        Print sr+Cround(re2(n),places);Tab(tb);si+Cround(im2(n),places); " i"
    Next n
    Return ""
End If
End Operator


Operator *(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)'n
End Operator

Operator +(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re+n2.re,n1.im+n2.im)'n
End Operator

Operator -(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re-n2.re,n1.im-n2.im)
End Operator

Operator /(n1 As complex,n2 As complex) As complex
Dim As Double d = n2.re*n2.re+n2.im*n2.im
Return Type<complex>((n1.re*n2.re+n1.im*n2.im)/d,(n1.im*n2.re - n1.re*n2.im)/d)
End Operator

Function CExp( Z As Complex) As Complex
    Return Type<Complex> (Exp(z.re) * Cos(Z.im), Exp(z.re) * Sin(Z.im))
End Function 

Function Clog(z As complex) As complex
    Dim As Double r=Sqr(z.re*z.re+z.im*z.im)
    Dim As Double theta=Atan2(z.im,z.re)
    Return  Type<complex>(Log(r),theta)
End Function

Operator ^(a As complex,b As complex) As complex
If a.re = 0 And a.im = 0 Then
    If b.re = 0 And b.im = 0 Then
        Return Type<Complex> (1, 0) 
    Else
        Return Type<Complex> (0, 0)
    End If  
End If
Return CExp(B*CLog(A))
End Operator

Function cosh(Byval x As Double) As Double
cosh = (Exp(x) + Exp(-x)) / 2
End Function

Function sinh(Byval x As Double) As Double
 sinh = (Exp(x) - Exp(-x)) / 2 
End Function

Function Ccos(z As complex) As complex
    Return Type<complex>(Cos(z.re)*cosh(z.im),-Sin(z.re)*sinh(z.im))
End Function

Function Ccosh(z As complex) As complex
    Dim As complex i=Type(0,1)
    Return ccos(i*z)
End Function

Function Carccosh(Byval z As complex) As complex
    Dim As complex ans,tmp
    Ans= clog(z+ (z+Type(1,0))^Type(.5,0) *   (z-Type(1,0))^Type(.5,0))'main return
    ans.MoreSolutions=true   'always two solutions for this function
    If ans.MoreSolutions Then
        Redim ans.re2(1 To 1):Redim ans.im2(1 To 1)
        z.im=-z.im 'do what is needed for solution 2
        tmp= clog(z+ (z+Type(1,0))^Type(.5,0) *   (z-Type(1,0))^Type(.5,0))
        ans.re2(1)=tmp.re:ans.im2(1)=tmp.im
    End If
    Return ans
End Function

Function Roots(Z As Complex,N As Long) As Complex
    Dim As complex p
    If N <= 0 Then
        Return Type<Complex> (0, 0)
    End If
    If (Sqr(Z.re*Z.re+Z.im*Z.im)^(1/N))=0 Or (Sqr(Z.re*Z.re+Z.im*Z.im)^(1/N))<0 Then 
        Return Type<complex>(0,0)
    End If
    if n=1 then return z
    Dim As Double pi2=8*Atn(1)
    'the primary root
    p= Type((Sqr(Z.re*Z.re+Z.im*Z.im)^(1/n))*Cos(((Atan2(Z.im,Z.re))/n)),(Sqr(Z.re*Z.re+Z.im*Z.im)^(1/n))*Sin(((Atan2(Z.im,Z.re)/n))))

    Redim p.re2(1 To n-1):Redim p.im2(1 To n-1)
    p.moresolutions=true 'the other roots
    For k As Long=1 To n-1
        p.re2(k)=(Sqr(Z.re*Z.re+Z.im*Z.im)^(1/N))*Cos(((Atan2(Z.im,Z.re)+K*pi2)/N))
        p.im2(k)=(Sqr(Z.re*Z.re+Z.im*Z.im)^(1/N))*Sin(((Atan2(Z.im,Z.re)+K*pi2)/N))
    Next k
    Return p
End Function
'==============

Dim As complex z=Type(-1,-2)


Print "number "
Print z
'print Ccosh(z)

Print
Print "answers"
var a= cArccosh(ccosh(z))
print a 
Print "Primary answer parts"
print a.re,a.im
print "invert test"
Print ccosh(carccosh(z))
Print
Print
print "cube roots of "
print type<complex>(-1,-2)
var r= roots(Type(-1,-2),3)
print "are"
print r
print
print "Primary cube root parts"
print r.re,r.im

Sleep


  

[srvaldez]

dodicat, shouldn't the two answers be 1 + 2i and -1 - 2i ?
excuse me if I am wrong, don't know much about complex numbers. :-)

[dodicat]

from the c++ code
http://en.cppreference.com/w/cpp/numeric/complex/acosh
Putting in the complex cosh of (-1,-2)
that is:
-0.64214812471552 , 1.068607421382778 i
and then
-0.64214812471552 , -1.068607421382778 i

gives
1+2i
1-2i

That's what I am going by.
changing the sign of the .im part gives the other answer.(seems!)
But I am no expert in the inverse complex hyperbolic cosine, never used it until now.
And with any luck, it'll stay this way.