PiplusErationals/Pipluserationals.mac at master · T3gu1a/PiplusErationals · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
/************************************************************************************************************

This Maxima package provides two functions "pirationals" and "pipluserationals" that compute very good rational approximations to Pi and Pi+e, respectively. The output is a list of pairs of integers where the first component represents the numerator and the second represents the denominator.

Eg. Once the package is load to a Maxima environement (like wxMaxima), the code can be used as follows:

(%i1) pirationals(1,500);

(%o1) [[22,7],[44,14],[355,113]]

(%i2) pipluserationals(1,1000);

(%o2) [[920,157]]

920/157, and {22/7, 355/113} are particular convergents of the continued fraction of pi+e, and pi respectively.
***************************************************************************************************************/

Copyright():=block([],
    print("Copyright (c) May 2020 Bertrand Teguia T."),
    print("https://www.bertrandteguia.com"),
    print("bteguia@mathematik.uni-kassel.de"),
    print("University of Kassel, Germany"),
    print("")
)$

/* The code for Pi */

sumpartpiestim[0]:4$

sumpartpiestim[n]:=sumpartpiestim[n-1]+binomial(2*n,n)/4^(n-1)/(2*n+1)$

piestim(q,j):=block([],
    1/q-sumpartpiestim[j]/2/(j+1)
)$

pirationals(minnum,maxnum):=catch(block(
    [lobound,upbound,Listrat,n,j,q,evlobound,evpiestimn,evupbound,evpiestimq,k],
    if(maxnum<6) then throw([]),
    lobound:binomial(2*(n+1),n+1)/2/(n+1)/4^n,
    upbound:binomial(2*(n+1),n+1)/(2*n+1)/4^n,
    Listrat:[],
    for j:max(minnum,6) thru maxnum do (
        evlobound:ev(lobound,n=j),
        evpiestimn:piestim(q,j),
        evupbound:ev(upbound,n=j),
        k:floor(j/3),
        evpiestimq:ev(evpiestimn,q=k),
        while(evpiestimq<=evupbound and k>1) do (
            if(evpiestimq>=evlobound) then Listrat:endcons([j+1,k],Listrat),
            k:k-1,
            evpiestimq:ev(evpiestimn,q=k)
        )

    ),
    Listrat

))$

/* The code for Pi+e */

sumparteestim[0]:1$

sumparteestim[n]:=sumparteestim[n-1]+1/n!$

pipluseestim(q,j):=block([],
    1/q-(sumpartpiestim[j]/2+sumparteestim[j])/(j+1)
)$


pipluserationals(minnum,maxnum):=catch(block(
    [lobound,upbound,Listrat,n,j,q,evlobound,evpipluseestimj,evupbound,
    evpipluseestimq,k],
    if(maxnum<10) then throw([]),
    lobound:binomial(2*(n+1),n+1)/2/(n+1)/4^n,
    upbound:binomial(2*(n+1),n+1)/(2*n+1)/4^n+1/(n+1)!/n,
    Listrat:[],
    for j:max(minnum,10) thru maxnum do (
        evlobound:ev(lobound,n=j),
        evpipluseestimj:pipluseestim(q,j),
        evupbound:ev(upbound,n=j),
        k:floor(j/5),
        evpipluseestimq:ev(evpipluseestimj,q=k),
        while(evpipluseestimq<=evupbound and k>1) do (
            if(evpipluseestimq>=evlobound) then Listrat:endcons([j+1,k],Listrat),
            k:k-1,
            evpipluseestimq:ev(evpipluseestimj,q=k)
        )
    ),
    Listrat
))$

Copyright();