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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
#################################################################################
# This script requires Math::GMPq, Math::GMPz, and Math::MPFR. #
# It calculates the euler number e (2.7182818...), correct to $ARGV[0] bits. #
# The calculated value is displayed unless $ARGV[1] is both provided and false. #
# With each iteration of the for{} loop (below) we get closer and closer to #
# the actual value of e. Furthermore, with successive iterations of the for{} #
# loop, the values alternate between "less than e" and "greater than e". #
# Hence the actual (irrational) value of e is always between the values #
# calculated by successive iterations of the for{} loop. #
# #
# Of course, the simplest and most efficient way to get the value of e, to #
# $ARGV[0] bits is simply to do: #
# Rmpfr_exp($rop, Math::MPFR->new(1), MPFR_RNDN) #
# where $rop is a $ARGV[0]-bit precision Math::MPFR object. #
# But doing it that way is a bit less interesting. #
# #
# The same for{} loop can be also used to calculate the exact probabilities of #
# "winning" at a simplistic solitaire-type card game. See demos/solitaire.p #
# in the Math::GMPz source distro. #
#################################################################################
use strict;
use warnings;
use Math::GMPz qw(:mpz);
use Math::GMPq qw(:mpq);
use Math::MPFR qw(:mpfr);
die "Usage: perl euler.pl bits [True|False]" unless @ARGV;
my $bits = shift;
Rmpfr_set_default_prec($bits);
my $display_value;
$display_value = defined($ARGV[0]) ? shift : 1;
#################################################################################
# For the sanity checks (below), set $e_big_p to e, correct to $bits+100 bits. #
# Then convert $e_big_p exactly to a rational, $e_q (a Math::GMPq object). #
#################################################################################
my $e_q = Math::GMPq->new(); #
my $e_big_p = Rmpfr_init2($bits + 100); #
Rmpfr_exp($e_big_p, Math::MPFR->new(1), MPFR_RNDN); #
Rmpfr_get_q($e_q, $e_big_p); #
#################################################################################
#################################################################
# Create some variables, and assign some initial values #
#################################################################
my $first = Math::GMPz->new(1); #
my $second = Math::GMPz->new(0); #
my $current_items = 2; #
my $factorial = Math::GMPz->new(1); #
my $e_check = Math::GMPq->new(); #
my ($e, $e_first_fr, $e_second_fr) = (Math::MPFR->new(), #
Math::MPFR->new(), #
Math::MPFR->new(), #
); #
my $chance; # becomes a MATH::GMPz object on assignment #
my $e_first = Math::GMPq->new(3); #
my $e_second = Math::GMPq->new(); #
my $count = 0; #
my $t = Math::GMPq->new(); #
my $save = Math::GMPq->new(4); #
#################################################################
#################################################################
# Set $e to a $bits-bit approximation of the euler number, #
# rounded to nearest. #
# This should exactly equal the number that we calculate. #
#################################################################
Rmpfr_exp($e, Math::MPFR->new(1), MPFR_RNDN); #
#################################################################
#########################
# Do the calculations #
#########################
for(;;) {
$count++;
Rmpz_mul_ui($factorial, $factorial, $current_items); #$factorial *= $current_items;
$chance = ($current_items - 1) * ($first + $second);
#########################################################
# In this block we just perform some sanity checks. #
# This block plays no part in the calculation of the #
# actual value. #
# Assign the calculated rational value to $e_check. #
# Check that for every 2nd iteration, $e_check > $e_q #
# and that for every other iteration, $e_check < $e_q #
# Also check that, with each iteration, we get closer #
# to the value of e (ie closer to the value of $e_q) #
#########################################################
Rmpq_set_num($e_check, $factorial); #
Rmpq_set_den($e_check, $chance); #
Rmpq_canonicalize($e_check); # gcd(num, den) == 1 #
if($count % 2) { #
unless($e_check < $e_q) {die "$count: >="} #
} #
else { #
unless($e_check > $e_q) {die "$count: <="} #
} #
Rmpq_sub($t, $e_q, $e_check); #
if(abs($t) < $save) {Rmpq_set($save, abs($t))} #
else {die "$count: No closer to e"} #
#########################################################
Rmpq_set_num($e_second, $factorial);
Rmpq_set_den($e_second, $chance);
Rmpq_canonicalize($e_second); # gcd(num, den) == 1
Rmpfr_set_q($e_first_fr, $e_first, MPFR_RNDN);
Rmpfr_set_q($e_second_fr, $e_second, MPFR_RNDN);
#########################################################
# Exit the loop when $e_first_fr == $e_second_fr #
# as this equivalence indicates that both variables #
# contain the euler number, correct to $bits bits. #
#########################################################
last if Rmpfr_equal_p($e_first_fr, $e_second_fr); #
#########################################################
Rmpz_set($first, $second);
Rmpz_set($second, $chance);
Rmpq_set($e_first, $e_second);
$current_items++;
}
if($e == $e_first_fr) {
print "Iterations: $count ok\n";
if($display_value) {print "$e_first_fr\n"}
}
else {print print "Iterations: $count not ok\n$e\n$e_first_fr\n"}
__END__
The sequence:
With 1st iteration, e = 2! divided by 1 (ie divided by 1 * (1 + 0 ))
With 2nd iteration, e = 3! divided by 2 (ie divided by 2 * (0 + 1 ))
With 3rd iteration, e = 4! divided by 9 (ie divided by 3 * (1 + 2 ))
With 4th iteration, e = 5! divided by 44 (ie divided by 4 * (2 + 9 ))
With 5th iteration, e = 6! divided by 265 (ie divided by 5 * (9 + 44 ))
With 6th iteration, e = 7! divided by 1854 (ie divided by 6 * (44 + 265 ))
With 7th iteration, e = 8! divided by 14833 (ie divided by 7 * (265 + 1854 ))
With 8th iteration, e = 9! divided by 133496 (ie divided by 8 * (1854 + 14833 ))
With 9th iteration, e = 10! divided by 1334961 (ie divided by 9 * (14833 + 133496))
|
