Hello, on theirs paper "Fast, minimum storage ray/triangle intersection", the algorithm is presented as having two slight variations, one that can be used to cull back face triangles, postponing the use of the division operation and the other, that cannot do so and works for both kind of triangles.

On my Ray-Tracer, I have to use the one that postpone the divisions, intended to be used with front faced triangles, to calculate the position and use the one that cannot postpone, intended to be used with two sided triangles, to calculate the shadows. On that case I get perfect image as can be seen on this link: https://postimg.org/image/d8rrpu0fv/

Now if I invert the order (not postponed division to determine triangles position, postponed for shadows) I get no shadows and this light blue strip region:  https://postimg.org/image/6ibagdq4r/

If I use the postponed, front face in both cases I get lack of shadows: https://postimg.org/image/wqmf5rpnv/

And if a use the non postponed back faced version in both occasions I get the light blue strip: https://postimg.org/image/r2g4euy63/

Here is the code for the one that postpones the division operation, followed by the one that don't.

 Vector edge1 = s.vert2 - s.vert1;
   Vector edge2 = s.vert3 - s.vert1;
   Vector s1;
   s1.x = ray.dir.y * edge2.z - ray.dir.z * edge2.y;
   s1.y = ray.dir.z * edge2.x - ray.dir.x * edge2.z;
   s1.z = ray.dir.x * edge2.y - ray.dir.y * edge2.x;
   float det = edge1 * s1;
   if (det<0.000001)
      return false;
   Vector distance = ray.origin - s.vert1;
   float barycCoord_1 = distance * s1;
   if (barycCoord_1<0.0 || barycCoord_1>det)
      return false;
   Vector s2;
   s2.x = distance.y * edge1.z - distance.z * edge1.y;
   s2.y = distance.z * edge1.x - distance.x * edge1.z;
   s2.z = distance.x * edge1.y - distance.y * edge1.x;
   float barycCoord_2 = ray.dir * s2;
   if (barycCoord_2 < 0.0 || (barycCoord_1 + barycCoord_2) > det)
      return false;
   float intersection = edge2 * s2;
   float invDet = 1/det;
   intersection *= invDet;
   barycCoord_1 *= invDet;
   barycCoord_2 *= invDet;
   if (0.1f<intersection && intersection<t) {
      t=intersection;
      return true;
   }
   return false;

 Vector edge1 = s.vert2 - s.vert1;
   Vector edge2 = s.vert3 - s.vert1;
   Vec3f edge11 (s.vert22 - s.vert11);
   Vec3f edge22 (s.vert33 - s.vert11);
   Vector s1;
   s1.x = ray.dir.y * edge2.z - ray.dir.z * edge2.y;
   s1.y = ray.dir.z * edge2.x - ray.dir.x * edge2.z;
   s1.z = ray.dir.x * edge2.y - ray.dir.y * edge2.x;
   float divisor =  s1 * edge1;
   if (divisor==0.0)
      return false;
   float invDivisor = 1/divisor;
   Vector distance = ray.origin - s.vert1;
   float barycCoord_1 = distance * s1 * invDivisor;
   if (barycCoord_1 < 0.0 || barycCoord_1 > 1.0)
      return false;
   Vector s2;
   s2.x = distance.y * edge1.z - distance.z * edge1.y;
   s2.y = distance.z * edge1.x - distance.x * edge1.z;
   s2.z = distance.x * edge1.y - distance.y * edge1.x;
   float barycCoord_2 = ray.dir * s2 * invDivisor;
   if (barycCoord_2 < 0.0 || (barycCoord_1 + barycCoord_2) > 1.0)
      return false;
   float intersection = edge2 * s2 * invDivisor;
   if (0.1f<intersection && intersection<t) {
      t=intersection;
      return true;
   }
   return false;

I have no idea of what is going wrong and would appreciate suggestion about what can be going wrong, and if that is not possible, perhaps suggestions of what to do to try to figure out what is happening.

Thanks in advance.