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.


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$ – Lior Silberman May 23 '20 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 '20 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$ – Lior Silberman May 28 '20 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 '20 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$ – Spencer Bryngelson May 15 '20 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$ – Wolfgang Bangerth May 16 '20 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$ – Spencer Bryngelson May 16 '20 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$ – Wolfgang Bangerth May 16 '20 at 15:34
  • $\begingroup$ agreed. That was the spirit of my comment about oscillatory integration methods, which do just that. $\endgroup$ – Spencer Bryngelson May 17 '20 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 '20 at 7:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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