st: Precision issues with Mata. Where's my error?

[Marta García-Granero <mgarciagranero@gmail.com>]

Hi everybody:

I have burned my eyelashes trying to spot my error, but I can't find it.

Situation: I am translating to Stata/mata an old SPSS/MATRIX code, that 
I succeeded in translating into an EXCEL VBA-UDF. It computes confidence 
limits for Newcombe's method 10 for paired proportions (I have found and 
installed an ado for independent samples, but I could not find anything 
for paired samples, I'd be glad if anyone had something to share).

This is the original Excel VBA function (really close to mata syntax, 
that's why I thought the translation would be a piece of cake)

' R.G. Newcombe. Stats in Medicine 1998, 17, 2635-2650. Method 10
Public Function Newcombe10p(posAposB As Integer, posAnegB As Integer, 
negAposB As Integer, negAnegB As Integer, alfa As Double)
Dim res(1 To 3) As Variant
Z = Application.WorksheetFunction.NormSInv(1 - alfa / 2)
ZSquare = Z ^ 2
x5 = posAposB + posAnegB
x6 = negAposB + negAnegB
x7 = posAposB + negAposB
x8 = posAnegB + negAnegB
x9 = x7 + x8
x10 = (posAnegB - negAposB) / x9
x11 = 2 * x5 + ZSquare
x12 = Z * (ZSquare + 4 * x5 * x6 / x9) ^ 0.5
x13 = 2 * (x9 + ZSquare)
x14 = (x11 + x12) / x13
x15 = (x11 - x12) / x13
x16 = x5 / x9 - x15
x17 = x14 - x5 / x9
x21 = 2 * x7 + ZSquare
x22 = Z * (ZSquare + 4 * x7 * x8 / x9) ^ 0.5
x24 = (x21 + x22) / x13
x25 = (x21 - x22) / x13
x26 = x7 / x9 - x25
x27 = x24 - x7 / x9
x29 = x5 * x6 * x7 * x8
x30 = 1
If (x29 = 0) Then x30 = 0
x31 = posAposB * negAnegB - posAnegB * negAposB
x32 = 0
If (x31 > 0) Then x32 = 1
x33 = x31 - x9 / 2
x35 = Application.WorksheetFunction.max(x33, 0)
x36 = x32 * x35 + (1 - x32) * x31
x37 = x30 * x36
x38 = x30 * x29 ^ 0.5 + (1 - x30)
x39 = x37 / x38
x40 = x16 * x16 - 2 * x39 * x16 * x27 + x27 * x27
x41 = x17 * x17 - 2 * x39 * x17 * x26 + x26 * x26
x42 = x10 - (x40) ^ 0.5
x43 = x10 + (x41) ^ 0.5
res(1) = x10
res(2) = x42
res(3) = x43
Newcombe10p = res
End Function

And this is my mata translation (still very preliminary, in due time, 
when I acquire certain knowledge, I'll make an ado like: pairpcii 16 13 
1 20, level(95) ).

mata
// Input data (modified by user)
papb=16
panb=13
napb=1
nanb=20
alpha=0.05
// Calculations start here
z=invnormal(1-alpha/2)
zsq=z^2
x5=papb+panb
x6=napb+nanb
x7=papb+napb
x8=panb+nanb
x9=x7+x8
x10=(panb-napb)/x9
x11=2*x5+zsq
x12=z*sqrt(zsq+4*x5*x6/x9)
x13=2*(x9+zsq)
x14=(x11+x12)/x13
x15=(x11-x12)/x13
x16=x5/x9 - x15
x17=x14-x5/x9
x21=2*x7+zsq
x22=z*sqrt(zsq+4*x7*x8/x9)
x24=(x21+x22)/x13
x25=(x21-x22)/x13
x26=x7/x9-x25
x27=x24-x7/x9
x29=x5*x6*x7*x8
x30=(x29==0 ? 0 : 1)
x31=papb*nanb-panb*napb
x32=(x31>0 ? 1 :0)
x33=x31-x9/2
x35=(x33>0 ? x33 : 0)
x36=x32*x35+(1-x32)*x31
x37=x30*x36
x38=x30*sqrt(x29)+(1-x30)
x39=x37/x38
x40=x16^2-2*x39*x16*(x27^3)
x41=x17^2-2*x39*x17*(x26^3)
x42=x10-sqrt(x40)
x43=x10+sqrt(x41)
printf("%6.2f", 100*x10)
printf("%6.2f", 100*x42)
printf("%6.2f", 100*x43)
end

The problem is that I expected x42 and x43 to be 9.95 and 36.34 (checked 
with other software, like Winpepi, or SPSS 21 NPTESTS command), and I 
get 10.36 and 36.55.

I can't find my error.

Any hints?

Thanks a lot in advance,
Marta GG
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/