# Implement the Riemann R function

The Riemann R function is as follows:

$$R (x)=\sum _{n=1}^{\infty } \frac{\mu (n) \text{li}\left(x^{1/n}\right)}{n}.$$

This uses the Möbius function as well as the logarithmic integral.

From Wikipedia, the Möbius function is defined so that for any positive integer $$\n\$$, $$\μ(n)\$$ is the sum of the primitive nth roots of unity. It has values in $$\\{−1, 0, 1\}\$$ depending on the factorization of $$\n\$$ into prime factors:

• $$\μ(n) = +1\$$ if $$\n\$$ is a square-free positive integer with an even number of prime factors.
• $$\μ(n) = −1\$$ if $$\n\$$ is a square-free positive integer with an odd number of prime factors.
• $$\μ(n) = 0\$$ if $$\n\$$ has a squared prime factor.

The logarithmic integral is defined as:

$$\operatorname{li}(x) = \int_2^x \frac{dt}{\log t}.$$

An alternative way to compute the Riemann R function is via the Gram series. That is:

$$R(x) = 1 + \sum_{k=1}^{\infty}\frac{(\ln x)^k}{k k!\zeta(k+1)}.$$

The function $$\\zeta()\$$ is the Riemann zeta function.

# Challenge

Write code that computes $$\R(x)\$$ for $$\x\$$ up to $$\10^{31}\$$. The output may have a margin of error of up to 1.

# Test cases

Here are the answers for $$\10^i\$$ for $$\i = \{1, 2, \dots, 31\}\$$.

1 4.56458314100509023986577469558
2 25.6616332669241825932267979404
3 168.359446281167348064913310987
4 1226.93121834343310855421625817
5 9587.43173884197341435161292391
6 78527.399429127704858870292141
7 664667.447564747767985346699887
8 5761551.86732016956230886495973
9 50847455.4277214275139488757726
10 455050683.306846924463153241582
11 4118052494.63140044176104610771
12 37607910542.2259102347456960174
13 346065531065.82602719789292573
14 3204941731601.68903475050075412
15 29844570495886.9273782222867278
16 279238341360977.187230253927299
17 2623557157055978.00387546001566
18 24739954284239494.4025216514448
19 234057667300228940.234656688561
20 2220819602556027015.40121759224
21 21127269485932299723.733864044
22 201467286689188773625.159011875
23 1925320391607837268776.08025287
24 18435599767347541878146.803359
25 176846309399141934626965.830969
26 1699246750872419991992147.22186
27 16352460426841662910939464.5782
28 157589269275973235652219770.569
29 1520698109714271830281953370.16
30 14692398897720432716641650390.6
31 142115097348080886394439772958.0


Your code doesn't need to be fast, but ideally it should complete in under a minute.

# Precision

You will need 128 bit floats to represent the output. In C __float128 from quadmath.h is the simplest way (long double will most likely be 80 bits). Other languages may have standard libraries to support 128 bit floats (e.g Decimal in Python). challenges are judged per language so there is no penalty in using whatever is needed for your favorite langauge.

• at least one decimal place looks too high, requiring custom type(higher than long double)
– l4m2
Commented Dec 28, 2022 at 12:13
• It's still out of long double. Usually a relative error bound is provided for such kind
– l4m2
Commented Dec 28, 2022 at 12:20
• Saying the solution shouldn't be fast but needs to not timeout on TIO is a bad rule in my opinion. If someone comes up with a valid, well-golfed solution that is slow, it shouldn't be invalid because it times out on TIO Commented Dec 28, 2022 at 20:57
• I agree with @l4m2, I think a relative error is much better than a blanket +-1. Otherwise the required accuracy of the solution increases with $i$ Commented Dec 28, 2022 at 21:04
• This challenge made me write this thing to avoid when writing challenges: real-valued output with bad specification. Another thing to avoid: requiring time limits. The site culture is that theoretically correct solutions are fine, even if you can't see the results, as long as explanation is provided. Commented Dec 29, 2022 at 6:19

# PARI/GP, 40 bytes

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


Attempt This Online!

Using the Gram series.

• An excellent answer benefitting from pari’s strengths.
– Simd
Commented Dec 29, 2022 at 6:26

# Wolfram Language (Mathematica), 8 bytes

RiemannR


Try it online!

• @Joao-3 That request was edited in after Adám's answer was submitted. Commented Dec 28, 2022 at 14:57
• @chunes I understand, but LOOK! Commented Dec 28, 2022 at 19:21
• The precision of this solution isn't high enough at 10^31. Something like Round@*RiemannR would satisfy the precision requirement. Commented Jan 5, 2023 at 9:39
• No, it's not a printing issue. RiemannR[10.^31] returns the machine-precision number 1.421150973480811*^29, which is far off the mark (the error is about 2×10^14). You'd have to call it with an arbitrary-precision argument, like RiemannR[1099^31], to get the correct answer. I'm not sure that's allowed by the rules. Commented Jan 5, 2023 at 20:53
• @Roman I don't see why requiring that input form would be against default rules.
Commented Jan 6, 2023 at 7:40

# 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
}
}

• Wow. This is amazing.
– Simd
Commented Jan 2, 2023 at 5:57
– Simd
Commented Jan 2, 2023 at 7:16
• It doesn't help with the golfing but it seems 99 can be reduced to 23.
– Simd
Commented Jan 2, 2023 at 8:36
• The answer for 10^1 is not accurate to the required precision.
– Simd
Commented Jan 6, 2023 at 11:04
• the algo for 999 cicles of the loop not go out than 4 digits after the "," example 16 279238341360977.187230253927299 10^16 279238341360977.187244069148541155 Commented Feb 25 at 11:01

# 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.

