Plotting Riemann-Siegel Formula with J

The Riemann-Siegel formula, Z(t) has been used to find non-trivial zeros of Riemann's zeta function.
wikipedia

The following source code gives a numerical calculation of \[Z(t) \approx 2 \sum_{n=1}^{N}n^{-\frac{1}{2}}\cos [\theta(t) - t\log n] + R,\]

where \[\theta(t) = \Im \log \Gamma (1 + \frac{it}{2}- \frac{3}{4}) - \frac{t}{2}\log\pi . \]  This can be simplified to an easily calculatable approximation given by:
\[\theta(t) \approx \frac{t}{2}\log (\frac{t}{2\ pi}) - \frac{t}{2} - \frac{\pi}{8} + \frac{1}{48t }+ \frac{7}{5760t^{3}}.\]

The final term, [\ R,\] in the original sum, is the remainder term, which is complicated to derive but  can be reduced to a sum,[
\[ R \approx (-1)^{N-1}(\frac{t}{2 \pi})^{-1/4}[C_{0} + C_{1}(\frac{t}{2 \pi})^{-1/2} +C_{2}(\frac{t}{2 \pi})^{-1} +

C_{3}(\frac{t}{2 \pi})^{-3/2} + C_{4}(\frac{t}{2 \pi})^{-2}].\]

We will use this approximation in the calculation.

NB. Calculation of Z(x), estimate of zeta(x) for x = a +ib where b = 0.5  
NB. Riemann-Siegel formula. (from Riemann Zeta Function, H.M Edwards)      

load 'plot'  
load 'numeric trig'  
    
C0 =: 0.38268343236508977173, 0.43724046807752044936, 0.13237657548034352333, _0.01360502604767418865, _0.01356762197010358088, _0.00162372532314446528, 0.00029705353733379691,0.00007943300879521469, 0.00000046556124614504, _0.00000143272516309551, _0.00000010354847112314, 0.00000001235792708384, 0.00000000178810838577, _0.00000000003391414393, _0.00000000001632663392  

C1 =: 0.02682510262837535, _0.01378477342635185, _0.03849125048223508, _0.00987106629906208, 0.00331075976085840, 0.00146478085779542, 0.00001320794062488,_0.00005922748701847, _0.00000598024258537, 0.00000096413224562, 0.00000018334733722
  
C2 =: 0.005188542830293, 0.000309465838807, _0.011335941078299, 0.002233045741958,0.005196637408862, 0.000343991440762, _0.000591064842747, _0.000102299725479, 0.000020888392217,0.000005927665493, _0.000000164238384, _0.000000151611998  

C3 =: 0.0013397160907, _0.0037442151364, 0.0013303178920, 0.0022654660765,_0.0009548499998, _0.0006010038459, 0.0001012885828, 0.0000686573345, _0.0000005985366,_0.0000033316599, _0.0000002191929, 0.0000000789089, 0.0000000094147 
 
C4 =: 0.00046483389, _0.00100566074, 0.00024044856, 0.00102830861,_0.00076578609, _0.00020365286, 0.00023212290, 0.00003260215, _0.00002557905, _0.00000410746,0.00000117812, 0.00000024456  
    
pi2 =: +: o. 1  
    
theta =: 3 : 0  
assert. y > 0  
p1 =. ( * -:@:^.@:(%&pi2)) y  
p2 =. --: y  
p3 =. -(o. 1) % 8  
p4 =. 1 % (48 * y)  
p5 =. 7 % (5760 * y ^ 3)  
p1 + p2 + p3 + p4 + p5  
)  
    
summation =: 3 : 0 "1  
   
