Point And Complex

wiki

Points and Complex numbers in Lua

It is not difficult to redefine the meaning of the usual arithmetic operators for user-defined types. Normally, this will be a run-time error; if I have a table t then t + 1 will give "attempt to perform arithmetic on global `t' (a table value)". However, if that table has a metatable, then Lua will see if a function __add is defined, and use it instead.

The main question with operator overloading is whether it will make using your objects easier for other programmers. It is similar to user-interface decisions; it really is a good idea to keep to common user-interface conventions for your platform. In interfaces, dull is good, and the unexpected will cause misunderstanding. For instance, you can make the addition operator concatenate lists, but as Paul Graham observed [1], one tends to misread such expressions as arithmetic.

I'm going to present two cases where redefining operators makes perfect sense, because they are both generalizations of the usual real number arithmetic we use. To add two points with p1 + p2 is not only convenient, but mathematically correct.

Using this Point class, you will be able to express vector algebra in something like the usual notation. For example,

x = ((p1^p2)..q)*q

-- point.lua
-- A class representing vectors in 3D
-- (for class.lua, see SimpleLuaClasses)
require 'class'

Point = class(function(pt,x,y,z)
   pt:set(x,y,z)
 end)

local function eq(x,y)
  return x == y
end

function Point.__eq(p1,p2)
  return eq(p1[1],p2[1]) and eq(p1[2],p2[2]) and eq(p1[3],p2[3])
end

function Point.get(p)
  return p[1],p[2],p[3]
end

-- vector addition is '+','-'
function Point.__add(p1,p2)
  return Point(p1[1]+p2[1], p1[2]+p2[2], p1[3]+p2[3])
end

function Point.__sub(p1,p2)
  return Point(p1[1]-p2[1], p1[2]-p2[2], p1[3]-p2[3])
end

-- unitary minus  (e.g in the expression f(-p))
function Point.__unm(p)
  return Point(-p[1], -p[2], -p[3])
end

-- scalar multiplication and division is '*' and '/' respectively
function Point.__mul(s,p)
  return Point( s*p[1], s*p[2], s*p[3] )
end

function Point.__div(p,s)
  return Point( p[1]/s, p[2]/s, p[3]/s )
end

-- dot product is '..'
function Point.__concat(p1,p2)
  return p1[1]*p2[1] + p1[2]*p2[2] + p1[3]*p2[3]
end

-- cross product is '^'
function Point.__pow(p1,p2)
   return Point(
     p1[2]*p2[3] - p1[3]*p2[2],
     p1[3]*p2[1] - p1[1]*p2[3],
     p1[1]*p2[2] - p1[2]*p2[1]
   )
end

function Point.normalize(p)
  local l = p:len()
  p[1] = p[1]/l
  p[2] = p[2]/l
  p[3] = p[3]/l
end

function Point.set(pt,x,y,z)
  if type(x) == 'table' and getmetatable(x) == Point then
     local po = x
     x = po[1]
     y = po[2]
     z = po[3]
  end
  pt[1] = x
  pt[2] = y
  pt[3] = z 
end

function Point.translate(pt,x,y,z)
   pt[1] = pt[1] + x
   pt[2] = pt[2] + y
   pt[3] = pt[3] + z 
end

function Point.__tostring(p)
  return string.format('(%f,%f,%f)',p[1],p[2],p[3])
end

local function sqr(x) return x*x end

function Point.len(p)
  return math.sqrt(sqr(p[1]) + sqr(p[2]) + sqr(p[3]))
end

Point is a simple class which can be constructed using call notation (see SimpleLuaClasses):

> p1 = Point(10,20,30)
> p2 = Point(1,2,3)
> = p1
(10.000000,20.000000,30.000000)
> = p1 + p2
(11.000000,22.000000,33.000000)
> = 2*p1
(20.000000,40.000000,60.000000)

