--[[---------------------------------------------------------------------------- LIBRARY FOR COMPLEX NUMBERS Philippe CASTAGLIOLA Université de Nantes & IRCCyN UMR CNRS 6597 WWW: http://philippe.castagliola.free.fr/ --]]---------------------------------------------------------------------------- local error_=error local function error(fmt,...) error_(string.format(fmt,...)) end local assert_=assert local function assert(b,fmt,...) return fmt and assert_(b,string.format(fmt,...)) or assert_(b) end type_=type -------------------------------------------------------------------------------- local math_={} for _,v in pairs{"abs","acos","asin","atan","cos","exp","log","pow","sin","sqrt","tan"} do math_[v]=math[v] end -------------------------------------------------------------------------------- function math.acosh(z) -------------------------------------------------------------------------------- return math_.log(z+math_.sqrt(z*z-1)) end -------------------------------------------------------------------------------- function math.asinh(z) -------------------------------------------------------------------------------- return math_.log(z+math_.sqrt(z*z+1)) end -------------------------------------------------------------------------------- function math.atanh(z) -------------------------------------------------------------------------------- return math_.log((1+z)/(1-z))/2 end for _,v in pairs{"acosh","asinh","atanh"} do math_[v]=math[v] end -- SHOULD BE REMOVED IN LUA 5.2 -- START -------------------------------------------------------------------------------- function math.cosh(z) -------------------------------------------------------------------------------- return (math_.exp(z)+math_.exp(-z))/2 end -------------------------------------------------------------------------------- function math.sinh(z) -------------------------------------------------------------------------------- return (math_.exp(z)-math_.exp(-z))/2 end -------------------------------------------------------------------------------- function math.tanh(z) -------------------------------------------------------------------------------- local e=math_.exp(2*z) return (e-1)/(e+1) end for _,v in pairs{"cosh","sinh","tanh"} do math_[v]=math[v] end table.unpack=unpack -- SHOULD BE REMOVED IN LUA 5.2 -- END -------------------------------------------------------------------------------- complex={} complex.__index=complex complex.version="2011-01-05" math.i=setmetatable({r=0,i=1},complex) -------------------------------------------------------------------------------- local is={} function is.complex(c) return getmetatable(c)==complex end -------------------------------------------------------------------------------- function type(x) return getmetatable(x)==complex and "complex" or type_(x) end -------------------------------------------------------------------------------- local check={} function check.complex(c,i) local t=type(c) assert(getmetatable(c)==complex,"bad argument #%d (complex number expected, got %s)",i,t) end function check.number(n,i) assert(tonumber(n),"bad argument #%d (number expected, got %s)",i,type(n)) end -------------------------------------------------------------------------------- local string_={} string_.format=string.format function string.format(fmt,...) local arg_={} for i=1,#arg do if is.complex(arg[i]) then table.insert(arg_,arg[i].r) table.insert(arg_,arg[i].i) else table.insert(arg_,arg[i]) end end return string_.format(fmt,table.unpack(arg_)) end -------------------------------------------------------------------------------- local function hypot(xr,xi) -------------------------------------------------------------------------------- local xr=math_.abs(xr) local xi=math_.abs(xi) if xr==0 then return xi end if xi==0 then return xr end if xr>xi then return xr*math_.sqrt(1+(xi/xr)^2) else return xi*math_.sqrt(1+(xr/xi)^2) end end -------------------------------------------------------------------------------- local function loghypot(xr,xi) -------------------------------------------------------------------------------- local xr=math_.abs(xr) local xi=math_.abs(xi) if xr==0 then return math_.log(xi) end if xi==0 then return math_.log(xr) end if xr>xi then return math_.log(xr)+math_.log(1+(xi/xr)^2)/2 else return math_.log(xi)+math_.log(1+(xr/xi)^2)/2 end end -------------------------------------------------------------------------------- local function sqrthypot(xr,xi) -------------------------------------------------------------------------------- local xr=math_.abs(xr) local xi=math_.abs(xi) if xr==0 then return math_.sqrt(xi) end if xi==0 then return math_.sqrt(xr) end if xr>xi then return math_.sqrt(xr*math_.sqrt(1+(xi/xr)^2)) else return math_.sqrt(xi*math_.sqrt(1+(xr/xi)^2)) end end -------------------------------------------------------------------------------- function complex.abs(x) -------------------------------------------------------------------------------- check.complex(x,1) return hypot(x.r,x.i) end -------------------------------------------------------------------------------- function math.abs(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then return math_.abs(z) elseif is.complex(x) then return complex.abs(x) else error("bad argument #1 to 'abs' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.__add(x1,x2) -------------------------------------------------------------------------------- local x if type_(x1)=="number" then x={r=x1+x2.r,i=x2.i} elseif type_(x2)=="number" then x={r=x1.r+x2,i=x1.i} else check.complex(x1,1) check.complex(x2,2) x={r=x1.r+x2.r,i=x1.i+x2.i} end return setmetatable(x,complex) end -------------------------------------------------------------------------------- function complex.acos(x) -------------------------------------------------------------------------------- check.complex(x,1) local x1r=1-x.r*x.r+x.i*x.i local x1i=-2*x.r*x.i local a=math.atan2(x1i,x1r)/2 local m=sqrthypot(x1r,x1i) local x2r=m*math_.cos(a)-x.i local x2i=m*math_.sin(a)+x.r local y={r=math.pi/2-math.atan2(x2i,x2r),i=loghypot(x2r,x2i)} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.acos(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then if -1<=z and z<=1 then return math_.acos(z) else x={r=z,i=0} setmetatable(x,complex) end end if is.complex(x) then return complex.acos(x) else error("bad argument #1 to 'acos' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.acosh(x) -------------------------------------------------------------------------------- check.complex(x,1) local x1r=x.r*x.r-x.i*x.i-1 local x1i=2*x.r*x.i local a=math.atan2(x1i,x1r)/2 local m=sqrthypot(x1r,x1i) local x2r=m*math_.cos(a)+x.r local x2i=m*math_.sin(a)+x.i local y={r=loghypot(x2r,x2i),i=math.atan2(x2i,x2r)} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.acosh(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then if 1<=z then return math_.acosh(z) else x={r=z,i=0} setmetatable(x,complex) end end if is.complex(x) then return complex.acosh(x) else error("bad argument #1 to 'acosh' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.arg(x) -------------------------------------------------------------------------------- check.complex(x,1) return math.atan2(x.i,x.r) end -------------------------------------------------------------------------------- function complex.asin(x) -------------------------------------------------------------------------------- check.complex(x,1) local x1r=1-x.r*x.r+x.i*x.i local x1i=-2*x.r*x.i local a=math.atan2(x1i,x1r)/2 local m=sqrthypot(x1r,x1i) local x2r=m*math_.cos(a)-x.i local x2i=m*math_.sin(a)+x.r local y={r=math.atan2(x2i,x2r),i=-loghypot(x2r,x2i)} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.asin(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then if -1<=z and z<=1 then return math_.asin(z) else x={r=z,i=0} setmetatable(x,complex) end end if is.complex(x) then return complex.asin(x) else error("bad argument #1 to 'asin' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.asinh(x) -------------------------------------------------------------------------------- check.complex(x,1) local x1r=x.r*x.r-x.i*x.i+1 local x1i=2*x.r*x.i local a=math.atan2(x1i,x1r)/2 local m=sqrthypot(x1r,x1i) local x2r=m*math_.cos(a)+x.r local x2i=m*math_.sin(a)+x.i local y={r=loghypot(x2r,x2i),i=math.atan2(x2i,x2r)} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.asinh(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then return math_.asinh(z) elseif is.complex(x) then return complex.asinh(x) else error("bad argument #1 to 'asinh' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.atan(x) -------------------------------------------------------------------------------- check.complex(x,1) local x1r=1+x.i local x1i=-x.r local x2r=1-x.i local x2i=x.r local y={r=(math.atan2(x2i,x2r)-math.atan2(x1i,x1r))/2, i=(loghypot(x1r,x1i)-loghypot(x2r,x2i))/2} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.atan(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then return math_.atan(z) elseif is.complex(x) then return complex.atan(x) else error("bad argument #1 to 'atan' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.atanh(x) -------------------------------------------------------------------------------- check.complex(x,1) local x1r=1+x.r local x1i=x.i local x2r=1-x.r local x2i=-x.i local m2=x2r*x2r+x2i*x2i local x3r=(x1r*x2r+x1i*x2i)/m2 local x3i=(x1i*x2r-x1r*x2i)/m2 local y={r=loghypot(x3r,x3i)/2,i=math.atan2(x3i,x3r)/2} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.atanh(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then if -1=1 and math.floor(n)==n,"bad argument #2 (must be an integer >= 1, got %g)",n) local r={} local a=math.atan2(x.i,x.r) local m=(hypot(x.r,x.i))^(1/n) for k=0,n-1 do table.insert(r,{r=m*math_.cos((a+2*k*math.pi)/n), i=m*math_.sin((a+2*k*math.pi)/n)}) setmetatable(r[k+1],complex) end return r end -------------------------------------------------------------------------------- function complex.sin(x) -------------------------------------------------------------------------------- check.complex(x,1) local e1=math_.exp(-x.i) local e2=math_.exp(x.i) local y={r=(e1+e2)*math_.sin(x.r)/2, i=(e2-e1)*math_.cos(x.r)/2} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.sin(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then return math_.sin(z) elseif is.complex(x) then return complex.sin(x) else error("bad argument #1 to 'sin' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.sinh(x) -------------------------------------------------------------------------------- check.complex(x,1) local e1=math_.exp(x.r) local e2=math_.exp(-x.r) local y={r=(e1-e2)*math_.cos(x.i)/2, i=(e1+e2)*math_.sin(x.i)/2} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.sinh(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then return math_.sinh(z) elseif is.complex(x) then return complex.sinh(x) else error("bad argument #1 to 'sinh' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.sqrt(x) -------------------------------------------------------------------------------- check.complex(x,1) local a=math.atan2(x.i,x.r)/2 local m=sqrthypot(x.r,x.i) local y={r=m*math_.cos(a),i=m*math_.sin(a)} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.sqrt(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then if z>=0 then return math_.sqrt(z) else x={r=z,i=0} setmetatable(x,complex) end end if is.complex(x) then return complex.sqrt(x) else error("bad argument #1 to 'sqrt' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.__sub(x1,x2) -------------------------------------------------------------------------------- local x if type_(x1)=="number" then x={r=x1-x2.r,i=-x2.i} elseif type_(x2)=="number" then x={r=x1.r-x2,i=x1.i} else check.complex(x1,1) check.complex(x2,2) x={r=x1.r-x2.r,i=x1.i-x2.i} end return setmetatable(x,complex) end -------------------------------------------------------------------------------- function complex.tan(x) -------------------------------------------------------------------------------- check.complex(x,1) local e=math_.exp(-2*x.i) local x1r=e*math_.cos(2*x.r)-1 local x2r=x1r+2 local x1i=e*math_.sin(2*x.r) local m2=x2r*x2r+x1i*x1i local y={r=(x1i*x2r-x1r*x1i)/m2,i=-(x1r*x2r+x1i*x1i)/m2} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.tan(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then return math_.tan(z) elseif is.complex(x) then return complex.tan(x) else error("bad argument #1 to 'tan' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.tanh(x) -------------------------------------------------------------------------------- check.complex(x,1) local e=math_.exp(2*x.r) local x1r=e*math_.cos(2*x.i)-1 local x2r=x1r+2 local x1i=e*math_.sin(2*x.i) local m2=x2r*x2r+x1i*x1i local y={r=(x1r*x2r+x1i*x1i)/m2, i=(x1i*x2r-x1r*x1i)/m2} return setmetatable(y,complex) end -------------------------------------------------------------------------------- function math.tanh(x) -------------------------------------------------------------------------------- local z=tonumber(x) if z then return math_.tanh(z) elseif is.complex(x) then return complex.tanh(x) else error("bad argument #1 to 'tanh' (real or complex number expected, got %s)",type(x)) end end -------------------------------------------------------------------------------- function complex.__tostring(self) -------------------------------------------------------------------------------- if self.i==0 then return tostring(self.r) elseif self.r==0 then if self.i==1 then return "i" elseif self.i==-1 then return "-i" else return string.format("%si",tostring(self.i)) end else if self.i==1 then return string.format("%s+i",tostring(self.r)) elseif self.i==-1 then return string.format("%s-i",tostring(self.r)) elseif self.i>0 then return string.format("%s+%si",tostring(self.r),tostring(self.i)) else return string.format("%s%si",tostring(self.r),tostring(self.i)) end end end -------------------------------------------------------------------------------- function complex.__unm(x) -------------------------------------------------------------------------------- check.complex(x,1) return setmetatable({r=-x.r,i=-x.i},complex) end