v =. >0{y  
m =. >1{y  
thetav =. theta v  
ms =. >: i. m  
cos =. 2&o.  
s =. +:@:(^&_0.5) * cos@:(thetav&-@:(v&*@:^.))  
+/ s ms  
)  
    
remainderterms =: 3 : 0  
im =. y  
v =. %: pi2 %~ y  
N =. <. v  
p =. v - N  
q =. 1 - 2 * p  
final =. 0  
   
mv0 =. +/ (({&C0) * (q&^@:+:)) i. # C0  
final =. mv0  
  
mv1 =. +/ (({&C1) * (q&^@:+:@:>:)) i. # C1  
mv1 =. mv1 * ((im % pi2) ^ _0.5)  
final =. final + mv1  
   
mv2 =. +/ (({&C2) * (q&^@:+:)) i. # C2  
mv2 =. mv2 * ((im % pi2) ^ _1)  
final =. final + mv2  
    
mv3 =. +/ (({&C3) * (q&^@:+:@:>:)) i. # C3  
final =. final + mv3  
    
mv4 =. +/ (({&C4) * (q&^@:+:)) i. # C4  
final =. final + mv4  
   
    
final * (_1 ^ <: N) * ((im % pi2) ^ _0.25)  
)  
   
    
z =: 3 : 0 "0  
im =. y  
v =. %: pi2 %~ y  
N =. <. v  
p =. v - N  
v1 =. summation im ; N  
v2 =. remainderterms im  
v1 + v2  
)     
   
x =: steps 0.1 100 3500  
plot x; z x  

Graphical plot

The lines:

 x =: steps 0.1 100 3500  
 plot x; z x  

create a plot of the function from 0.1 to 100 in 3500 steps, as shown below:

You can see the zeros of the Zeta Function from this graph. The first zero is at about x = 14.29.

Plotting Zeta function on the complex plane

 zeta3d =: 3 : 0  
 if. (0{+.y) > 0.5 do.  
 I =. >: i. 100  
 +/I^-y  
 elseif. (0{+.y) < 0.5 do.  
 z =. zeta3d 1-y  
 g =. ! 1-y  
 s =. 1 o. -:(o.1)*y  
 c =. ((o.1)^(y-1))*2^y  
 r=.z*g*s*c  
 if. 40 < |r do.  
 r =. 40  
 end.  
 r  
 elseif. 1 do.  
 0  
 end.  
 )  
   
 absz =: |@:zeta3d  
 realz=: 0&{@:+.@:zeta3d  
 imz=: 1&{@:+.@:zeta3d  
 D =.steps _15 15 40  
 'surface; viewpoint 1 _1 2' plot imz"0 j./~ D NB. use realz for real graph.

Here is a plot of the imaginary parts:

....and the real parts.

Constructing Steiner Triple Systems with OCaml

A steiner system, S(a,b,c), where a,b,c are all positive integers with a <= b <= c, is defined as a set of c elements, a collection of subsets of this set, each of size b, such that every a elements exist in exactly one subset.

e.g. S(2,3,7)

We need a set of 7 elements,

S = {0,1,2,3,4,5,6},

and we need to create subsets of S, each of size 3, where each pair of elements exist in exactly one subset:

b1 = {0,1,2}

b2 = {0,3,4}

b3 = {0,5,6}

b4 = {1,4,5}

b5 = {1,3,6}

b6 = {2,4,6}

b7 = {2,3,5}

From observation, it is clear that any two elements exist in exactly one subset b1,...,b7. We shall refer to each subset as a block.

Clearly, a steiner system does not necessarily exist for any given a,b,c. For example, there is no S(2,7,8) system. You cannot generate subsets of {0,1,2,3,4,5,6,7}, each of size 7, such that any two elements are in exactly one subset.

An S(a,b,c) is called a steiner triple system if a = 2, b = 3, and can be denoted STS(n). It can be easily proven that a S(2,3,n) exists for any n such that

\[n\cong1,3\mod6\]
But, this fact doesn't actually give us the triples of the STS(n), which is what we are interested in. The above example's 7 blocks were found by trial and error, and this isn't really a scalable method for constructing triple systems. It is not immediately clear how we can construct an STS(n) for a given n. This post creates the algorithm(s) necessary to create a STS(n) for an any applicable n.

The algorithm
We first split the problem into two parts. We will solve the problem for the case of \(n\cong 1\mod 6\) and \(n\cong3\mod6\) separately.

Case 1: \[ n\cong3\mod6\] . This is sometimes known as the Bose Construction. I won't go into details but will just present the OCaml source code, as it is not particularly long or complicated:

 let construct_for_residue_3 v =  
     let res = (v - 3) / 6 in  
     let arr = ref [] in  
     let sz = (2*res) + 1 in  
     let mult = (sz + 1) / 2 in  
     for i = 0 to sz - 1 do  
         arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
     done;  
     for j = 0 to sz - 2 do  
         for k = (j+1) to sz - 1 do  
             arr := (List.append (!arr) [[[j;0];[k;0];[(mult*(j + k)) mod sz; 1]];  
                             [[j;1];[k;1];[(mult*(j+k)) mod sz; 2]];  
                             [[j;2];[k;2];[(mult*(j+k)) mod sz; 0]]])  
         done  
     done;  
     !arr  

As an example,  let's try it for the case where n = 9.

 let l = construct_for_residue_3 9;;  

which gives

 [[[0; 0]; [0; 1]; [0; 2]]; [[1; 0]; [1; 1]; [1; 2]]; [[2; 0]; [2; 1]; [2; 2]];                         [[0; 0]; [1; 0]; [2; 1]]; [[0; 1]; [1; 1]; [2; 2]]; [[0; 2]; [1; 2]; [2; 0]];                         
   [[0; 0]; [2; 0]; [1; 1]]; [[0; 1]; [2; 1]; [1; 2]]; [[0; 2]; [2; 2]; [1; 0]];   
   [[1; 0]; [2; 0]; [0; 1]]; [[1; 1]; [2; 1]; [0; 2]]; [[1; 2]; [2; 2]; [0; 0]]]  

We have a List of 12 elements, where each element contains 3 sublists of ordered number pairs. e.g. the first element is

[[0; 0]; [0; 1]; [0; 2]];

Note that any two number pairs [i; j], [k;l] (where i,k is in range [0,2] and j,l is in range [0,2]),  exist in exactly one sublist.

Case 2:\(n\cong1\mod6\)
This case is slightly more complicated, but the approach is the same. The source code is given below:

 let construct_for_residue_1 v =  
     let res = (v - 1) / 6 in  
     let arr = ref [] in  
     for i = 0 to res - 1 do  
         arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
     done;  
     for j = 0 to 2 * res - 2 do  
         for k = j + 1 to 2 * res - 1 do  
             let z = qgf j k (2 * res) in  
             arr := (List.append (!arr) [[[j;0];[k;0];[z;1]];  
                             [[j;1];[k;1];[z;2]];  
                             [[j;2];[k;2];[z;0]]])  
         done;  
     done;  
     for l = 0 to res - 1 do  
         arr := (List.append (!arr) [[[(-1);(-1)];[l+res;0];[l;1]];  
                         [[(-1);(-1)];[l+res;1];[l;2]];  
                         [[(-1);(-1)];[l+res;2];[l;0]]])  
     done;  
     !arr  

Here, the [i,j] pairs range over different integer values but the triples are the same, i.e. a triple contains three [i,j] pairs. And every two pairs exist in exactly one triple.

For example lets try it with n = 7.

 let l = construct_for_residue_1 7;;  

This gives

  [[[0; 0]; [0; 1]; [0; 2]]; [[0; 0]; [1; 0]; [1; 1]]; [[0; 1]; [1; 1]; [1; 2]];

  [[0; 2]; [1; 2]; [1; 0]]; [[-1; -1]; [1; 0]; [0; 1]]; [[-1; -1]; [1; 1]; [0; 2]];                       
   [[-1; -1]; [1; 2]; [0; 0]]]  

Here, we can see 7 triples and it is easy to check any [[i;j];[k;l]] pair exists in one triple. We can map pairs to the b's in the example at the top of the page. For example

[0;0] ->   b1
[0;1] ->   b2
[0;2] ->   b3
[1;0] ->   b4
[1;1] ->   b5
[1;2] ->   b6
[-1;-1] -> b7

A good question at this point might be, "What's the point in this?", to which the answer is that steiner systems and combinatorial designs in general, have a wide range of connections to different areas of mathematics -  for example:

Projective and Affine Geometry (in fact the above STS(7) steiner system can be considered a projective plane);

Coding theory (e.g. Golay codes);

Group theory (Not going to go into it here, but the above STS(7) has 168 symmetries, and its group of
symmetries is isomorphic to PGL(3,2));

Experimental designs.

It's also worth noting that there are many more interesting steiner systems, beyond the triple system examples here. Probably the most noteworthy is S(5,8,24). Just a quick recap on what that means: 24 points, grouped into subsets, each of size 8, such that any given 5 points exist in exactly 1 subset. That's the definition, but the interesting point is that this object is immensely symmetric, and has over 140 million symmetries, and what's more, any 5 points can be mapped on to any other five points and still preserve symmetry.

source:

1:  (* steiner triple systems *)  
2:  open Core.Std  
3:  exception Impossible_construction  
4:    
5:    
6:  (* construct STS for vertex size 3 modulo 6 *)  
7:  let construct_for_residue_3 v =  
8:      let res = (v - 3) / 6 in  
9:      let arr = ref [] in  
10:      for i = 1 to (2 * res + 1) do  
11:          arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
12:      done;  
13:        
14:      for j = 1 to (2 * res + 1) do  
15:          for k = (j+1) to (2 * res + 1) do  
16:              arr := (List.append (!arr) [[[j;1];[k;1];[(j + k) / 2; 1]];  
17:                              [[j;2];[k;2];[(j+k) / 2; 2]];  
18:                              [[j;0];[k;0];[(j+k) / 2;0]]])  
19:          done  
20:      done;  
21:      !arr  
22:    
23:  (* quasigroup function -- just a function used for Skolem's construction *)  
24:  let qgf left right max =  
25:      let z = (left + right) mod max in  
26:      if z mod 2 = 0 then  
27:          z / 2  
28:      else  
29:          (z + max - 1) / 2  
30:    
31:  (* Skolem's construction - for vertex size 1 modulo 6 *)  
32:  let construct_for_residue_1 v =  
33:      let res = (v - 1) / 6 in  
34:      let arr = ref [] in  
35:      for i = 0 to res - 1 do  
36:          arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
37:      done;  
38:    
39:      for j = 0 to 2 * res - 2 do  
40:          for k = j + 1 to 2 * res - 1 do  
41:              let z = qgf j k (2 * res) in  
42:              arr := (List.append (!arr) [[[j;0];[k;0];[z;1]];  
43:                              [[j;1];[k;1];[z;2]];  
44:                              [[j;2];[k;2];[z;0]]])  
45:          done;  
46:      done;  
47:    
48:      for l = 0 to res - 1 do  
49:          arr := (List.append (!arr) [[[(-1);(-1)];[l+res;0];[l;1]];  
50:                          [[(-1);(-1)];[l+res;1];[l;2]];  
51:                          [[(-1);(-1)];[l+res;2];[l;0]]])  
52:      done;  
53:    
54:      !arr  
55:    
56:  (* find all blocks intersecting the element *)  
57:  let filter_element triplesystem e =  
58:      List.filter triplesystem (fun i -> List.mem i e)