(*Reciprocal square root with a quadratic approximation followed by a 2nd-order Householder iteration.*)

fi[x_]:=pi.NestList[#*x&,1,Length[pi]-1]

f[x_]:=1/Sqrt[(x+1)/2]

xi = N[(Table[-Cos[i*π/3], {i, 0,3}]*3 + 1)/4]

pi={p0,p1,p2}/.(Solve[({p0,p1,p2}.{1,#,#^2}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ}*(f/@xi),{p0,p1,p2,ϵ}])[[1]]

{1.437799046117536`,-0.823394375837328`,0.4096419668459485`}

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/f[x],x]/D[(fi[x]-f[x])/f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x]- f[x])/f[x],{x,-0.5,1}]

Round[pi*16384]

mq14[a_,b_]:=Floor[a*b*(1./2^14)]

mq15[a_,b_]:=Floor[a*b*(1./2^15)]

rsqrt[x_]:=Module[{n,r,y},
n=x-32768;
r=23557+mq15[n,-13490+mq15[n,6713]];
y=2*(mq15[mq15[r,r],n]+mq15[r,r]-16384);
r+mq15[r,mq15[y,-16384+mq15[y,12288]]]
]

rsqrt/@{16384,16385,16386,16387,16388,16389,16390,16391}

rsqrt/@{65535,65534,65533,65532,65531,65530,65529,65528}

Plot[(rsqrt[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536]),{x,16384,65535}]

Sqrt[Mean[Table[N[(rsqrt[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]^2,{x,16384,65535}]]]

Max[Table[Abs[N[(rsqrt[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]],{x,16384,65535}]]

Max[Table[Abs[N[rsqrt[x]-16384/Sqrt[x/65536]]],{x,16384,65535}]]

(*Reciprocal square root with a linear approximation followed by two Newton iterations. This is the same amount of work, but is not as accurate as the above.*)

xi = N[(Table[-Cos[i*π/2], {i, 0,2}]*3 + 1)/4]

pi={p0,p1}/.(Solve[({p0,p1}.{1,#}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ}*(f/@xi),{p0,p1,ϵ}])[[1]]

{1.5234089594072642`,-0.6093635837629056`}

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/f[x],x]/D[(fi[x]-f[x])/f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x]- f[x])/f[x],{x,-0.5,1}]

Round[pi*16384]

rsqrt[x_]:=Module[{n,r,y},
n=x-32768;
r=24960+mq15[n,-9985];
y=mq15[mq15[r,r],n]+mq15[r,r]-16384;
r=r-mq15[r,y];
y=mq15[mq15[r,r],n]+mq15[r,r]-16384;
r=r-mq15[r,y]
]

rsqrt/@{16384,16385,16386,16387,16388,16389,16390,16391}

rsqrt/@{65535,65534,65533,65532,65531,65530,65529,65528}

Plot[(rsqrt[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536]),{x,16384,65535}]

Sqrt[Mean[Table[N[(rsqrt[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]^2,{x,16384,65535}]]]

Max[Table[Abs[N[(rsqrt[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]],{x,16384,65535}]]

Max[Table[Abs[N[rsqrt[x]-16384/Sqrt[x/65536]]],{x,16384,65535}]]

(*Reciprocal square root with pure polynomial approximations. These can be less work, but are nowhere near as accurate, even when using the same number of mulitplies as above.*)

xi = N[(Table[-Cos[i*π/5], {i, 0,5}]*3 + 1)/4]

pi={p0,p1,p2,p3,p4}/.(Solve[({p0,p1,p2,p3,p4}.{1,#,#^2,#^3,#^4}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ,-ϵ}*(f/@xi),{p0,p1,p2,p3,p4,ϵ}])[[1]]

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/(f[x]),x]/D[(fi[x]-f[x])/(f[x]),x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/f[x], {x, -0.5, 1}]

Round[16384*pi]

rsqrtp[x_]:=Module[{n},
n=x-32768;
23122+mq15[n,-11508+mq15[n,9847+mq15[n,-8939+mq15[n,3896]]]]
]

Plot[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536]),{x,16384,65535}]

Sqrt[Mean[Table[N[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]^2,{x,16384,65535}]]]

Max[Table[Abs[N[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]],{x,16384,65535}]]

Max[Table[Abs[N[rsqrtp[x]-16384/Sqrt[x/65536]]],{x,16384,65535}]]

xi = N[(Table[-Cos[i*π/6], {i, 0,6}]*3 + 1)/4]

pi={p0,p1,p2,p3,p4,p5}/.(Solve[({p0,p1,p2,p3,p4,p5}.{1,#,#^2,#^3,#^4,#^5}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ,-ϵ,ϵ}*(f/@xi),{p0,p1,p2,p3,p4,p5,ϵ}])[[1]]

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/f[x],x]/D[(fi[x]-f[x])/f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/f[x], {x, -0.5, 1}]

Round[16384*pi]

rsqrtp[x_]:=Module[{n},
n=x-32768;
23168+mq15[n,-11458+mq15[n,8685+mq15[n,-8731+mq15[n,7832+mq15[n,-3123]]]]]
]

Plot[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536]),{x,16384,65535}]

Sqrt[Mean[Table[N[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]^2,{x,16384,65535}]]]

Max[Table[Abs[N[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]],{x,16384,65535}]]

Max[Table[Abs[N[rsqrtp[x]-16384/Sqrt[x/65536]]],{x,16384,65535}]]

xi = N[(Table[-Cos[i*π/7], {i, 0, 7}]*3 + 1)/4]

pi={p0,p1,p2,p3,p4,p5,p6} /.(Solve[({p0,p1,p2,p3,p4,p5,p6}.{1,#,#^2,#^3,#^4,#^5,#^6}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ,-ϵ,ϵ,-ϵ}*(f/@xi),{p0,p1,p2,p3,p4,p5,p6,ϵ}])[[1]]

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]}, #-((D[(fi[x]-f[x])/f[x],x]/D[(fi[x]-f[x])/f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x]- f[x])/f[x],{x,-0.5,1}]

Round[16384*pi]

rsqrtp[x_]:=Module[{n},
n=x-32768;
23175+mq15[n,-11565+mq15[n,8473+mq15[n,-7378+mq15[n,8103+mq15[n,-6969+mq15[n,2549]]]]]]
]

Plot[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536]),{x,16384,65535}]

Sqrt[Mean[Table[N[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]^2,{x,16384,65535}]]]

Max[Table[Abs[N[(rsqrtp[x]-16384/Sqrt[x/65536])/(16384/Sqrt[x/65536])]],{x,16384,65535}]]

Max[Table[Abs[N[rsqrtp[x]-16384/Sqrt[x/65536]]],{x,16384,65535}]]

(*Polynomial approximation of a square root.*)

f[x_]:=Sqrt[(x+1)/2]

xi = N[(Table[-Cos[i*π/5], {i, 0,5}]*3 + 1)/4]

pi={p0,p1,p2,p3,p4}/.(Solve[({p0,p1,p2,p3,p4}.{1,#,#^2,#^3,#^4}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ,-ϵ}*(f/@xi),{p0,p1,p2,p3,p4,ϵ}])[[1]]

{0.7072403335147993`,0.35278823789299407`,-0.09189057085305771`,0.051860055091630206`,-0.020249053852030158`}

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/(f[x]),x]/D[(fi[x]-f[x])/(f[x]),x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/f[x], {x, -0.5, 1}]

Round[32768*pi]

sqrtp[x_]:=Module[{n},
n=x-32768;
23175+mq15[n,11561+mq15[n,-3011+mq15[n,1699+mq15[n,-664]]]]
]

Plot[(sqrtp[x]-Sqrt[x/65536]*32768)/(Sqrt[x/65536]*32768),{x,16384,65535}]

Sqrt[Mean[Table[N[(sqrtp[x]-Sqrt[x/65536]*32768)/(Sqrt[x/65536]*32768)]^2,{x,16384,65535}]]]

Max[Table[Abs[N[(sqrtp[x]-Sqrt[x/65536]*32768)/(Sqrt[x/65536]*32768)]],{x,16384,65535}]]

Max[Table[Abs[N[(sqrtp[x]-Sqrt[x/65536]*32768)]],{x,16384,65535}]]

(*JM's original square root approximation for comparison.*)

sqrtp[x_]:=Module[{n},
n=x-32768;
23174+mq15[n,11584+mq15[n,-3011+mq15[n,1570+mq15[n,-557]]]]
]

Plot[(sqrtp[x]-Sqrt[x/65536]*32768)/(Sqrt[x/65536]*32768),{x,16384,65535}]

Sqrt[Mean[Table[N[(sqrtp[x]-Sqrt[x/65536]*32768)/(Sqrt[x/65536]*32768)]^2,{x,16384,65535}]]]

Max[Table[Abs[N[(sqrtp[x]-Sqrt[x/65536]*32768)/(Sqrt[x/65536]*32768)]],{x,16384,65535}]]

Max[Table[Abs[N[(sqrtp[x]-Sqrt[x/65536]*32768)]],{x,16384,65535}]]

(*Polynomial approximation of cosine.*)

f[x_]:=Cos[Sqrt[x]*π/2]

xi = N[Table[-Cos[i*π/4], {i, 0,4}]+1]/2

pi={p0,p1,p2,p3}/.(Solve[({p0,p1,p2,p3}.{1,#,#^2,#^3}&/@xi)-(f/@xi)=={-ϵ,ϵ,-ϵ,ϵ,-ϵ},{p0,p1,p2,p3,ϵ}])[[1]]

{0.9999932952821674`,-1.233484503785825`,0.2525802391345619`,-0.01909573534873664`}

(fi/@xi)-(f/@xi)

xi=Join[{xi[[1]]},#-((D[fi[x]-f[x],x]/D[fi[x]-f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[fi[x] - f[x], {x, 0, 1}]

Round[32768*pi]+{0,32768,0,0}

mq15p[a_,b_]:=Round[a*b*(1./2^15)]

cospi2p[x_]:=Module[{x2},
x2=Max[mq15p[x,x]-1,0];
32767-x2+mq15[x2,-7650+mq15[x2,8277+mq15[x2,-626]]]
]

Plot[(cospi2p[x]-Cos[x*π/65536]*32768),{x,0,32767}]

Sqrt[Mean[Table[N[(cospi2p[x]-Cos[x*π/65536]*32768)]^2,{x,0,32767}]]]

Max[Table[Abs[N[(cospi2p[x]-Cos[x*π/65536]*32768)]],{x,0,32767}]]

(*This is JM's original approximation, which was pretty good, if slightly more complicated.*)

cospi2p[x_]:=Module[{x2},
x2=mq15p[x,x];
1+Min[32766,32767-x2+mq15p[x2,-7651+mq15p[x2,8277+mq15p[x2,-626]]]]
]

Plot[(cospi2p[x]-Cos[x*π/65536]*32768),{x,0,32767}]

Sqrt[Mean[Table[N[(cospi2p[x]-Cos[x*π/65536]*32768)]^2,{x,0,32767}]]]

Max[Table[Abs[N[(cospi2p[x]-Cos[x*π/65536]*32768)]],{x,0,32767}]]

(*Polynomial approximation of the binary logarithm.*)

f[x_]:=Log[2,(2*x+3)/4]

(*8-bit approximation*)

xi = N[(Table[-Cos[i*π/4], {i, 0,4}])/2]

pi={p0,p1,p2,p3}/.(Solve[({p0,p1,p2,p3}.{1,#,#^2,#^3}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ}*(2+f/@xi),{p0,p1,p2,p3,ϵ}])[[1]]

{-0.4144541824871411`,0.9590923197873218`,-0.3395129038105771`,0.16541096501128538`}

((fi/@xi)-(f/@xi))/(2+(f/@xi))

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/(2+f[x]),x]/D[(fi[x]-f[x])/(2+f[x]),x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/(2+f[x]), {x, -0.5, 0.5}]

Round[pi*{16384,8192,4096,2048}]

Round[pi*16384]

log2p[x_]:=Module[{n},
n=x-32768-16384;
Floor[(-6758+mq15[n,15715+mq15[n,-5563+mq15[n,2708]]])/64]
]

Plot[(log2p[x] - Log[2, x/65536]*256)/((2 + Log[2, x/65536])*256), {x,
   32768, 65535}]

Sqrt[Mean[
  Table[N[((log2p[x] - 
         Log[2, x/65536]*256)/((2 + Log[2, x/65536])*256))^2], {x, 
    32768, 65535}]]]

Max[Table[
  Abs[N[(log2p[x] - 
       Log[2, x/65536]*256)/((2 + Log[2, x/65536])*256)]], {x, 32768, 
   65535}]]

Max[Table[Abs[N[log2p[x] - Log[2, x/65536]*256]], {x, 32768, 65535}]]

(*10-bit approximation*)

xi = N[(Table[-Cos[i*π/5], {i, 0,5}])/2]

pi={p0,p1,p2,p3,p4}/.(Solve[({p0,p1,p2,p3,p4}.{1,#,#^2,#^3,#^4}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ,-ϵ}*(2+f/@xi),{p0,p1,p2,p3,p4,ϵ}])[[1]]

{-0.41509302963303146`,0.9609890551383969`,-0.31836011537636605`,0.15530808010959576`,-0.08556153059057618`}

((fi/@xi)-(f/@xi))/(2+(f/@xi))

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/(2+f[x]),x]/D[(fi[x]-f[x])/(2+f[x]),x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/(2+f[x]), {x, -0.5, 0.5}]

Round[pi*{16384,8192,4096,2048,1024}]

Round[pi*16384]

log2p[x_]:=Module[{n},
n=x-32768-16384;
Floor[(-6793+mq15[n,15746+mq15[n,-5217+mq15[n,2545+mq15[n,-1401]]]])/16]
]

Plot[(log2p[x]-Log[2,x/65536]*1024)/((2+Log[2,x/65536])*1024),{x,32768,65535}]

Sqrt[Mean[Table[N[((log2p[x]-Log[2,x/65536]*1024)/((2+Log[2,x/65536])*1024))^2],{x,32768,65535}]]]

Max[Table[Abs[N[(log2p[x]-Log[2,x/65536]*1024)/((2+Log[2,x/65536])*1024)]],{x,32768,65535}]]

Max[Table[Abs[N[log2p[x]-Log[2,x/65536]*1024]],{x,32768,65535}]]

(*This was JM's original approximation, which was pretty good except for the lack of rounding bias at the end.*)

log2p[x_]:=Module[{n},
n=x-32768-16384;
Floor[(-6791+mq14[n,7872+mq14[n,-1392+mq14[n,319]]])/64]
]

Plot[(log2p[x]-Log[2,x/65536]*256)/((2+Log[2,x/65536])*256),{x,32768,65535}]

Sqrt[Mean[Table[N[((log2p[x]-Log[2,x/65536]*256)/((2+Log[2,x/65536])*256))^2],{x,32768,65535}]]]

Max[Table[Abs[N[(log2p[x]-Log[2,x/65536]*256)/((2+Log[2,x/65536])*256)]],{x,32768,65535}]]

Max[Table[Abs[N[log2p[x]-Log[2,x/65536]*256]],{x,32768,65535}]]

(*Polynomial approximation of the binary exponential.*)

f[x_]:=2^x

xi = N[(Table[-Cos[i*π/4], {i, 0,4}])+1]/2

pi={p0,p1,p2,p3}/.(Solve[({p0,p1,p2,p3}.{1,#,#^2,#^3}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ}*(f/@xi),{p0,p1,p2,p3,ϵ}])[[1]]

{0.9999252185627103`,0.6958335404948245`,0.2260671554272475`,0.07802452264063828`}

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/f[x],x]/D[(fi[x]-f[x])/f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/f[x], {x, 0,1}]

Round[pi*16384]

Round[pi*{16384,32768,65536,65536*2}]

exp2p[x_]:=
16383+mq15[x,22804+mq15[x,14819+mq15[x,10204]]]

Plot[(exp2p[x]-2^(x/16384)*16384)/(2^(x/16384)*16384),{x,0,16383}]

Sqrt[Mean[Table[N[((exp2p[x]-2^(x/16384)*16384)/(2^(x/16384)*16384))^2],{x,0,16383}]]]

Max[Table[Abs[N[(exp2p[x]-2^(x/16384)*16384)/(2^(x/16384)*16384)]],{x,0,16383}]]

Max[Table[Abs[N[exp2p[x]-2^(x/16384)*16384]],{x,0,16384}]]

(*This was JM's original approximation.*)

exp2p[x_]:=
16384+mq14[x,11356+mq14[x,3726+mq14[x,1301]]]

Plot[(exp2p[x]-2^(x/16384)*16384)/(2^(x/16384)*16384),{x,0,16383}]

Sqrt[Mean[Table[Abs[N[((exp2p[x]-2^(x/16384)*16384)/(2^(x/16384)*16384))^2]],{x,0,16383}]]]

Max[Table[Abs[N[(exp2p[x]-2^(x/16384)*16384)/(2^(x/16384)*16384)]],{x,0,16383}]]

Max[Table[Abs[N[exp2p[x]-2^(x/16384)*16384]],{x,0,16384}]]

(*Polynomial approximation of the reciprocal.*)

f[x_]:=4/(x+3)

xi = N[(Table[-Cos[i*π/5], {i, 0,5}])]

pi={p0,p1,p2,p3,p4}/.(Solve[({p0,p1,p2,p3,p4}.{1,#,#^2,#^3,#^4}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ,-ϵ,ϵ,-ϵ}*(f/@xi),{p0,p1,p2,p3,p4,ϵ}])[[1]]

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/f[x],x]/D[(fi[x]-f[x])/f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/f[x], {x, -1,1}]

Round[pi*16384]

rcpp[x_]:=Module[{n},
n=x-65536-32768;
21845+mq15[n,-7249+mq15[n,2416+mq15[n,-935+mq15[n,315]]]]
]

Plot[(rcpp[x]-(65536/x)*32768)/((65536/x)*32768),{x,65536,131071}]

Sqrt[Mean[Table[Abs[N[((rcpp[x]-(65536/x)*32768)/((65536/x)*32768))^2]],{x,65536,131071}]]]

Max[Table[Abs[N[(rcpp[x]-(65536/x)*32768)/((65536/x)*32768)]],{x,65536,131071}]]

Max[Table[Abs[N[rcpp[x]-(65536/x)*32768]],{x,65536,131071}]]

(*This was JM's original approximation, which was quite good. The major discrepancy is that it optimized for absolute, not relative error.*)

rcpp[x_]:=Module[{n},
n=x-65536-32768;
21848+mq15[n,-7251+mq15[n,2403+mq15[n,-934+mq15[n,327]]]]
]

Plot[(rcpp[x]-(65536/x)*32768)/((65536/x)*32768),{x,65536,131071}]

Sqrt[Mean[Table[Abs[N[((rcpp[x]-(65536/x)*32768)/((65536/x)*32768))^2]],{x,65536,131071}]]]

Max[Table[Abs[N[(rcpp[x]-(65536/x)*32768)/((65536/x)*32768)]],{x,65536,131071}]]

Max[Table[Abs[N[rcpp[x]-(65536/x)*32768]],{x,65536,131071}]]

(*Reciprocal approximation using linear interpolation followed by two rounds of Newton's method. This gets us to the level of truncation error with just 1 more multiply than the polynomial approximation.*)

f[x_]:=2/(x+1)

xi = N[(Table[-Cos[i*π/2], {i, 0,2}])+1]/2

pi={p0,p1}/.(Solve[({p0,p1}.{1,#}&/@xi)-(f/@xi)=={ϵ,-ϵ,ϵ}*(f/@xi),{p0,p1,ϵ}])[[1]]

{1.8823529411764706`,-0.9411764705882353`}

((fi/@xi)-(f/@xi))/(f/@xi)

xi=Join[{xi[[1]]},#-((D[(fi[x]-f[x])/f[x],x]/D[(fi[x]-f[x])/f[x],x,x]/.x->#)&/@#)&[Take[xi,{2,-2}]],{xi[[-1]]}]

Plot[(fi[x] - f[x])/f[x], {x, 0,1}]

Round[pi*16384]

rcp[x_]:=Module[{n,r},
n=x-32768;
r=30840+mq15[-15420,n];
r=r-mq15[r,mq15[n,r]+(r-32768)];
r=r-(mq15[r,mq15[n,r]+(r-32768)]+1)
]

Plot[(rcp[x]-(32768/x)*32768)/((32768/x)*32768),{x,32768,65535}]

Sqrt[Mean[Table[Abs[N[((rcp[x]-(32768/x)*32768)/((32768/x)*32768))^2]],{x,32768,65535}]]]

Max[Table[Abs[N[(rcp[x]-(32768/x)*32768)/((32768/x)*32768)]],{x,32768,65535}]]

Max[Table[Abs[N[rcp[x]-(32768/x)*32768]],{x,32768,65535}]]