Compute direction of gravity vector
First, a distinction between &quot;gravity&quot; and &quot;gravitation&quot; as used in geophysics. &quot;Gravitation&quot; pertains to Newton's law of gravitation. (Nobody uses general relativity to model gravitation for small, lumpy masses such as the Earth, Mars, or the Moon.) &quot;Gravity&quot; pertains to how things appear to fall from the perspective of an observer fixed with respect to the rotating Earth. Thus gravity includes gravitational acceleration (inward, more or less toward the center of the Earth) and centrifugal acceleration (outward, away from the Earth's rotation axis). This question is asking about gravity rather than gravitation.

As I understand to compute vector of gravity not saficient to compute normal to elipsoid, but we need to computer normal to geoid (by definition geoid is a surface to which the force of gravity is everywhere perpendicular).

It's the other way around. The geoid is a calculated surface. Before the satellite era, one of the key sets of inputs to the calculation of the geoid were observations of deviations of local gravity, particularly the deflection of the vertical. Which way gravity pointed at several locations gave clues as to the local shape of the geoid.
Observing the orbits of satellites provides a global measure Earth's gravitational field. Satellites have been specially constructed for this job, most recently, GRACE and GOCE. The Earth's gravitational field is published in terms of spherical harmonics coefficients. The gravitational vector for a point in space on or above the Earth's surface can be calculated directly from these coefficients. Vectorially adding in the centrifugal acceleration due to the Earth's rotation results in the gravity vector. The position also yields the nominal gravity vector assuming an ellipsoidal Earth.

So I have some questions:

How to compute normal to geoid?


As noted above, a geoid is not needed. Modern models of the geoid are calculated from those same spherical harmonic coefficients (plus Earth rotation) that used needed to calculate the gravitational acceleration and gravity.


What is the maximum discrepancy of angle between elipsoid normal and geoid normal?


The technical term is &quot;deflection of the vertical&quot; (with variations). Per Hirt et al., it's up to 100 seconds of arc, about 10 kilometers south of the peak Annapurna II. This is a calculated value based on various satellite models (which are a bit coarse) coupled with digital terrain maps coupled with some more hairy mathematics to create fine-scale models.
C. Hirt, et al., “New ultrahigh‐resolution picture of Earth's gravity field,” Geophysical Research Letters, 40.16 (2013): 4279-4283.