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.


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.

Related challenges


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.

  • 3
    \$\begingroup\$ at least one decimal place looks too high, requiring custom type(higher than long double) \$\endgroup\$
    – l4m2
    Commented Dec 28, 2022 at 12:13
  • 3
    \$\begingroup\$ It's still out of long double. Usually a relative error bound is provided for such kind \$\endgroup\$
    – l4m2
    Commented Dec 28, 2022 at 12:20
  • 4
    \$\begingroup\$ 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 \$\endgroup\$
    – jezza_99
    Commented Dec 28, 2022 at 20:57
  • 5
    \$\begingroup\$ 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\$ \$\endgroup\$
    – jezza_99
    Commented Dec 28, 2022 at 21:04
  • 6
    \$\begingroup\$ 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. \$\endgroup\$
    – Bubbler
    Commented Dec 29, 2022 at 6:19

9 Answers 9


PARI/GP, 40 bytes


Attempt This Online!

Using the Gram series.

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

Wolfram Language (Mathematica), 8 bytes


Try it online!

  • 6
    \$\begingroup\$ @Joao-3 That request was edited in after Adám's answer was submitted. \$\endgroup\$
    – chunes
    Commented Dec 28, 2022 at 14:57
  • 1
    \$\begingroup\$ @chunes I understand, but LOOK! \$\endgroup\$
    – Joao-3
    Commented Dec 28, 2022 at 19:21
  • 1
    \$\begingroup\$ The precision of this solution isn't high enough at 10^31. Something like Round@*RiemannR would satisfy the precision requirement. \$\endgroup\$
    – Roman
    Commented Jan 5, 2023 at 9:39
  • 2
    \$\begingroup\$ 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[10`99^31], to get the correct answer. I'm not sure that's allowed by the rules. \$\endgroup\$
    – Roman
    Commented Jan 5, 2023 at 20:53
  • 1
    \$\begingroup\$ @Roman I don't see why requiring that input form would be against default rules. \$\endgroup\$
    – Adám
    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.


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

Slightly golfed less:

typedef _Float128 X;
X logq();
r(X x,X*y){
  for(X a=*y=k=1,q;i=999,k<i;){
    /* estimate the Riemann zeta function */
  • \$\begingroup\$ Wow. This is amazing. \$\endgroup\$
    – Simd
    Commented Jan 2, 2023 at 5:57
  • \$\begingroup\$ How about _Float128 instead of __float128? \$\endgroup\$
    – Simd
    Commented Jan 2, 2023 at 7:16
  • \$\begingroup\$ It doesn't help with the golfing but it seems 99 can be reduced to 23. \$\endgroup\$
    – Simd
    Commented Jan 2, 2023 at 8:36
  • \$\begingroup\$ The answer for 10^1 is not accurate to the required precision. \$\endgroup\$
    – Simd
    Commented Jan 6, 2023 at 11:04
  • \$\begingroup\$ the algo for 999 cicles of the loop not go out than 4 digits after the "," example 16 279238341360977.187230253927299 10^16 279238341360977.187244069148541155 \$\endgroup\$
    – Rosario
    Commented Feb 25 at 11:01

Vyxal, 44 bytes


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.


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ƛ   ;          # 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.
  • \$\begingroup\$ The link seems to time out . \$\endgroup\$
    – Simd
    Commented Jan 1, 2023 at 13:46
  • \$\begingroup\$ @simd it never will not time out. \$\endgroup\$
    – lyxal
    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.

  • \$\begingroup\$ What do you get if you implement zeta too? \$\endgroup\$
    – Simd
    Commented Jan 1, 2023 at 19:00
  • \$\begingroup\$ @Simd sorry, I don't know much about SymPy. What do mean by that? \$\endgroup\$
    – The Thonnu
    Commented Jan 1, 2023 at 19:01
  • \$\begingroup\$ You are using the function zeta from sympy. I was just wondering how much longer the code would be if you implemented zeta yourself. \$\endgroup\$
    – Simd
    Commented Jan 1, 2023 at 19:06
  • \$\begingroup\$ @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 \$\endgroup\$
    – The Thonnu
    Commented Jan 1, 2023 at 19:31
  • \$\begingroup\$ That looks cool. Could you add that as an alternative answer? \$\endgroup\$
    – Simd
    Commented Jan 1, 2023 at 19:39

Julia, 141 138 133 125 81 bytes


-44 thanks to @MarcMush

Try it online!

Slightly less golfed.

function h(x)
  for K=1:200
    for n=1:99

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)
  for K=1:187
  • \$\begingroup\$ Now Rust please …. \$\endgroup\$
    – Simd
    Commented Jan 5, 2023 at 9:24
  • 1
    \$\begingroup\$ @Simd I haven’t quite worked out how to use f128 in TIO or ATO. \$\endgroup\$
    – ceilingcat
    Commented Jan 5, 2023 at 15:43
  • \$\begingroup\$ This is not right for 10^1 . \$\endgroup\$
    – Simd
    Commented Jan 6, 2023 at 11:07
  • \$\begingroup\$ The new version works! \$\endgroup\$
    – Simd
    Commented Jan 9, 2023 at 7:55
  • \$\begingroup\$ 81 bytes, mainly by replacing the for loops with sums Try it online! \$\endgroup\$
    – MarcMush
    Commented Feb 13, 2023 at 23:45

Python + mpmath, 28 bytes

from mpmath import*

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!

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

APL(NARS), 330 chars


r←p z q;n;t;w;m;d;y;h;e

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


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)



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


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)


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
 R r/=1-cpowq(2,1+t);
 LD Rr(LD q)
{LD r,l,d,e,m,t,n;
 R r;

Try it online!

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.