3
$\begingroup$

From "The First 50 Million Prime Numbers" by Don Zagier: primes are integral roots of$$ 1-\frac{\sin(\frac{\pi\Gamma(s)}s)}{\sin(\frac\pi s)}. $$ The graph of this function looks like

enter image description here

I would like to produce the sequence of all positive real roots of this function in a given interval (say, from 10 to 12), with some given precision (say, $10^{-10}$). What is an efficient way to do this?

$\endgroup$
3
  • 2
    $\begingroup$ What about NSolve[{1 == Sin[Pi Gamma[s]/s]/Sin[Pi/s], 10 < s < 12}, s, Reals,WorkingPrecision->10]? $\endgroup$ Oct 24 '19 at 21:06
  • $\begingroup$ @UlrichNeumann Looks impressive, but 11 must be among the roots, and with me it goes like ..., 10.99682955, 10.99852262, 11.00128148, 11.00147394, ... $\endgroup$ Oct 24 '19 at 21:14
  • 1
    $\begingroup$ Since you are only interested in integral roots: Solve[{1 == Sin[Pi Gamma[s]/s]/Sin[Pi/s], 10 < s < 12}, s, Integers] $\endgroup$
    – Bob Hanlon
    Oct 24 '19 at 22:16
7
$\begingroup$
Clear["Global`*"]

f[s_] = 1 - Sin[Pi Gamma[s]/s]/Sin[Pi/s];

Since you are interested in integral roots

sol = s /. Solve[{f[s] == 0, 2 <= s < 100}, s, Integers]

{* {2, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 
    61, 67, 71, 73, 79, 83, 89, 97} *)

The integral roots are prime

And @@ PrimeQ[sol]

(* True *)

These are all of the primes in the interval

Rest@sol == Prime@Range[Length@sol - 1]

(* True *)

Rest@sol == NestList[NextPrime, 2, Length@sol - 2]

(* True *)

EDIT: To find the real roots use high precision for calculations, then reduce the precision for display.

solr = NSolve[{f[s] == 0, 10 < s < 12}, s, WorkingPrecision -> 100,
   VerifySolutions -> True];

Verifying the solutions

And @@ (f[s] == 0 /. solr)

(* True *)

The roots are dense

Length@solr

(* 1490 *)

Nonetheless, NSolve misses the integral root.

solr[[740 ;; 750]] // N

(* {{s -> 10.9867}, {s -> 10.9912}, {s -> 10.9924}, {s -> 10.9929}, {s -> 
   10.9968}, {s -> 10.9985}, {s -> 11.0013}, {s -> 11.0015}, {s -> 
   11.0024}, {s -> 11.0032}, {s -> 11.0033}} *)

f[11]

(* 0 *)

Combining the real and integral solutions

sol = Join[solr,
    Solve[{f[s] == 0, 10 < s < 12}, s, Integers]] //
   SortBy[#, Last] &;

sol[[740 ;; 750]] /. x_Real :> N[x]

(* {{s -> 10.9867}, {s -> 10.9912}, {s -> 10.9924}, {s -> 10.9929}, {s -> 
   10.9968}, {s -> 10.9985}, {s -> 11}, {s -> 11.0013}, {s -> 11.0015}, {s -> 
   11.0024}, {s -> 11.0032}} *)

Alternatively, do a search with FindRoot

solf = Union[
   FindRoot[f[s] == 0, {s, #},
      WorkingPrecision -> 100] & /@
    Range[10, 12, 10^-4],
   SameTest ->
    (Abs[#1[[1, -1]] - #2[[1, -1]]] < 10^-4 &)];

This is much slower but identifies many more roots

Length@solf

(* 12440 *)

including the integral root

solf[[6200 ;; 6205]] // N

(* {{s -> 10.9996}, {s -> 10.9998}, {s -> 11.}, {s -> 11.0002}, {s -> 
   11.0003}, {s -> 11.0005}} *)

solf[[6202]]

(* {s -> 11.000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000} *)
$\endgroup$
0
1
$\begingroup$

Because of the highly oscillatory nature of the function under consideration, an NDSolve[]-based method with event location (discussed by Daniel in this answer) seems appropriate, as it is more likely to be careful about traversing the function under consideration. Of course, the price to pay for this carefulness is that the evaluations take quite a while.

zagier[u_] := 1 - Sin[π Gamma[u]/u]/Sin[π/u]

With[{u0 = 10, u1 = 12}, 
     AbsoluteTiming[rts = Reap[NDSolve[{y'[u] == zagier'[u], y[u0] == zagier[u0], 
                                        WhenEvent[y[u] == 0, Sow[u],
                                                  "LocationMethod" -> "Brent"]},
                                       {}, {u, u0, u1}, MaxSteps -> ∞,
                                       WorkingPrecision -> 20]][[-1, 1]];]]
   {36267.7, Null}

Length[rts]
   3290112

Norm[zagier[Take[rts, 10]], ∞]
   1.95*10^-9

Norm[zagier[Take[rts, -10]], ∞]
   0.0013240266

We see that the accuracy has degraded for the larger roots; one can still use these results as seeds for FindRoot[] to polish, however.

Nearest[rts, 11]
   {10.999999999995645962}
$\endgroup$

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.