One can now define higher level operations. For instance, here is a useful function for finding the minimum distance between a point and a line (very useful if you're doing an editable graphics program):

-- given a point q, where does the perp cross the line (p1,p2)?
function perp_to_line(p1,p2,q)
  local diff = p2 - p1
  local x = ((q - p1)..diff)/(diff..diff)
  return p1 + x*diff
end

-- minimum distance between q and a line (p1,p2)
function min_dist_to_line(p1,p2,q)
  local perp = perp_to_line(p1,p2,q)
  return Point.len(perp-q)
end

There is a 'gotcha' which you should keep in mind; Lua objects are always passed by reference, so beware of modifying points passed to functions. Use the copy constructor provided, e.g. say local pc = Point(p) to make a local copy of a point passed as an argument.

Complex Numbers

Complex numbers are a generalization of real numbers, so they understand all the usual operations, plus some more. If z is complex, then 1.5 + z and z + 1.5 are both complex expressions, so __add must handle the case where one of the arguments is a plain number.

Please note that this Complex Number class is an example, it should not be used for serious applications. It has some problems (which could, in principle, be fixed): division suffers from loss of precision, modulus overflows on some reasonable values, square root is careless with cuts and does not work for real values, pow only computes positive integer powers (and then slowly).

-- complex.lua
require 'class'

Complex = class(function(c,re,im)
                 if type(re) == 'number' then 
                   c.re = re
                   c.im = im
                 else
                   c.re = re.re
                   c.im = re.im
                 end
          end)

Complex.i = Complex(0,1)

local sqrt = math.sqrt
local cos = math.cos
local sin = math.sin
local exp = math.exp

local function check(z1,z2)
  if     type(z1) == 'number' then return Complex(z1,0),z2
  elseif type(z2) == 'number' then return z1,Complex(z2,0) 
  else return z1,z2
  end
end

-- redefine arithmetic operators!
function Complex.__add(z1,z2)
  local c1,c2 = check(z1,z2)
  return Complex(c1.re + c2.re, c1.im + c2.im)
end

function Complex.__sub(z1,z2)
  local c1,c2 = check(z1,z2)
  return Complex(c1.re - c2.re, c1.im - c2.im)
end

function Complex:__unm()
  return Complex(-self.re, -self.im)
end

function Complex.__mul(z1,z2)
  local c1,c2 = check(z1,z2)
  return Complex(c1.re*c2.re - c1.im*c2.im, c1.im*c2.re + c1.re*c2.im)
end

function Complex.__div(z1,z2)
  local c1,c2 = check(z1,z2)
  local a = c1.re
  local b = c1.im
  local c = c2.re
  local d = c2.im
  local lensq = c*c + d*d
  local ci = (a*c + b*d)/lensq
  local cr = (b*c + a*d)/lensq
  return Complex(cr,ci)
end

function Complex.__pow(z,n)
  local res = Complex(z)
  for i = 1,n-1 do res = res*z end
  return res  
end  

-- this is how our complex numbers will present themselves!
function Complex:__tostring()
  return self.re..' + '..self.im..'i'
end

-- operations only valid for complex numbers
function Complex.conj(z)
  return Complex(z.re,-z.im)
end

function Complex.mod(z)
  return sqrt(z.re^2 + z.im^2)
end

-- generalizations of sqrt() and exp()
function Complex.sqrt(z)
  local y = sqrt((Complex.mod(z)-z.re)/2)
  local x = z.im/(2*y)
  return Complex(x,y)
end

function Complex.exp(z)
  return exp(z.re)*Complex(cos(z.im),sin(z.im))
end

But is it Fast?

Obviously vector and complex arithmetic is going to be fairly slow in Lua, but that's relative. I can create 100,000 points in about two seconds; same time for adding that number. If your program works with a few thousand points, it's going to be fast enough. Programs that need to work with much larger number of points (such as GIS systems) would simply not use this representation, even in C++. Instead, one could write a Lua extension for arrays of points, and manipulate them in groups.

Comments

The expression "x = ((p1^p2)..q)*q" was not at first obvious to me. When once sees "^" and ".." they think exponentiation and concatenation, both of which don't apply here. Exponentiation does not apply to vectors of length >= 2, so maybe that reduces the chance of confusion, but it does apply to square matrices, which are not much different from vectors and are often used in the same expression, such as A^(-1)Av. ".." is odd too, though I admit the two dots somewhat resemble dot product. Lua's automatic string conversion might increase the change for error--if one mistakenly writes "(p1..p2)..3" the ".." is silently taken as regular string concatenation. "*" of course could be understood as scalar product, dot product, or cross product. Unfortunately, the number of operators we have to use is limited. These issues detract from the original argument of the article, but I agree that "+", "-", and such make perfect sense for vectors. --DavidManura

I cannot agree with you, David. Anybody who has done maths up to first year undergraduate level would recognize '^' as vector product for 3-dimensional vectors, or as exterior product more generally. The '..' for inner product is unfortunate, I agree. Alas, the historical quirks of mathematical notation really do not fit well with computers. The best programming language I have come across for accommodating mathematical notation is Gofer (not Hugs or Haskell proper - the standardized prelude screws up, because it was designed by people without sufficient mathematical background). -- GavinWraith

I think "^" (or actually \(\wedge\)) is less common. I and someone else in the US hadn't seen it used this way, but I asked someone in the UK who said usage of "^" and "x" at the ugrad level was about fifty-fifty, though after learning exterior calculus he preferred "x" for cross product to keep the two ideas separate, something about them essentially being related by Hodge star. Here's another reference that applies "^" as cross-product though [2]. In some contexts, being able to have user-defined operators would be useful, mainly for the ability to express them in infix notation [3]. See also CustomOperators. --DavidManura

A potential area of confusion is that ^ is right-associative in Lua but cross-products are not associative. So, a^b^c would mean a^(b^c). --DavidManura

I stumbled onto this page in an unrelated search, but seeing the sentence "To add two points with p1 + p2 is not only convenient, but mathematically correct." just about made me cry. There is absolutely no geometrical meaning to adding two points. After skimming the article it was fairly clear that the 'Point' class is also intended to be used for vectors, for some reason? I'm not to change anything here, but I thought it should at least be pointed out. --anonymous

"Vector" would likely be more suitable under the given operations (e.g. normalize and cross product). However, the interchange in terminology above was likely because points can be represented as [position vectors]. --DavidManura