Look at the limitations.
You can have X^N=998^998. That is 9,98e1000.
Look at fundamental data types range. double is something like not greater (-1e309;1e309) with 15 digits after point .  That is much smaller than you need.  So,  in function pow there is an overflow.  That means double LOSES PRECISION.  And if actually X^N mod M == Y for some big X^N them your program can't see it.  Because last digits of X^N are lost due to an overflow.

And in netman's code,  there is such a function,  that keeps ONLY the last digits.