code-generation lazy-evaluation raku eulers-number

Calculating the e number using Raku

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?

Solution

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

Super short solution¹

say e; # 2.718281828459045

Existing solutions (longer than a single character)

Though they have tried valiantly², no one will ever beat Damian Conway at writing e-is-for-extraordinarily entertaining technical articles, and his To compute a constant of calculus; (A treatise on multiple ways) is a classic e-is-for-exemplar.

Putting aside the fun he shares, and focusing for a moment on an appetizer of his article's technical substance, the following snippet is a quote and one-liner from about halfway through his 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 }

Using Raku's `...` sequence operator

You very reasonably decided to use Raku's ... sequence operator. (It's a great choice for... a sequence!) But something in your formulation doesn't add up, and we need to analyze what went wrong.

The way sequences work is natural when described in English -- "give me the next value each time I want one until we're all done" -- but I think I need to fairly exhaustively explain the process in this section as preparation for analyzing your code.

Feel free to initially skip or at least skim this section, and then only return to read this exhaustive section if something still isn't adding up.

Let's deal with how and when a sequence concludes before anything starts.

A sequence concludes when the tentative next value of the sequence matches the specified sequence endpoint. For this case the endpoint is infinite (Inf or *) so none of this matters, so feel free to skip this sub-section.

There are three cases leading to conclusion of a sequence:

The endpoint matches and includes whatever value matches it (i.e. the operator is written as a plain ...). In this case the tentative next value is included as the final value in the sequence.
The endpoint matches but excludes whatever value matches it (i.e. the right hand side of the sequence operator ends with a ^, i.e. it's written ...^ endpoint). In this case the tentative next value is not included in the final sequence.
The tentative next value exceeds the endpoint. In this case it is not included.

The sequence operator (...) yields a sequence as follows:

The values listed at the start of the list on its left, if any, are the initial values of the sequence. For example, the initial values in the following sequence are 1, 2, and 3:
```
1, 2, 3, { generator code goes here } ... *
```
If/when any of these values match or exceed the endpoint, or there are no values left and there is no generator routine or closure listed to generate more, the sequence is concluded. Otherwise, if/when the initial values are exhausted, then remaining values in the sequence are generated by the generator routine/closure that is the last value in the list to the left of the .... The generation process is explained in the next step.
Values are lazily generated by calling the routine/closure on the operator's immediate left. This is done once for each new value in the sequence, until the conclusion of the sequence as described earlier.

Each time the routine/closure is called, it is passed N positional arguments:

The arguments are the previous N terms in the sequence (either literally listed or generated).
N is the maximum of the number of previous terms already available, and the routine's/closure's count (the maximum number of positional arguments it will accept at a time).

A routine/closure has a signature, either implicit or explicit, which determines its count. For your closure there is no explicit signature, so it gets the default implicit one (which takes a single optional positional argument named $_, the "topic" variable).

So for your closure N is 1. (First, there's already a 1 literally listed at the start of the sequence on the left of the ..., so there's at least 1 argument available to pass. Second, the closure will accept (a maximum of) 1 argument each time, so it will accept the 1 and then accept whatever was the value it last generated as the next value it's given as $_.)

Analyzing your code

Here's your first line, generating the series:

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

The closure ({ code goes here }) computes the next term in the sequence.

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, then 1/1.

The topic in the second call is the value of the previous term, 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.

My way based on your way (#1)

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. Back to the drawing board!

My way based on your way (#2)

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 anonymous 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 both of its operands are 100% accurate number types (so they're not the fallback generic Num type, which is a float, and thus approximate), the / operator normally returns a 100% accurate Rat (a limited precision rational).

But if the denominator of the result overflows 64 bits then, by default, that result is converted to a Num (a float). This maintains high performance but at the cost of approximation.

We need to control this tradeoff.

The ideal way to control overflow behavior, as locally or globally as you want in some code, is by using $*RAT-OVERFLOW.

But if you just want to control an individual division, then just make sure that all values and (values returned by function calls) are individually integers or rationals, and then make the numerator a FatRat, eg:

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

All the 1s are integers; the ++$ will always yield an integer; therefore the 1 .. ++$ will too; therefore the [×] multiplication will too; and the numerator 1 is converted to a FatRat in order to ensure the value of the overall expression into a FatRat.

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 π). So 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 (or maybe that should be discredit for using invisible characters):

sub infix:<⁢> (\left, \right) is tighter(&infix:<**>) { left * right };
#          ⬑ Invisible character between `<` and `>`.

# Raku doesn't have built in symbolic math so use approximate equal 
say e**i⁢π + 1 ≅ 0; # True
#       ⬑ Invisible character between `i` and `π`.

² Damian's article is a must read, as always, as a tour de force homage to Raku's bicarbonate reincarnation of the TIMTOWTDI philosophy espoused by Larry Wall.³ 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.