I’ve been trying to implement/improve some functions in the math library and I’ve find it useful to use the free cloud version of Mathematica to generate test cases. I taught myself basics Mathematica/Wolfram language in the process and it is well possible I am missing something basic. Are there Stan users with Mathematica expertise that could double check that I did not mess up?

Tagging users that have talked about Mathematica here recently:

@vincent-picaud, @benwhite59, @sakrejda

The most polished case I have now is for `lbeta`

(logarithm of Beta function). The code is below, but can also be run from the notebook at: https://www.wolframcloud.com/obj/martin.modrak/Published/lbeta.nb

Here, I provide some additional comments to clarify what I intended. Overall I am generating C++ code that gets pasted into the test file. I want gold standard for both function values and derivatives w.r.t all arguments, which I hope Mathematica’s arbitrary precision math can give me.

```
lbeta[x_,y_]:= LogGamma[x] + LogGamma[y] - LogGamma[x + y]
dlbetadx[x_,y_]= D[lbeta[x,y],x];
dlbetady[x_,y_]= D[lbeta[x,y],y];
singleTest[x_,y_]:= Module[{val, cdx,cdy, big},{
val = N[lbeta[x,y],24];
(* Evaluating derivatives for large arguments exceeds memory, avoid it *)
cdx = If[x > 10^6 || y > 10^6,"NaN", CForm[N[dlbetadx[x,y],24]]];
cdy = If[x > 10^6 || y > 10^6,"NaN", CForm[N[dlbetady[x,y],24]]];
WriteString[out," {",CForm[x],",",CForm[y],",",
CForm[val],",",cdx,",",cdy,"},"]
}];
(* Wolfram Cloud limits the number of elements you can output, so writing to file. *)
out = OpenWrite["lbeta_test.txt"]
xs= SetPrecision[{0.00000001,0.004,1,1.00000001,5.0,23.0, 19845.0}, Infinity];
ys = SetPrecision[{0.00000000001,0.00002,1.000000000001,0.5,2.0,1624.0}, Infinity];
WriteString[out, "std::vector<TestValue> testValues = {"];
For[i = 1, i <= Length[xs], i++, {
For[j = 1, j <= Length[ys], j++, {
singleTest[xs[[i]], ys[[j]]];
}]
}]
WriteString[out,"};"];
Close[out];
(* I couldn't find a way to access the file on cloud, so I write it to output
which avoids most of the fragmentaiton*)
FilePrint[%]
```

The most tricky part is that I want to compute in arbitrary precision, but I need to start with a representable floating point number. This is to avoid false positives in test when the C++ implementation actually calculated correct high-precision output but the difference between the exact number and nearest representable float is big enough to be non-negligible. An example to illustrate the issue:

```
y = 1;
N[lbeta[1 + 10^-12,y], 22]
-9.999999999995000000000*^-13
N[lbeta[1.000000000001,y], 22]
-1.00009*^-12
N[lbeta[SetPrecision[1.000000000001,Infinity],y], 22]
-1.000088900581840922709*^-12
```

I now test for differences up to `10^-16`

so the first one triggered a false positive in the test.

Thanks for any input on the code.