# Following Bailey, Borwein, Plouffe (p. 3), we set up the integral for
# Pi with undetermined coefficients:
> g := (b1 + b2*y + b3*y^2 + b4*y^3 + b5*y^4 + b6*y^5 +
> b7*y^6)/(1-y^8/16);

                             2       3       4       5       6
             b1 + b2 y + b3 y  + b4 y  + b5 y  + b6 y  + b7 y
        g := -------------------------------------------------
                                          8
                                1 - 1/16 y

> pib := int(g, y=0..1);

        1/2          1/2            1/2              1/2
pib := 2    %1 b5 + 2    %2 b5 + 2 2    %1 b7 + 1/2 2    %1 b3

            1/2              1/2            1/2
     - 1/2 2    %2 b3 + 1/4 2    %2 b1 - 2 2    %2 b7

            1/2
     + 1/4 2    %1 b1 + 1/2 ln(5) b4 + 1/4 ln(3) b2 - 1/4 ln(5) b3

     + ln(5) b7 - 1/2 ln(3) b4 + ln(3) b6 - 1/2 ln(5) b5

     + 1/8 ln(5) b1 + 2 arctan(2) b6 - 2 arctan(2) b7 - arctan(2) b5

     + 1/2 arctan(2) b3 + 1/4 arctan(2) b1 - 1/2 arctan(2) b2 - Pi b6

     + 1/4 Pi b2

                   1/2
%1 := arctanh(1/2 2   )

                  1/2
%2 := arctan(1/2 2   )

# The above integral shall be an integral multiple of pi for selected
# bi's.  We construct a lattice with approximate coefficients for the
# bi's shifted by 2^k:
> latt := proc(k)
>    local approxpi, v1, v2, v3, v4, v5, v6, v7, v8;
>    approxpi := evalf(pib, k+2);
>    print(approxpi);
>    v1 := [1,0,0,0,0,0,0,0,trunc(evalf(coeff(approxpi,b1)*10^k,k))];
>    v2 := [0,1,0,0,0,0,0,0,trunc(evalf(coeff(approxpi,b2)*10^k,k))];
>    v3 := [0,0,1,0,0,0,0,0,trunc(evalf(coeff(approxpi,b3)*10^k,k))];
>    v4 := [0,0,0,1,0,0,0,0,trunc(evalf(coeff(approxpi,b4)*10^k,k))];
>    v5 := [0,0,0,0,1,0,0,0,trunc(evalf(coeff(approxpi,b5)*10^k,k))];
>    v6 := [0,0,0,0,0,1,0,0,trunc(evalf(coeff(approxpi,b6)*10^k,k))];
>    v7 := [0,0,0,0,0,0,1,0,trunc(evalf(coeff(approxpi,b7)*10^k,k))];
>    v8 := [0,0,0,0,0,0,0,1,trunc(evalf(Pi*10^k,k+1))];
>    RETURN([v1,v2,v3,v4,v5,v6,v7,v8]);
> end;
> 

latt := proc(k)
local approxpi, v1, v2, v3, v4, v5, v6, v7, v8;
    approxpi := evalf(pib, k + 2);
    print(approxpi);
    v1 := [1, 0, 0, 0, 0, 0, 0, 0,
        trunc(evalf(coeff(approxpi, b1)*10^k, k))];
    v2 := [0, 1, 0, 0, 0, 0, 0, 0,
        trunc(evalf(coeff(approxpi, b2)*10^k, k))];
    v3 := [0, 0, 1, 0, 0, 0, 0, 0,
        trunc(evalf(coeff(approxpi, b3)*10^k, k))];
    v4 := [0, 0, 0, 1, 0, 0, 0, 0,
        trunc(evalf(coeff(approxpi, b4)*10^k, k))];
    v5 := [0, 0, 0, 0, 1, 0, 0, 0,
        trunc(evalf(coeff(approxpi, b5)*10^k, k))];
    v6 := [0, 0, 0, 0, 0, 1, 0, 0,
        trunc(evalf(coeff(approxpi, b6)*10^k, k))];
    v7 := [0, 0, 0, 0, 0, 0, 1, 0,
        trunc(evalf(coeff(approxpi, b7)*10^k, k))];
    v8 := [0, 0, 0, 0, 0, 0, 0, 1, trunc(evalf(Pi*10^k, k))];
    RETURN([v1, v2, v3, v4, v5, v6, v7, v8])
end

> L3 := latt(30);

.20500255763642353394415033621848 b5

     + .1472019346726350271955981288534 b7

     + .33923024524519907251558532577876 b3

     + 1.0071844764146762286447601474504 b1

     + .25541281188299534160275704815185 b4

     + .50647687666743048095593942496125 b2

     + .1713170706664974589667327740000 b6


L3 := [[1, 0, 0, 0, 0, 0, 0, 0, 1007184476414676228644760147450],

    [0, 1, 0, 0, 0, 0, 0, 0, 506476876667430480955939424961],

    [0, 0, 1, 0, 0, 0, 0, 0, 339230245245199072515585325779],

    [0, 0, 0, 1, 0, 0, 0, 0, 255412811882995341602757048152],

    [0, 0, 0, 0, 1, 0, 0, 0, 205002557636423533944150336218],

    [0, 0, 0, 0, 0, 1, 0, 0, 171317070666497458966732774000],

    [0, 0, 0, 0, 0, 0, 1, 0, 147201934672635027195598128853],

    [0, 0, 0, 0, 0, 0, 0, 1, 3141592653589793238462643383280]]

# Next, we perform lattice reduction on L3 to obtain canditates for the
# bi's:
> readlib(lattice):
> lattice(L3);

[[-4, 0, 0, 2, 1, 1, 0, 1, 2], [0, -8, -4, -4, 0, 0, 1, 2, 1],

    [1744, 352, -9056, 12300, 2277, -45740, 378, 1690, 12068],

    [-116, 15655, -34105, 8724, -33082, 10163, 19406, 1182, 1910],

    [26048, 22830, -47673, 3474, 24875, 13607, -9911, -9067, 33918],

    [15089, 7207, -44923, 29231, 37260, 356, 23394, -7072, -14326],

    [-6870, 606, -13075, -1863, -2657, -17305, -64791, 7821, -5804],

    [43110, -34168, -3663, 67855, -15158, 48067, 4404, -15272, 9543]]

# The first short vector is Bailey et al., the second is ours.  We
# conclude by proving our formula:
> subs({b1=0, b2=8, b3=4, b4=4, b5=0, b6=0, b7=-1}, pib);

                                 2 Pi