Constructing a straight line from two points not working in integer modulo group

I have a pair of coordinates of which all values belong to $Z_{p}$ where $p$ is a prime. I want to construct a straight line that goes through those two coordinates. Then I want to generate two more random points on that line. I also need a function that returns the value of y given any x on that line. For that I can treat it as a linear diophantine equation and generate any number of random integer coordinates on that line.

So first I find out the $a, b, c$ in the equation of $ax + by = c$ where $ a = (x_{2} - x_{1}) ,\quad b = (y_{1} - y_{2}),\quad c = (x_{2}y_{1} - y_{2}x_{1}) $ and then I compute the value of $y^{*} = f(x^{*}) = \frac{c - ax^{*}}{b}$. I compute $f(x_{1}) = y_{1}$ to check if everything is correct or not.

It all works out if I don't perform these operations in group. But when I do it in an integer modulo group group then I get $f(x_{1}) \neq y_{1}$. Surprisingly I am getting $c = 0$ even though $y_{2}x_{1} \neq x_{2}y_{1}$.

So, there is something different that has to be done in case of modular arithmetic. I checked another question that has a similar objective.

There fist $y = \frac{ax + c}{b}$ is formulated and then it says for $y$ to be an integer $ax + c \equiv 0 \; (mod\; b)$. If I translate it to here then it becomes $ax + c \equiv 0 \; (mod\; b)\; (mod\; p)$. But what does that mean and how to solve that ?

I am using Crypto++ and following is my code.

CryptoPP::Integer x1("59155865408868259425557571612471620428297988408165845346132791388203026665780027334884289790496698975366550859751679098576474914656477108493456088740714638392460074168696957786644452257687431017811952658833178670113578592426217762888790213090753444458649558440508900270736330461377938922007295905864081619510.");
CryptoPP::Integer x2("49171792422277804394667843819141913847322600280665849545551047558807915445878344561751437746234579327835202256618036154257179856084340066038426734660095847831493104314549773675359068536806511510632346258713769955285829778738241140307603773114365677760401574036178965775449693351366955171071039150496892512626..");
CryptoPP::Integer y1("15532588564541259482922018944097853897300168981677954188716754703307497158337714880433137687641278802841087340039578684036658054493398074873834397457078926775924304480658287770519079275439081927738422285768406410397661454394363161024640626036869355676305517960711933863623200507534787097131264154704514178528.");
CryptoPP::Integer y2("18681104435936101884448227316030762542107116083698320896806859540547443331269148757807169307897834606930640490118708410979463001754418260497415780748530755543566898480965940648878220684347582565134185256730690865433076765778856737642958127745233063853975442413445206113883902055488910897572165540719794965920.");
CryptoPP::Integer p ("171320888899847002522538286010717970050712035465289873249959090208041304926005539317533010051124825774622387763181990139801020807777709472681854339326432257665967930555817418179389068496976851304848858267411231278762526840486673505596715269737069851957355166043141770807768677357720176332759237252035377645721");

std::cout << "x1: " << x1 << std::endl << "x2: " << x2 << std::endl << "y1: " << y1 << std::endl << "y2: " << y2 << std::endl;
std::cout << std::endl;

assert(x1 < p);
assert(x2 < p);
assert(y1 < p);
assert(y2 < p);

{
    std::cout << "Group: " << std::endl;
    CryptoPP::ModularArithmetic Gp(p);

    CryptoPP::Integer dx  = Gp.Subtract(x2, x1), dy = Gp.Subtract(y2, y1), mdy = Gp.Subtract(y1, y2);
    CryptoPP::Integer c   = Gp.Subtract(Gp.Multiply(x2, y1), Gp.Multiply(y2, x1));

    std::cout << "Gp.Multiply(x2, y1): " << Gp.Multiply(x2, y1) << std::endl
              << "Gp.Multiply(y2, x1): " << Gp.Multiply(y2, x1) << std::endl
              << "c: " << c << std::endl
              << "Gp.Multiply(x2, y1) - Gp.Multiply(y2, x1): " << (Gp.Multiply(x2, y1) - Gp.Multiply(y2, x1)) << std::endl;

    CryptoPP::Integer a = mdy, b = dx;

    std::cout << Gp.Divide(Gp.Subtract(c, Gp.Multiply(x1, a)), b) << std::endl;
    std::cout << Gp.Add( Gp.Multiply(a, x1), Gp.Multiply(b, y1) ) << " " << c << std::endl;
}

{
    std::cout << "Not in Group: " << std::endl;
    CryptoPP::Integer dx = x2 - x1, dy = y2 - y1, mdy = y1 - y2;
    CryptoPP::Integer c  = (x2 * y1) - (y2 * x1);

    CryptoPP::Integer a = mdy, b = dx;

    std::cout << (c - (x1 * a)) / b << std::endl;

}

Define $x_\Delta=x_2-x_1$ and $y_\Delta=y_2-y_1$. Assume one of $x_\Delta\not\equiv0\pmod p$ or $y_\Delta\not\equiv0\pmod p$, so that the line is defined [or define that in this case solutions are the $p^2$ points $(x,y)\in\mathbb Z_p\times\mathbb Z_p$ ]

The equation of the line is $(x-x_1)y_\Delta\equiv(y-y_1)x_\Delta\pmod p$.

When $x_\Delta\equiv0\pmod p$, there are $p$ solutions $(x_1,y)$ with $y\in\mathbb Z_p$.

When $x_\Delta\not\equiv0\pmod p$, since $p$ is prime, that $x_\Delta$ has well-defined multiplicative inverse $x_\Delta^{-1}\bmod p$. By definition, that's the integer in $[0,p)$ such that $x_\Delta\,x_\Delta^{-1}\equiv1\pmod p$. Therefore there are $p$ solutions $\bigl(x,((x-x_1)\,y_\Delta\,x_\Delta^{-1}+y_1\bmod p)\bigr)$ with $x\in\mathbb Z_p$.

$x_\Delta^{-1}\bmod p$ can be computed by the (half) extended Euclidean algorithm (with a variant using only non-negative integers there), and a number of other methods. Crypto++ has MultiplicativeInverse. Alternatively, since $p$ is prime, $x_\Delta^{-1}\bmod p$ can be computed as ${x_\Delta}^{p-2}\bmod p$.