Calculating the e number using Raku

Question

I'm trying to calculate the e constant (AKA Euler's Number) by calculating the formula

In order to calculate the factorial and division in one shot, I wrote this:

my @e = 1, { state $a=1; 1 / ($_ * $a++) } ... *;
say reduce  * + * , @e[^10];

But it didn't work out. How to do it correctly?

Answer 1

I analyze your code in the section Analyzing your code. Before that I present a couple fun sections of bonus material.

One liner One letter¹

say e; # 2.718281828459045

"A treatise on multiple ways"²

Click the above link to see Damian Conway's extraordinary article on computing e in Raku.

The article is a lot of fun (after all, it's Damian). It's a very understandable discussion of computing e. And it's a homage to Raku's bicarbonate reincarnation of the TIMTOWTDI philosophy espoused by Larry Wall.³

As an appetizer, here's a quote from about halfway through the article:

Given that these efficient methods all work the same way—by summing (an initial subset of) an infinite series of terms—maybe it would be better if we had a function to do that for us. And it would certainly be better if the function could work out by itself exactly how much of that initial subset of the series it actually needs to include in order to produce an accurate answer...rather than requiring us to manually comb through the results of multiple trials to discover that.

And, as so often in Raku, it’s surprisingly easy to build just what we need:

sub Σ (Unary $block --> Numeric) {
  (0..∞).map($block).produce(&[+]).&converge
}

Analyzing your code

Here's the first line, generating the series:

my @e = 1, { state $a=1; 1 / ($_ * $a++) } ... *;

The closure ({ code goes here }) computes a term. A closure has a signature, either implicit or explicit, that determines how many arguments it will accept. In this case there's no explicit signature. The use of $_ (the "topic" variable) results in an implicit signature that requires one argument that's bound to $_.

The sequence operator (...) repeatedly calls the closure on its left, passing the previous term as the closure's argument, to lazily build a series of terms until the endpoint on its right, which in this case is *, shorthand for Inf aka infinity.

The topic in the first call to the closure is 1. So the closure computes and returns 1 / (1 * 1) yielding the first two terms in the series as 1, 1/1.

The topic in the second call is the value of the previous one, 1/1, i.e. 1 again. So the closure computes and returns 1 / (1 * 2), extending the series to 1, 1/1, 1/2. It all looks good.

The next closure computes 1 / (1/2 * 3) which is 0.666667. That term should be 1 / (1 * 2 * 3). Oops.

Making your code match the formula

Your code is supposed to match the formula:

In this formula, each term is computed based on its position in the series. The kth term in the series (where k=0 for the first 1) is just factorial k's reciprocal.

(So it's got nothing to do with the value of the prior term. Thus $_, which receives the value of the prior term, shouldn't be used in the closure.)

Let's create a factorial postfix operator:

sub postfix:<!> (\k) { [×] 1 .. k }

(× is an infix multiplication operator, a nicer looking Unicode alias of the usual ASCII infix *.)

That's shorthand for:

sub postfix:<!> (\k) { 1 × 2 × 3 × .... × k }

(I've used pseudo metasyntactic notation inside the braces to denote the idea of adding or subtracting as many terms as required.

More generally, putting an infix operator op in square brackets at the start of an expression forms a composite prefix operator that is the equivalent of reduce with => &[op],. See Reduction metaoperator for more info.

Now we can rewrite the closure to use the new factorial postfix operator:

my @e = 1, { state $a=1; 1 / $a++! } ... *;

Bingo. This produces the right series.

... until it doesn't, for a different reason. The next problem is numeric accuracy. But let's deal with that in the next section.

A one liner derived from your code

Maybe compress the three lines down to one:

say [+] .[^10] given 1, { 1 / [×] 1 .. ++$ } ... Inf

.[^10] applies to the topic, which is set by the given. (^10 is shorthand for 0..9, so the above code computes the sum of the first ten terms in the series.)

I've eliminated the $a from the closure computing the next term. A lone $ is the same as (state $), an anonynous state scalar. I made it a pre-increment instead of post-increment to achieve the same effect as you did by initializing $a to 1.

We're now left with the final (big!) problem, pointed out by you in a comment below.

Provided neither of its operands is a Num (a float, and thus approximate), the / operator normally returns a 100% accurate Rat (a limited precision rational). But if the denominator of the result exceeds 64 bits then that result is converted to a Num -- which trades performance for accuracy, a tradeoff we don't want to make. We need to take that into account.

To specify unlimited precision as well as 100% accuracy, simply coerce the operation to use FatRats. To do this correctly, just make (at least) one of the operands be a FatRat (and none others be a Num):

say [+] .[^500] given 1, { 1.FatRat / [×] 1 .. ++$ } ... Inf

I've verified this to 500 decimal digits. I expect it to remain accurate until the program crashes due to exceeding some limit of the Raku language or Rakudo compiler. (See my answer to Cannot unbox 65536 bit wide bigint into native integer for some discussion of that.)

Footnotes

¹ Raku has a few important mathematical constants built in, including e, i, and pi (and its alias π). Thus one can write Euler's Identity in Raku somewhat like it looks in math books. With credit to RosettaCode's Raku entry for Euler's Identity:

# There's an invisible character between <> and i⁢π character pairs!
sub infix:<⁢> (\left, \right) is tighter(&infix:<**>) { left * right };

# Raku doesn't have built in symbolic math so use approximate equal 
say e**i⁢π + 1 ≅ 0; # True

² Damian's article is a must read. But it's just one of several admirable treatments that are among the 100+ matches for a google for 'raku "euler's number"'.

³ See TIMTOWTDI vs TSBO-APOO-OWTDI for one of the more balanced views of TIMTOWTDI written by a fan of python. But there are downsides to taking TIMTOWTDI too far. To reflect this latter "danger", the Perl community coined the humorously long, unreadable, and understated TIMTOWTDIBSCINABTE -- There Is More Than One Way To Do It But Sometimes Consistency Is Not A Bad Thing Either, pronounced "Tim Toady Bicarbonate". Strangely enough, Larry applied bicarbonate to Raku's design and Damian applies it to computing e in Raku.

Answer 2

There is fractions in $_. Thus you need 1 / (1/$_ * $a++) or rather $_ /$a++.

By Raku you could do this calculation step by step

1.FatRat,1,2,3 ... *   #1 1 2 3 4 5 6 7 8 9 ...
andthen .produce: &[*] #1 1 2 6 24 120 720 5040 40320 362880
andthen .map: 1/*      #1 1 1/2 1/6 1/24 1/120 1/720 1/5040 1/40320 1/362880 ...
andthen .produce: &[+] #1 2 2.5 2.666667 2.708333 2.716667 2.718056 2.718254 2.718279 2.718282 ...
andthen .[50].say      #2.71828182845904523536028747135266249775724709369995957496696762772