I think the R implementations of the cdfs on the log scale are more stable, but all their code’s GPL-ed and I’ve been afraid to use it. I don’t think you could be behind Stan. Stan doesn’t have an mvnormal or even binomral cdf, nor does it have any implementations of which I’m aware on the log scale. And absolutely none that derive a good algorithm for the derivatives beyond autodiffing the algorithm for the values. I think a lot of the expertise in writing special functions died out decades ago in the applied math community.