APL(NARS), 155 chars

⎕FPC←108⋄F←{0⌈⌊(⎕FPC-{1>∣⍵:0⋄1+2⍟∣⍵}⍵)÷3.323v}

r←R q;L;n;d;m;e
e←÷10*F q⋄L←⍟q×d←r←1+n←0v
r+←m←(÷n×{⍎1⌽'v1z',⍕⍵}n+1)×d×←L÷n+←1⋄→2×⍳e<∣m
    
⍪{(F⍵)⍕R⍵}¨10*1..31v

//46+16+26+46+21

{⍎1⌽'v1z',⍕⍵} should be the Riemann Zeta function that use one Nars builtin for build constant type 'Zeta'. R q would be the Riemann R function of question of, input q, that stop computation after F q decimal digits.

It seems to me:

1. ⎕fpc it is the lenght in bit of the float point number of type "NumbeRv"
2. {1>∣⍵:0⋄1+2⍟∣⍵}q it should be the lenght in base 2 of the integer part of number q; so the lenght in base 2 digit of the fractional part of q would be 0⌈⎕FPC-{1>∣⍵:0⋄1+2⍟∣⍵}q, and the lenght of fractional part of q in base 10 would be {0⌈⌊(⎕FPC-{1>∣⍵:0⋄1+2⍟∣⍵}⍵)÷3.323v}q that is F q
3. F q it should be the lenght in decimal digit of fractional part of each float type "NumbeRv" q
4. are enought 108 bits float for doing this coputation, ⎕fpc←108, float 128 bits are over sized.

For test:

  ⍪{(F⍵)⍕R⍵}¨10*1..31v 
4.5645831410050902398657746955840 
25.661633266924182593226797940356 
168.35944628116734806491331098673 
1226.9312183434331085542162581721 
9587.431738841973414351612923908  
78527.39942912770485887029214096  
664667.4475647477679853466998875  
5761551.867320169562308864959734  
50847455.42772142751394887577256  
455050683.3068469244631532415820  
4118052494.631400441761046107709  
37607910542.22591023474569601743  
346065531065.8260271978929257301  
3204941731601.689034750500754116  
29844570495886.92737822228672779  
279238341360977.1872302539272988  
2623557157055978.003875460015662  
24739954284239494.40252165144480  
234057667300228940.2346566885611  
2220819602556027015.401217592243  
21127269485932299723.73386404400  
201467286689188773625.1590118748  
1925320391607837268776.080252873  
18435599767347541878146.80335902  
176846309399141934626965.8309690  
1699246750872419991992147.221858  
16352460426841662910939464.57821  
157589269275973235652219770.5691  
1520698109714271830281953370.160  
14692398897720432716641650390.61  
142115097348080886394439772958.2  
             

the calculation that takes 1'' here... it is possible change the number of bit of big float in 64 bits for to have

  ⎕FPC←64⋄⍪{(F⍵)⍕R⍵}¨10*1..3v
4.56458314100509024 
25.6616332669241826 
168.359446281167348 

C (gcc), 354 bytes

#include <quadmath.h>
#define D __float128
#define U unsigned
D Z(U p,D q){D t=-q,y,d,r,w,e,h;U n,m;for(r=n=0,w=1,e=w/powq(10,p);(w*=0.5q)>e;++n,r+=w*y)for(y=d=m=1,h=n;m<=n;--h,++m)y+=powq(m+1,t)*(d*=-h/m);return r/=1-powq(2,1+t);}D Rr(D q){D r,l,d,e,m,n;for(m=d=r=n=1,l=logq(q),e=r/powq(10,31);e<fabsq(m);++n,r+=m)m=(d*=l/n)*(1/(n*Z(31,n+1)));return r;}

Try it online!

This should be the port with some change of APL NARS answer... a little long. Here printf can not print the number so it was need PP() for print in the screen the number. I minimize the function Rr and Z, skipping all ceck for arguments, but it seems take more than 1 minute so no result showed from the link...

This below is the version ungolfed, more long (only for 128 bit floats too) that ceck arguments but faster than above, and use the function Zeta could be right even for complex not only for the real part. 10 second is the run, so it is showed the result.

C (gcc), 616 bytes

#include <quadmath.h>
#define U  unsigned
#define LD  __float128
#define CP __complex128
#define FM  FLT128_MAX
#define R   return
#define F   for
CP zeta(U p,CP q)
{CP t=-q,y,r,s[130],u;
 LD w,e,d,h;
 U  n,m,k;
 if(q==1||p>38)R FM;
 p+=p==0;k=2+p*3.322; // 38*3.322+2=128
 F(m=1,s[0]=0;m<k;++m)s[m]=cpowq(m,t);
 F(r=n=0,w=1,e=w/powq(10,p);(w*=0.5q)>e;++n,r+=w*y)
         F(y=d=m=1,h=n;m<=n;--h,++m)y+=(m+1>=k?cpowq(m+1,t):s[m+1])*(d*=-h/m);
 R r/=1-cpowq(2,1+t);
}
    
 LD Rr(LD q)
{LD r,l,d,e,m,t,n;
 F(m=d=r=n=1,l=logq(q),e=r/powq(10,31);e<fabsq(m);++n,r+=m) 
           m=(d*=l/n)*(1/(n*zeta(31,n+1)));
 R r;
}

Try it online!

