Consider the following integral:


And consider the functions :



I want to compare them with each other (at least numerically for a large interval of values)

If graph for very large intervals (up to at least $10^4$) possible please add (please add three graphs in one axis system, so that I can compare)

(Does numerics suggest $S(q) \sim R(q)$ or $T(q)$? )

Related: Towards a new proof of infinitude of primes ( with possible unified application to other primes of special forms whose Infinitude is unknown) on Mathematics SE.

Note: Can't calculate the first integral on Mathematica for large values.

See this MSE post and this MSE post for more details.


3 Answers 3


There is no need for numerical computation here.

First, $T(q)$ is a well-known function, the logarithmic integral. Repeated integration by parts gives an asymptotic expansion

$$\mathrm{Li}(q) = \frac{q}{\log q}\sum_{k=0}^{K-1} \frac{k!}{\log^k q} + O\left(\frac{q}{\log^{K+1}q}\right).$$

There's also a fairly rapidly convergent representation due to Ramanujan which you can find on Wikipedia

Second, regarding the main integral, it has a different asymptotic. First, $\sin^2\theta = \frac{1-\cos(2\theta)}{2}$, we remove the $\frac{q-2}{2}$ coming from the constant and concentrate on getting cancellation in the oscillatory part. Second, recall the digamma function (the logarithmetic derivative of the gammafunction) $\digamma(x) = \frac{\Gamma'(x)}{\Gamma(x)}$ which satisfies $\log x - \frac{1}{x}\leq\digamma(x) \leq \log x - \frac{1}{2x}$ and $\digamma'(x) = \frac1x + \frac{1}{2x^2} +O(x^{-3})$

Letting $u = \frac{\Gamma(x)}x$ we have $\frac{du}{u} = d(\log u) = (\digamma(x)-\frac1x)dx$ so that

$$\begin{split}\frac{q-2}{2} - S(q) &= \int_{x=2}^{x=q} \cos\left(\pi\frac{\Gamma(x)}{2x}\right)dx \\ &= \int_{x=2}^{x=q} \frac{\cos(\pi u)}{u}\frac{du}{\digamma(x)-\frac1x} \end{split}$$

We now integrate by parts and get

$$\begin{split} &= \left[-\frac{\sin(\pi u)}{\pi u}\frac{1}{\digamma(x)-\frac1x}\right]_{x=2}^{x=q} - \int_{x=2}^{x=q} \frac{\sin(\pi u)}{\pi u^2}\frac{du}{\digamma(x)-\frac1x} \\ &- \int_{x=2}^{x=q} \frac{\sin(\pi u)}{\pi u}\frac{\digamma'(x)+\frac1{x^2}}{\left(\digamma(x)-\frac1x\right)^2} \frac{dx}{du}du \\ & = \left[-\frac{\sin(\pi u)}{\pi u}\frac{1}{\digamma(x)-\frac1x}\right]_{x=2}^{x=q} - \int_{x=2}^{x=q} \frac{\sin(\pi u)}{\pi u^2}\frac{du}{\digamma(x)-\frac1x} \\ &- \int_{x=2}^{x=q} \frac{\sin(\pi u)}{\pi u}\frac{\digamma'(x)+\frac1{x^2}}{\left(\digamma(x)-\frac1x\right)^3} du \end{split}$$

The first term is $O(1) + O\left(\frac{q}{\Gamma(q)}\right)$ and in particular is bounded. The second is similarly $O\left(\int_{x=2}^{x=q} \frac{du}{u^2}\right) = O(1)+ O\left(\frac{q}{\Gamma(q)}\right)$. For the last term divide the interval into two parts: up to $2\leq x\leq q^\delta$ and $q^\delta \leq x \leq q$ for some $\delta < 1$. On the first interval we use that $\digamma'(x)+\frac1{x^2} = O(\frac{1}{x}) = O(1)$ to bound the integral as $O(1)+O(\log u(q^\delta) = O(\log(\Gamma(q^\delta)) = O(\delta q^\delta \log q)$. On the second interval we have $\digamma'(x)+\frac1{x^2} = O(q^{-\delta})$ so the whole integral is $O(q^{1-\delta}\log q)$. Taking $\delta = \frac12$ we conclude that

$$ S(q) = \frac{q-2}{2} + O(q^{1/2}\log q)$$

And in particular has a different asymptotic.

Finally, a more careful analysis using the $\log^3x $ in the denominator of the third term would give the error term $O\left(\frac{q^{1/2}}{\log^2 q}\right)$.

  • $\begingroup$ I think multiple integrations by parts will give a $S(q) = \frac12 q + O_\epsilon(q^\epsilon)$ for any $\epsilon>0$. $\endgroup$ May 23, 2020 at 21:34
  • $\begingroup$ this is a great answer , thank you . I also added the links behind this question . Please comment your valuable advice . $\endgroup$
    – bambi
    May 28, 2020 at 13:58
  • 1
    $\begingroup$ The problem is that the connection to prints only occurs at integral values of $x$, and the function is so oscillatory that the sum over the integers is a terrible approximation to the integral. $\endgroup$ May 28, 2020 at 20:25
  • $\begingroup$ thank you for the reply . But you can see the post on growth condition on the functional, which seems optimistic. $\endgroup$
    – bambi
    May 29, 2020 at 20:56

You write $S(q)$ and $T(q)$ as integrals, but it is easier to think of them as solutions of ODEs: $$ S'(q) = \sin^2\left(\frac{π\Gamma(q)}{2q}\right) $$ with initial conditions $$ S(2) = 0, $$ and similarly for $T(q)$. You can then use any of the common ODE integrators in matlab, mathematica, maple, ..., to solve and plot the solutions so that you can compare them.

  • $\begingroup$ This is a great idea. Unfortunately, I don't think it solves the problem. See my edit above. $\endgroup$
    – user20857
    May 15, 2020 at 22:19
  • $\begingroup$ @SpencerBryngelson Why not? Maybe I don't see what your edit is, but I see no reason why you couldn't use this approach to compare the three functions. Can you explain what precisely doesn't work for you? $\endgroup$ May 16, 2020 at 3:45
  • 1
    $\begingroup$ I think the RHS has something like 10^5 oscillations by $q=50$. If you hope to resolve those, you’d need at least 10x more time steps than that. For larger $q$, the number of oscillations grows more rapidly still. I just don’t think it’s possible for the $q$ values OP mentioned. $\endgroup$
    – user20857
    May 16, 2020 at 7:02
  • 1
    $\begingroup$ @SpencerBryngelson -- well, but then no other integration method will be able to do it either unless you use particular properties of the integrand. For example, you could see whether you can approximating the integrand over one period of the sine and think about whether that helps you to approximate the integral by counting how many periods you cover from $q=2$ to wherever you need to go. $\endgroup$ May 16, 2020 at 15:34
  • $\begingroup$ agreed. That was the spirit of my comment about oscillatory integration methods, which do just that. $\endgroup$
    – user20857
    May 17, 2020 at 0:21

The only chance you stand to deal with this problem from a numerical perspective is oscillatory integration methods. Filon/Levin-type methods can sometimes handle problems like this, particularly when they are of $\sin$ or $\cos$ type, though the $\Gamma$ function and its run-away growth may be prohibitive for the $q$ you are hoping for. In any case, I was able to accurately evaluate the integral with Mathematica using Levin's method up to around $q = 50$ via

NIntegrate[Sin[Pi Gamma[x]/(2 x)]^2, {x, 2, 50}, MinRecursion -> 9, Method -> "LevinRule"]

  • If you want to analyze this from an analytical perspective, Math.SE is likely more appropriate
  • Mathematica can likely get you further if you can derive the amplitude/oscillation matrices and specify them in the function call. See, e.g. here. I'm not sure if this is possible in your case or not, but reading the classic papers by Levin from 1996/7 might be helpful.
  • Further questions regarding Mathematica implementation should probably be directed to Mathematica.SE


Wolfgang had the clever idea to just solve this in ODE form. However, I don't think solves the problem in the end. The integrand (or RHS in derivative form) is just too highly oscillatory for even modest $q$. For example, using standard methods:

NDSolve[{s'[x] == Sin[Pi Gamma[x]/(2 x)]^2, s[2] == 0}, s[x], {x, 2, q}, MaxSteps -> 10^7]

I am only able to reach $q=17$ in these $10^7$ time steps. The problem gets much worse for still larger $q$. A quick plot of number of oscillations per $q$ reveals this. I think somehow re-casting the problem in a non-oscillatory form (e.g. Levin's method) in the only real way forward, though I'm happy to learn otherwise.

  • $\begingroup$ I already posted the question on MSE but no answers $\endgroup$
    – bambi
    May 15, 2020 at 7:53

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.