Exact Normal Vectors On 3D Surfaces

[HazardHarry]

Hi everyone. I need a mathematician's help for an idea I had.
I've been casually writing a fractal generator for the past 5 years for fun, but lately I've been using it to teach myself 3D computer graphics:
http://giphy.com/gifs/3d-fractal-3oKIPpVsUFK3FcY65O/

I recently had an idea for a different kind of 3D fractal generator from the topological-style one I've been working with so far, but I've hit a bit of a snag at the beginning of the second iteration.

I have some ideas for other steps/variations I could add, but for now here's the basic gist:
(Forgive me if my math lingo is terrible. I'm doing my best. XD)

---
Given a Sphere (S0):

Iteration 1
For each Point (P) on S0:
1. Find the Unit Normal Vector (V) for P on S0.
2. Find the Complex number (Z) on the Riemann Sphere that has the Unit Normal Vector equal to V.
3. Perform some operation (e.g. Z^2 + u) to obtain a new Complex number (Z_Prime).
4. Perform some operation (e.g. Move [angle(Z, origin, x-axis)] units from P in the direction of V) to obtain a new Point (P_Prime).
The collection of P_Primes defines a new Shape (S1).

Iteration 2
For each Point (P) on S1:
1. ***Find the Unit Normal Vector (V) for P on S1.***
---

So here's where I'm stuck. I have no idea how to calculate the exact Unit Normal Vector on this new surface. I could estimate it by calculating extra points for S1, but it would be computationally intensive and not accurate. Accuracy is key because adjusting the vector a tiny bit could result in Z being off by a large margin if the change occurs with a vector that points towards the half of the Riemann Sphere that leads toward infinity.

If there isn't a way I'll proceed with the estimation method. Even if the end result is a bit wonky I still might make something cool and I'll definitely still be learning more about computer graphics. : )

[Tglad]

I think the technical terminology is that the point transformation is contravarient but the normal transformation is covariant.

In practice this means that any stretch in the spacing of neighbouring Ps is a compression for the normal on that axis. i.e. you use the inverse of the stretch of the transformation at each point (its Jacobian I think).

Simple case: if you transform a 1x1x1 sphere to a 1x2x3 ellipsoid then the corresponding normal is scaled by 1,1/2,1/3 on each axis.

[HazardHarry]

Thanks Tglad! I found a reference for what you were talking about and I believe it is exactly what I'm looking for.
http://math.etsu.edu/multicalc/prealpha/Chap3/Chap3-3/printversion.pdf

Man, I knew one day I'd regret not taking Calculus 3...