C (gcc), 193 183 173 165 154 145 144 bytes

typedef _Float128 X;X logq();i,k;r(X x,X*y){for(X a=*y=k=1,q;i=999,k<i;*y-=(ldexp(a*=logq(x)/k++,-k)-a)/q)for(q=0;--i;)q+=(i%2?k:-k)*pow(i,~k);}

Try it online!

This estimates the Riemann R function from the Gram series. Note that this uses a horribly obfuscated but faster converging implementation of the Riemann zeta function.

$$\zeta(s)=\frac{1}{1-2^{1-s}}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n^s}$$

-1 thanks to @Simd. -8 thanks to @c--.

Slightly golfed less:

typedef _Float128 X;
X logq();
i,k;
r(X x,X*y){
  for(X a=*y=k=1,q;i=999,k<i;){
    /* estimate the Riemann zeta function */
    for(q=0;--i;)
      q+=(i%2?k:-k)*pow(i,~k);
    *y-=(ldexp(a*=logq(x)/k++,-k)-a)/q
  }
}

Julia, 141 138 133 125 81 bytes

h(x,a=big(1.))=1-sum(k->(a*=log(x)/k)/sum(n->(-1)^n/n/n^k,1:99)/k*(1-2^-k),a:200)

-44 thanks to @MarcMush

Try it online!

Slightly less golfed.

function h(x)
  B=BigFloat
  a=y=B(1);
  for K=1:200
    q=0;
    k=B(K);
    for n=1:99
      q-=(-1)^n/n/n^k
    end;
    a*=log(x)/k;
    y+=a/q/k*(1-2^-k)
  end;
  y
end

Julia using zeta() builtin, 209 205 bytes

using Base.MPFR
function h(x)B=BigFloat;a=y=B(1);for K=1:187 k=B(K);a*=log(x)/k;z=B();ccall((:mpfr_zeta,:libmpfr),Int8,(Ref{BigFloat},Ref{BigFloat},Int8),z,k+1,Base.MPFR.ROUNDING_MODE[]);y+=a/k/z end;y end

Try it online!

Slightly golfed less.

using Base.MPFR
function h(x)
  B=BigFloat;
  a=y=B(1);
  for K=1:187
    k=B(K);
    a*=log(x)/k;
    z=B();
    ccall((:mpfr_zeta,:libmpfr),Int8,(Ref{BigFloat},Ref{BigFloat},Int8),z,k+1,Base.MPFR.ROUNDING_MODE[]);
    y+=a/k/z
  end;
  y
end

Python + SymPy, 103 99 bytes

lambda x:1+Sum(ln(x)**k/k/gamma(k+1)/zeta(k+1),(k,1,oo)).evalf(50)
from sympy import*;k=Symbol('k')

Attempt This Online! Takes ~30s on ATO for all 31 test cases.

• Thanks to Simd for posting this chat message
• -4 thanks to ceilingcat

Wolfram Language (Mathematica), 8 bytes

RiemannR

Try it online!

Python + mpmath, 28 bytes

from mpmath import*
riemannr

Attempt This Online!

Not using the builtin, 84 bytes

from mpmath import*
f=lambda x:1+nsum(lambda k:log(x)**k/k/fac(k)/zeta(k+1),[1,inf])

Attempt This Online!

Vyxal, 44 bytes

Þ∞KƛƛÞ∞K$›eĖṠ2l≬₈vḞ≈c*n¡*?∆Lne$/;∑;2l≬₈vḞ≈c›

Try it Online!

I think this is correct. Might be 38 if there's an exact limit convergance somehow.

Precision is met by having a) things evaluated to 256 decimal places when approximating and b) exact values used until an approximation is needed. Good luck getting this to return an actual answer in the time we have left in the universe. The algorithm should be correct though.

Explained

The main idea to find the sums of things with an infinite upper bound is to check every overlapping pair of items in an infinite list of the sum applied from 1..1, 1..2, 1..3, and so on until the pair has all the same items. This is basically checking for convergence manually.

Þ∞Kƛ...;2l≬₈vḞ≈c›
Þ∞Kƛ   ;          # Calculate the gram series for all possible upper bounds
        2l        # get all the overlapping pairs
          ≬₈vḞ≈c  #  and get the first where the items to 256 decimal places are the same
                › # increment

ƛÞ∞K$e›ĖṠ2l≬₈vḞ≈c*n¡*?∆Lne$/;∑  # Note that the context variable is set to whatever number in a prefix in the prefix is being gram seriesed is being zeta'd
 Þ∞K                            # To each prefix of an infinite list of positive integers
    $e›ĖṠ2l≬₈vḞ≈c                # Zeta function
                 *              # times k
                  n¡*           # times k!
                     ?∆Lne$/    # log(input) ** k divided by above
                            ;∑  # sum the result of apply to each k

$e›ĖṠ2l≬₈vḞ≈c # the top of the stack is the prefix list
$e›           # each number in each prefix to the power of k + 1
   Ė          # reciprocal of each number in each prefix
    Ṡ         # sum of each prefix
     2l       # overlapping pairs of sums
       ≬₈vḞ≈c # first item where pair is all the same to 256 decimal places.

PARI/GP, 40 bytes

x->1+suminf(k=1,log(x)^k/k/k!/zeta(k+1))

Attempt This Online!

Using the Gram series.