This is an extremely difficult non-linear integer programming problem. The difficulty consists in the restriction $h_1^{x_1}-h_2^{x_2}=kp$ which is hardly handled by known methods. I don't know any specialized soft to this end. Mathematica 14 easily solves it for

```
p = Prime[10^2]; SeedRandom[1234]; h1 = RandomChoice[Range[p - 1]];
h2 = RandomChoice[Range[p - 1]]; NMinimize[{(x1 - x2)^2,
Mod[h1^x1 - h2^x2, p] == 0 && x1 \[Element] Integers &&
x2 \[Element] Integers && 0 < x1 && x1 < p - 1 && x2 > 0 &&
x2 < p - 1}, {x1, x2}] // Timing
```

`{0.671875, {9., {x1 -> 192, x2 -> 195}}}`

with $p =541$.`SeedRandom[1234];`

guaranties the reproducibility of randomly chosen $h_1,h_2$ . Making use of options, one obtains

```
p = Prime[2*10^2]; SeedRandom[1234]; h1 = RandomChoice[Range[p - 1]];
h2 = RandomChoice[Range[p - 1]]; NMinimize[{(x1 - x2)^2,
Mod[h1^x1 - h2^x2, p] == 0 && x1 \[Element] Integers &&
x2 \[Element] Integers && 0 < x1 && x1 < p - 1 && x2 > 0 &&
x2 < p - 1}, {x1, x2}, Method -> {"DifferentialEvolution", "ScalingFactor" -> 1,
"SearchPoints" -> 150}, MaxIterations -> 2000] // Timing
```

`{133.375, {23409., {x1 -> 936, x2 -> 783}}}`

with $p=1223$. My comp is not powerful and fails with `p=Prime[10^3]`

(equals $7919$). The methods are described in the documentation.

4more comments