I have a lattice three dimensional lattice of size lx * ly * lz with periodic boundary conditions on all three sides. My planar symmetry cuts are two horizontal and vertical planes (slope 0 and inf) plus two diagonal cuts (slope -1 and +1) perpendicular to the xy plane, the yz plane and the xz plane. So I have a total of 9 types of planes, all of them passing through lattice points and not between them.

Each of those cuts can pass through at any place on the plane perpendicular to them. For example a slope 0 cut perpendicular to the xy plane can pass at any of the ly points. A slope inf cut perpendicular to the xy plane can pass at any of the lx points. A slope +1 cut perpendicular to the xy plane can pass through any of the lx intersects. A slope -1 cut perpendicular to the xy plane can pass through any of the ly intersects.

Now given a site i on the lattice between 0 and (lx * ly * lz)-1, I would like to reflect it about any one of the random cuts c between 1 and 9. What would be the fastest algorithm to calculate the reflected lattice site preferably in 'C'.

The algorithm should take three inputs , the site i , the type of cut c, and the intersect through which the cut passes, between 0 and lx-1 or 0 and ly-1 or 0 and lz-1, depending on the cut and output the reflected site.

Writing plane as an equation

Instead of “slopes” I'd rather speak about normals of these planes. So the way I understand your question, you have

xy plane, slope  0 ⇒ normal (0,  1,  0)
xy plane, slope  ∞ ⇒ normal (1,  0,  0)
xy plane, slope  1 ⇒ normal (1, -1,  0)
xy plane, slope -1 ⇒ normal (1,  1,  0)
xz plane, slope  0 ⇒ normal (0,  0,  1)
xz plane, slope  ∞ ⇒ normal (1,  0,  0)
xz plane, slope  1 ⇒ normal (1,  0, -1)
xz plane, slope -1 ⇒ normal (1,  0,  1)
yz plane, slope  0 ⇒ normal (0,  0,  1)
yz plane, slope  ∞ ⇒ normal (0,  1,  0)
yz plane, slope  1 ⇒ normal (0,  1, -1)
yz plane, slope -1 ⇒ normal (0,  1,  1)

So the 9 kinds of planes would correspond to the normal directions

(1, 0, 0), (0, 1, 0), (0, 0, 1),
(1, 1, 0), (1, 0, 1), (0, 1, 1),
(1, -1, 0), (1, 0, -1), (0, 1, -1)

For each of these directions, you can take the normal vector (a, b, c) and turn that into the equation of a plane:

a*x + b*y + c*z = d

but what values for d are premissible? For the first row above, the planes parallel to one of the coordinate planes, things are simple: for (a, b, c) = (1, 0, 0) you have 0 ≤ d < lx, and similar for the other two. For the diagonal planes, your (in my opinion strange) intercept rules hold. If I understand you right, the (1, -1, 0) planes could pass through any point on the x axis, leading to 0 ≤ d < lx again. The (1, 1, 0) planes could pass through any point on the y axis, so you'd have 0 ≤ d < ly. For the other diagonals, please work out the bounds for d yourself.

Reflecting in such a plane

So now you have the equation of a plane, and want to reflect in that plane. The link Woodface provided is essentially the right thing to consider here, but you might prefer a simpler formulation of that idea. Start by rewriting the equation of the plane to

a*x + b*y + c*z - d = 0

If the left hand side is not zero, then a given point does not lie on the plane. In that case, the value you obtain is proportional to the distance of that point from the plane. Assume for the moment that a²+b²+c²=1. In that case, the value of the left hand side would indeed be the distance form the plane. Multiplying that number by the normal vector (a, b, c) you obtain a vector which points from the plane to the point in question. Multiplying the quantity by -(a, b, c) instead you get a vector pointing from the point to the plane, and multiplying by -2*(a, b, c) you get a vector pointing from the point to its mirror image.

But what if a²+b²+c²≠1? In that case, the value of the left hand side of the equation will be sqrt(a²+b²+c²) times the real distance. And you multiply that distance by a vector which has not unit length but length sqrt(a²+b²+c²) instead, so your final result vector would be too big by a factor of a²+b²+c² in total. So what you do is you scale by that factor.

To sum it up: to reflect a point (x, y, z) in the plane a*x + b*y + c*z = d you compute

(x, y, z) - 2/(a² + b² + c²)*(a*x + b*y + c*z - d)*(a, b, c)

or written in C code:

int f = 2/(a*a + b*b + c*c)*(a*x + b*y + c*z - d);
x = x - f*a;
y = y - f*b;
z = y - f*z;

You can use int here since for your normal vectors, a²+b²+c² will be either 1 or 2, so 2/(a*a + b*b + c*c) will always be an integer.