• The link seems to time out .
– Simd
Commented Jan 1, 2023 at 13:46
• @simd it never will not time out. Commented Jan 1, 2023 at 22:26

# 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.

• What do you get if you implement zeta too?
– Simd
Commented Jan 1, 2023 at 19:00
• @Simd sorry, I don't know much about SymPy. What do mean by that? Commented Jan 1, 2023 at 19:01
• You are using the function zeta from sympy. I was just wondering how much longer the code would be if you implemented zeta yourself.
– Simd
Commented Jan 1, 2023 at 19:06
• @Simd oh. I think we could have f=lambda a:sum(m**-a for m in range(1,10**6)) (45 bytes) - inspired by this answer Commented Jan 1, 2023 at 19:31
• That looks cool. Could you add that as an alternative answer?
– Simd
Commented Jan 1, 2023 at 19:39

# 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

– Simd
Commented Jan 5, 2023 at 9:24
• @Simd I haven’t quite worked out how to use f128 in TIO or ATO. Commented Jan 5, 2023 at 15:43
• This is not right for 10^1 .
– Simd
Commented Jan 6, 2023 at 11:07
• The new version works!
– Simd
Commented Jan 9, 2023 at 7:55
• 81 bytes, mainly by replacing the for loops with sums Try it online! Commented Feb 13, 2023 at 23:45

# 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!

• Also works with mp.dps=30
– Simd
Commented Jan 1, 2023 at 17:44
• What would you get it if you implemented zeta too?
– Simd
Commented Jan 1, 2023 at 17:44

# APL(NARS), 330 chars

F←{0⌈(⌊1+10⍟⍵)-⍨⌊3.322÷⍨⎕fpc-1}

r←p z q;n;t;w;m;d;y;h;e
→2×⍳∼q=1⋄r←∞⋄→0
t←q×¯1v⋄r←q-q⋄e←÷10x*p⋄n←0v⋄w←1v
→7×⍳e>w×←0.5v
y←d←1v⋄h←n+1⋄m←0v
h-←1⋄y+←(t*⍨m+1)×d×←¯1×h÷m+←1⋄→5×⍳m<n
r+←w×y⋄n+←1⋄→3
r÷←1-2*1+t

r←R q;L;n;d;m;e;y
y←q×1v⋄r←1v⋄L←⍟y⋄n←0⋄d←1v⋄e←÷10x*2+F q
r+←m←(÷n×(⌊⎕fpc÷3.5)z n+1)×d×←L÷n+←1⋄→2×⍳e<∣m

⍪{{⍵⍕⍨1⌈¯8+F⍵}R⍵}¨10v*1..31


32+ 24+16+33+14+18+38+15+10 +18+39+45 +28=330

z should be the Riemann Zeta function that stop calc until p digits after the decimal point, with input one number that should be integer, float, big float,complex, big float clomplex and one output that has to be big float or big float complex (I'm not sure of this speech 100%). It seems if the number is not a complex the calculation in z function is done not use complex operations but only float or big float operations. I'm not sure is safe to use ∞ as result for zeta in case input 1.

R it seems to be the Riemann R function cut on the possible digits it can reach from the instruction (⌊⎕fpc÷4.1) for speed reason, seen what the question ask too.

It seems:

1. ⎕fpc it is the lenght in bit of the float point number of type "NumbeRv"

2. ⌊1+10⍟⍵ it should be the lenght in decimal digit of the integer part of number in ⍵

3. ⌊3.322÷⍨⎕fpc-1 it should be the lenght in decimal digit of each float type "NumbeRv"

4. F ⍵ it should be the max lenght in decimal digit of fractional part number in ⍵, so it is the right place where end the computation

5. It seems is asked one decimal ¯8 + teoretical possible 128 bit precision for the calculus and for show the number

6. In function (r←p z q), if it was used the instruction r←q-q the calculation in z not use big complex float, but only big float

7. It seems 2+F q is safe as precision for R q (not for all other numbers return z function in R)

   ⍪{{⍵⍕⍨1⌈¯8+F⍵}R⍵}¨10v*1..31


4.56458314100509023986577469558
25.6616332669241825932267979404
168.359446281167348064913310987
1226.93121834343310855421625817
9587.43173884197341435161292391
78527.3994291277048588702921410
664667.447564747767985346699888
5761551.86732016956230886495974
50847455.4277214275139488757726
455050683.306846924463153241582
4118052494.63140044176104610771
37607910542.2259102347456960174
346065531065.826027197892925730
3204941731601.68903475050075412
29844570495886.9273782222867278
279238341360977.187230253927299
2623557157055978.00387546001566
24739954284239494.4025216514448
234057667300228940.234656688561
2220819602556027015.40121759224
21127269485932299723.7338640440
201467286689188773625.159011875
1925320391607837268776.08025287
18435599767347541878146.8033590
176846309399141934626965.830969
1699246750872419991992147.22186
16352460426841662910939464.5782
157589269275973235652219770.569
1520698109714271830281953370.16
14692398897720432716641650390.6
142115097348080886394439772958.2

one calculation that takes 15' here... it is possible change the number of bit of big float for to have

   ⎕fpc←64
⍪{{⍵⍕⍨F⍵}R⍵}¨10v*⍳3
4.56458314100509191
25.6616332669242044
168.359446281167561


with the last 3 or 4 decimal digits wrong it seems for the instruction
(⌊⎕fpc÷4.1) that should be (⌊⎕fpc÷3.322) for to have what it seems to me the exact solution (with about 2x the time for calculus here some second only for this)

   ⍪{{⍵⍕⍨F⍵}R⍵}¨10v*⍳3
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!

• Why is this so much less golfed than the other C answer?
– Simd
Commented Feb 23 at 11:44
• @ceilingcat ok very well Commented Feb 27 at 21:00
• 244 bytes Commented Mar 3 at 5:14