Wait the light to fall

计算微积分常数

焉知非鱼

𝑒 is for economics

第 21 周挑战的第一个任务是一个非常古老的任务:找出(最终被称为)欧拉数。

故事要追溯到 1683 年, 雅各布·伯努利(Jacob Bernoulli),以及他对高利贷数学的研究。例如,假设你提供一年期的 1000 美元贷款, 利息为 100%, 按年还款。很显然,到了年底, 签收客户要把 1000 美元还给你,再加上($1000×100%)的利息……所以你现在有 2000 美元。我能说什么呢?甜言蜜语!

但是, 你又会想: 假设每 6 个月收取一次利息呢? 在这种情况下, 6 个月后, 他们已经欠你 500 美元的利息($1000 x 100% × 6∕12), 你可以立即开始收取利息! 所以在最后 6 个月后, 他们现在已经欠你原来的 1000 美元加上前 6 个月的利息, 加上后 6 个月的利息, 再加上前 6 个月的利息上的后 6 个月的利息: $1000 + ($1000 × 50%) + ($1000 × 50%) + ($1000 × 50% × 50%)。 也就是 2250 美元。

这更香了。

当然, 用数学方法计算出他们欠你的钱更容易: $1000 × (1 + ½)²

附加的 ½ 来自于每半年收取一半的年利息。2 的幂来自于每年收取两次利息。

但为什么要止步于此呢? 为什么不按月收取利息呢? 那么你就会拿回原来的 1000 美元, 加上第一个月的 83.33 美元的利息(即: 1000 美元的 1∕12), 加上第二个月的 90.28 的利息(即: 1000 美元的 1∕12 加上 83.33 美元), 加上第三个月的 97.80 的利息(即: 1000 美元的 1∕12 加上 83.33 美元, 加上 90.28 美元), 等等。 换句话说: 1000 x (1 + ¹⁄₁₂)¹² 或 2613.04 美元。 非常好。

如果我们按天收取利息, 我们会得到 2714.57 美元(即 1000 x (1 + ¹⁄₃₆₅)³⁶⁵)。

如果我们按小时收取利息, 我们会得到 2718 美元(即 1000 x (1 + ¹⁄₈₇₆₀)⁸⁷⁶⁰)。

如果我们按分钟收取利息, 我们会得到 2718.27 美元(即 1000 x (1 + ¹⁄₅₂₅₆₀₀)⁵²⁵⁶⁰⁰)。

伯努利接着问了一个显而易见的问题:你能榨取出多少钱? 如果你连续收取利息, 你能收拿回多少钱呢?

换句话说,如果你每时每刻都收取和积累利息,你最初的1000美元的倍数是多少?或者说,从数学上讲, 当 x→∞ 时, (1 + ¹⁄ₓ)ᕽ 的极限是多少?

这个问题的答案是超自然常数 𝑒:

2.7182818284590452353602874713526624977572470936999595749669...

它是以欧拉而不是以雅各布伯努利的名字命名的, 因为:

(a) 欧拉是第一个在发表的科学论文中明确使用它的人, 以及 (b) 如果我们用欧拉的名字来命名所有的东西, 确实会让一切变得简单很多。

𝑒 is for elapsed

所以, 鉴于 𝑒 是 (1 + ¹⁄ₓ)ᕽ 的极限, 对于越来越大的 x 值, 我们只需用下面的代码就可以计算出逐渐精确的常数近似值:

for 1, 2, 4 ... ∞ -> \x {
    say (1 + x⁻¹) ** x
}

…它产生如下输出:

2
2.25
2.441406
2.565784514
2.6379284973666
2.676990129378183
2.697344952565099
2.7077390196880207
2.7129916242534344
2.7156320001689913
2.7169557294664357
2.7176184823368796
2.7179500811896657
2.718115936265797
2.7181988777219708
2.718240351930294
2.7182610899046034
2.7182714591093062
2.718276643766046
2.718279236108013
⋮

这是一种简单的解决办法, 但不是很好。经过 20 次迭代, 在 N=524288 的情况下, 仍然只能精确到小数点后五位(2.718279236108013)。

更重要的是, 这个方法很慢。有多慢呢? 让我们抡起一个集成的计时机制来告诉我们:

# Evaluate the block and add timing info to the result...
sub prefix:<timed> (Code $block --> Any) {
    # Represent added timing info...
    role Timed {
        # Store the actual timing...
        has Duration $.timing;

        # When stringified, append the timing...
        my constant FORMAT = "%-20s  [%.3f𝑠]";
        method Str  { sprintf FORMAT, callsame, $.timing }
        method gist { self.Str }
    }

    # Evaluate the block, mixing timing info into result...
    return $block() but Timed(now - ENTER now);
}

在这里,我们创建了一个新的高优先级的前缀运算符:timed。这是一个真正的运算符,尽管它的符号也是一个正确的标识符。如果我们愿意的话, 也可以用 ASCII 或 Unicode 符号来命名它:

# Vaguely like a clock face...
sub prefix:« (<) » (Code $block --> Any) {…}

# Exactly like a stopwatch...
sub prefix:« ⏱ »   (Code $block --> Any) {…}

即使没有花俏的符号,我们仍然需要一个运算符,而不是一个子例程,因为常规子例程调用的优先级太低。如果 timed 被声明为 sub timed, 我们就不能把定时块放到周围的参数列表中,因为对 timed 的调用会试图吞噬一切跟随的参数,然后发现它只允许一个参数:

say timed { expensive-func() }, " finished at ", now;
# Too many positionals passed; expected 1 argument but got 3
#   in sub timed at demo.p6 line 5

timed 成为一个前缀运算符可以让编译器知道它只能接受一个参数,所以像上面这样的内联调用就可以正常工作。

需要注意的是,由于我们今天觉得需要额外的约束,所以我们要利用 Raku 的类型系统,严格地对我们使用的每一个变量和参数进行类型化。当然,由于 Raku 的静态类型化是渐进式的,所以如果我们后来把这些类型声明中的每一个都去掉,代码的工作原理还是完全一样的……只是这样一来,编译器就不能很好地保护我们不受自己的愚蠢影响了。

timed 运算符将一个 Code 对象作为参数,执行它($block()),用额外的定时信息增强结果,并返回增强后的结果。可以是任何类型(--> Any)。

通过将名为 Timed 的通用类组件(一个角色)定义的一些额外的能力混合到返回的对象中,该定时信息被添加到块的结果中。该角色赋予了存储和访问持续时间的能力(has Duration $.timing),以及用于增强定时对象的字符串化和打印方式的方法(方法 Str 和方法 gist)。重载的 .Str 方法调用对象的同一方法的原始版本(callsame)来做实际工作,然后将定时信息附加到它上面,整齐地用 sprintf 格式化。覆盖的 .gist 只是重用了 .Str

最有趣的特点是,当我们将定时信息添加到块的结果中时(but Timed(...)),我们通过从周围子程序被输入的当前瞬间(ENTER now)减去块执行后的瞬间(now)来计算块的执行时间。在任何表达式前加上 ENTER,都会设置一个"phaser"(是的,我们知道:我们是无法治愈的怪胎)。phaser 是一个块或 thunk,当周围的块被输入时就会执行。然后,ENTER 会记住由 thunk 产生的值,并在周围的代码执行时对其进行求值:在本例中,当它试图从 now 的值中减去 ENTER 的值时。

您可以在 Raku 中的任何地方使用同样的技术,作为一种方便的方式来为任何特定的代码块计时。您也可以使用其他 phaser(如 LEAVE,当控制离开周围的块时执行)将整个操作内嵌到一个易于粘贴的片段中。例如,为了给 for 循环的每次迭代计时,你可以在循环块的开始处添加一行:

for @values -> $value {
    LEAVE say "$value took ", (now - ENTER now), " seconds";

    do-something-expensive-with($value);
    do-something-else-protracted-with($value);
    do-the-last-costly-thing-with($value);
}

所以,现在我们可以很容易地计时任何代码块,只需在块前面加上 timed 运算符。当我们在我们的计算循环中这样使用时:

for 1, 2, 4 ... ∞ -> \x {
    say timed { (1 + x⁻¹) ** x }
}

…我们发现,为了提高精确度,我们要等待成倍的时间:

⋮
2.7176184823368796    [0.001𝑠]
2.7179500811896657    [0.002𝑠]
2.718115936265797     [0.011𝑠]
2.7181988777219708    [0.051𝑠]
2.718240351930294     [0.241𝑠]
2.7182610899046034    [1.123𝑠]
2.7182714591093062    [4.802𝑠]
2.718276643766046     [21.468𝑠]
2.718279236108013     [106.543𝑠]
⋮

显然,我们需要一个更好的算法。

𝑒 is for extrema

初等数学的100个大问题中,海因里希·多里表明,正如 (1 + ¹⁄ₓ)ᕽ 是当 x→∞ 时 𝑒 的下界,所以 (1 + ¹⁄ₓ)ᕽ⁺¹ 是 𝑒 的上界。所以,在 N 值相同的情况下,我们可以通过取这两个约束值的平均值,得到一个更好的近似值 𝑒:

for 1, 2, 4 ... ∞ -> \x {
    say timed {
        ½ × sum (1 + x⁻¹)**(x),
                (1 + x⁻¹)**(x+1)
    }
}

它产生如下输出:

3                     [0.001𝑠]
2.8125                [0.000𝑠]
2.746582              [0.000𝑠]
2.72614604607         [0.000𝑠]
2.7203637629093063    [0.000𝑠]
2.7188181001497167    [0.001𝑠]
2.718417960007014     [0.001𝑠]
2.718316125233677     [0.000𝑠]
2.7182904360195543    [0.000𝑠]
2.718283984544156     [0.000𝑠]
2.7182823680062143    [0.000𝑠]
2.718281963411669     [0.001𝑠]
2.718281862205436     [0.005𝑠]
2.7182818368966726    [0.021𝑠]
2.718281830568581     [0.101𝑠]
2.718281828986445     [0.456𝑠]
2.7182818285908974    [2.171𝑠]
2.7182818284920085    [9.628𝑠]
2.718281828467286     [46.717𝑠]
2.718281828461105     [206.456𝑠]
⋮

在撞上指数墙之前, 已经有10 位正确的小数位(2.7182818284920085)。比较好,但还是不够好。

𝑒 is for evaluation

看起来我们需要在一个很宽的数值范围内尝试很多不同的技术。因此,如果能有一个更简单的方法来指定一系列类似的测试,并有一个更好的方法来查看它们的性能如何,那将会很方便。

所以我们要创建一个框架,让我们更简单地测试各种技术,就像这样:

#| Bernoulli's limit
assess -> \x { (1 + x⁻¹)**(x) }

#| Dörrie's bounds
assess -> \x {
    ½ × sum (1 + x⁻¹)**(x),
            (1 + x⁻¹)**(x+1)
}

或者,如果我们愿意的话,可以指定一个更合适的试验值范围,而不仅仅是 1..∞,像这样:

#| Dörrie's bounds
assess -> \x=(1,10,100...10⁶) {
    ½ × sum (1 + x⁻¹)**(x),
            (1 + x⁻¹)**(x+1)
}

无论哪种方式,assess 函数都会计算每一个结果,决定何时放弃缓慢的计算,将结果列表,对其准确性进行颜色编码,并打印出整齐的标签,就像这样:

[1] Bernoulli's limit (x = 1): 2                     [0.001𝑠]
                      (x = 2): 2.25                  [0.000𝑠]
                      (x = 4): 2.441406              [0.000𝑠]
                      (x = 8): 2.565784514           [0.000𝑠]
                     (x = 16): 2.6379284973666       [0.000𝑠]
                     (x = 32): 2.676990129378183     [0.000𝑠]
                     (x = 64): 2.697344952565099     [0.000𝑠]
                    (x = 128): 2.7077390196880207    [0.000𝑠]
                    (x = 256): 2.7129916242534344    [0.000𝑠]
                    (x = 512): 2.7156320001689913    [0.000𝑠]
                   (x = 1024): 2.7169557294664357    [0.000𝑠]
                   (x = 2048): 2.7176184823368796    [0.001𝑠]
                   (x = 4096): 2.7179500811896657    [0.003𝑠]
                   (x = 8192): 2.718115936265797     [0.011𝑠]
                  (x = 16384): 2.7181988777219708    [0.052𝑠]
                  (x = 32768): 2.718240351930294     [0.250𝑠]
                  (x = 65536): 2.7182610899046034    [1.175𝑠]
                 (x = 131072): 2.7182714591093062    [5.337𝑠]



   [4] Dörrie's bounds (x = 1): 3                     [0.001𝑠]
                      (x = 2): 2.8125                [0.000𝑠]
                      (x = 4): 2.746582              [0.000𝑠]
                      (x = 8): 2.72614604607         [0.000𝑠]
                     (x = 16): 2.7203637629093063    [0.000𝑠]
                     (x = 32): 2.7188181001497167    [0.000𝑠]
                     (x = 64): 2.718417960007014     [0.001𝑠]
                    (x = 128): 2.718316125233677     [0.001𝑠]
                    (x = 256): 2.7182904360195543    [0.000𝑠]
                    (x = 512): 2.718283984544156     [0.000𝑠]
                   (x = 1024): 2.7182823680062143    [0.000𝑠]
                   (x = 2048): 2.718281963411669     [0.001𝑠]
                   (x = 4096): 2.718281862205436     [0.004𝑠]
                   (x = 8192): 2.7182818368966726    [0.025𝑠]
                  (x = 16384): 2.718281830568581     [0.109𝑠]
                  (x = 32768): 2.718281828986445     [0.500𝑠]
                  (x = 65536): 2.7182818285908974    [2.369𝑠]
                 (x = 131072): 2.7182818284920085   [10.474𝑠]

为了实现这一切,我们从一个适合存储、审核和格式化任何我们收集的数据的通用类组件(即另一个角色)开始:

# Define a new type: Code that takes exactly one parameter...
subset Unary of Code where *.signature.params == 1;

# Reportable values, assessed against a target...
role Reportable [Cool :$target = ""]
{
    # Track reportable objects and eventually display them...
    state Reportable @reports;
    END @reports».display: @reports».width.max;

    # Add reportable objects to tracking list, and validate...
    submethod TWEAK (Unary :$check = {True}) {
        @reports.push: self;
        $check(self);
    }

    # Each report has a description and value...
    has Str  $.desc  handles width => 'chars';
    has Cool $.value handles *;

    # Display the report in two columns...
    method display (Int $width) {
        printf "%*s: %s\n", $width, $!desc, &.assess-value;
    }

    # Colour-code the value for display...
    method assess-value () {
        # Find leading characters that match the target...
        my Int $correct
            = ($!value ~^ $target).match( /^ \0+/ ).chars;

        # Split the value accordingly...
        $!value.gist ~~ /^ $<good> = (. ** {$correct})
                           $<bad>  = (\S*)
                           $<etc>  = (.*)
                        /;

        # Correct chars: blue; incorrect chars: red...
        use Terminal::ANSIColor;
        return colored($<good>.Str, 'blue')
             ~ colored($<bad>.Str,  'red')
             ~ $<etc>
    }
}

我们首先为内置的 Code 类型(即一般的块、lambda、子例程等类型)创建一个新的子类型(subset Unary),这个新的子类型要求 Code 对象必须只接受一个参数(where *.signature.params == 1)。因为在 Raku 中(几乎)所有的东西都是对象,所以语言很容易提供这些类型的关于值、变量、类型和代码的详细自省方法

接下来我们创建一个可插拔的组件(role Reportable),用于构建生成自组织报表的类。这个角色被参数化为接受一个命名的值(:$target),随后将用于评估被报告的各个值的准确性。该目标值被指定为 Cool 类型,允许它是字符串、数字、布尔值、持续时间、哈希值、列表或任何其他可自动转换为字符串或数字的类型。这个类型约束将确保以后报告的任何值能够被正确评估。:target 参数是可选的;默认 target 是空字符串。

Reportable 角色要为我们自动跟踪所有可报告的对象……然后在程序结束时生成汇总报告。所以我们声明一个共享变量来存放每一个报表(state Reportable @reports),并指定类型,要求数组中的每个元素都能执行 Reportable 角色。为了保证报表最终被打印出来,我们添加一个 END phaser:

END @reports».display: @reports».width.max;

在执行结束时,这个语句调用每个报表的 .display 方法(@reports».display),将任何报表描述的最大宽度(@reports».width.max)传递给它,这样它就可以将报表格式化为两列干净的报表。

为了积累这些报表,我们安排每个 Reportable 对象在被构造时自动将自己添加到共享的 @reports 数组中……通过声明一个 TWEAK 子方法。submethod 是一个类特有的、非继承的方法,适合于指定每个类的初始化器。TWEAK 子方法在对象被创建后会被自动调用:通常是为了调整它的初始值,或者以某种方式测试它的完整性。

在这里,我们同时做了这两件事:让子方法将每个对象添加到报告列表中(@reports.push: self),并应用传递给构造函数的任何 :check 代码($check(self))。这就允许用户在调用构造函数时传递一个任意的测试,并在对象被初始化后将其应用到对象上。我们很快就会看到这有多有用。

每个报告由一个字符串描述和一个可以是任何 Cool 类型的值组成,所以我们需要每个对象属性来存储它们。我们用 has 关键字和一个"点"辅助符号来声明它们,使它们成为公共的:

has Str  $.desc;
has Cool $.value;

我们还需要能够访问描述的字符串宽度。我们可以声明一个方法来提供这些信息:

method width { $.desc.chars }

但是,由于这只是将 Reportable 对象的请求转发到其内部的 $.desc 对象,所以有一个更简单的方法来实现同样的效果:我们只需告诉 $.desc 属性,它应该通过调用自己的 .chars 方法来处理所有对象级对 .width 的调用。像这样:

has Str $.desc handles width => 'chars';

更重要的是,我们还需要能够将方法转发到 $.value 属性上(例如,查询它的 .timing 信息)。但由于我们不知道值在每种情况下可能是什么样的对象(除了通用的是 Cool),我们无法提前知道可能需要转发哪些方法。所以我们只需要告诉 $.value 处理周围对象本身不能处理的东西,就像这样:

has Cool $.value handles *;

现在我们只需要实现 .display 方法,END phaser 将用于输出每个 Reportable 对象。该方法获取描述应该被格式化的宽度,并使用 printf 将其打印出来。它还会打印值,这个值是用 .assess-value 方法预先处理的。

.assess-value 通过在 ($!value ~^ $target) 两者之间进行字符 XOR 来计算值与目标值的匹配程度。对于每一个匹配的字符(或字符串化数字中的数字),XOR 将产生一个空字符('/0')。对于每一个不同的字符,产生的 XOR 字符将是其他的东西。所以我们可以通过计算 XOR 产生的前导 '\0' 字符的数量(.match( /^ \0+/ ).chars)来确定值的初始字符有多少与目标匹配。

然后,我们只需使用正则表达式将值串拆分成三个子串,使用 Terminal::ANSIColor 模块将前导匹配的字符染成蓝色,把尾部不匹配的字符染成红色。

一旦我们有了可用的角色,我们就可以建立几个适合我们实际数据的 Reportable 类。我们产生的每一个结果都需要根据 𝑒 的准确表示方式进行评估:

class Result
    does Reportable[:target<2.71828182845904523536028747135266>]
{}

而且我们还希望我们的报告在各种技术之间用空行来格式化,所以我们需要一个 Reportable 类,它的 .display 方法被重载,只打印空行:

class Spacer
    does Reportable
{
    method display (Int) { say "" }
}

最后我们建立 assess 子程序本身:

constant TIMEOUT = 3;  # seconds;

# Evaluate a block over a range of inputs, and report...
sub assess (Unary $block --> Any) {
    # Introspect the block's single parameter...
    my Parameter $param    = $block.signature.params[0];
    my Bool      $is-topic = $param.name eq '$_';

    # Extract and normalize the range of test values...
    my Any @values = do given try $param.default.() {
        when !.defined && $is-topic { Empty           }
        when !.defined              {    1, *×2 ... ∞ }
        when .?infinite             { .min, *×2 ... ∞ }
        default                     { .Seq            }
    }

    # Introspect the test description (from doc comments)...
    my Str $label = "[$block.line()] {$block.WHY // ''}".trim;

    # New paragraph in the report...
    Spacer.new;

    # Run all tests...
    for @values Z $label,"",* -> ($value, $label) {
        Result.new:
            desc  => "$label ($param.name() = $value)",
            value => timed { $block($value) },
            check => { last if .timing > TIMEOUT }
    }
    if !@values {
        Result.new:
            desc  => $label,
            value => timed { $block(Empty) }
    }
}

这个子程序只接受一个参数:一个 Code 对象,它本身接受一个参数(我们通过给参数类型 Unary 来执行)。

我们立即内省那一个参数($block.signature.params[0]),它(自然)是一个 Parameter 对象,并把它存储在一个合适的类型变量(my Parameter $param)中。我们还需要知道这个参数是否是隐式的主题变量(也就是 $_),所以我们也要测试一下。

一旦我们有了参数,我们需要确定调用者是否给了它一个默认值……这将代表我们要评估传递给评估块中的代码的一组值。换句话说,如果用户写的是:

assess -> \N=1..100 { some-function-of(N) }

…然后我们需要提取默认值 (1..100),这样我们就可以遍历这些值,并将每个值依次传递到块中。

我们可以通过调用相应的内省方法来提取参数的默认值(如果有的话!): $param.default,它将返回另一个 Code 对象,当调用时产生默认值。因此,要想得到实际的默认值,我们需要调用 .default,然后调用 .default 返回的代码。即: $param.default.()

当然,参数可能没有默认值,在这种情况下,.default 会返回一个未定义的值。在该未定义值上尝试第二次调用将是致命的,所以我们在 try 里面进行尝试,这将异常转换为另一个未定义值。

然后我们对提取的默认值进行测试,以确定它的含义:

when !.defined && $is-topic { Empty           }
when !.defined              {    1, *×2 ... ∞ }
when .?infinite             { .min, *×2 ... ∞ }
default                     { .Seq            }

如果它是未定义的,那么就没有指定默认值,所以我们要么使用一个空列表作为我们的测试值(如果参数只是隐式的主题,在这种情况下,它是一个无参数的一次性试验),否则我们就使用无休止的加倍序列 1,2,4 ... ∞。使用这个序列而不是 1..∞ 可以让我们在每一个数字数量级上都有合理的覆盖率,而不需要枯燥地尝试每一个可能的值。

如果默认是一个范围到无穷大(when .?infinite),那么我们将其调整为类似的翻倍值到无穷大的序列,从同一个下限(.min)开始。而如果是其他的值,我们就按原样使用,只是把它转换为一个序列(.Seq)。

一旦我们有了测试值,我们就需要对整个测试进行描述。我们本可以简单地为 assess 添加第二个参数,但 Raku 强大的内省功能提供了一个更有趣的选择……

我们需要传达三个信息:调用 assess 的行号,被评估的试验的描述,以及每个试验的试验值。但是,为了避免不伦不类的冗余,只需要在第一个试验中标明行号和描述,后续试验只需要显示下一个试验值。

我们可以通过内省代码块来获得行号: $block.line()。但是我们从哪里可以得到描述呢?

嗯……为什么不直接从注释中读取呢?

在 Raku 中,任何以竖条开头的注释(即 #| Your comment here)都是一种特殊形式的文档,被称为声明符块。当您指定这样的注释时,它的内容会自动附加到注释后面的第一个声明中。这可能是一个用 my 声明的变量,一个用 submulti 声明的子例程,或者(在这种情况下)一个在两个括号之间声明的匿名代码块。

换句话说,当我们这样写的时候:

#| Bernoulli's limit
assess -> \x { (1 + x⁻¹)**(x) }

…字符串 “Bernoulli’s limit” 会自动附加到传递给 assess 的块中。我们也可以将文档注释放在代码块之后,使用 #= 作为注释引入符,而不是 #|:

assess -> \x { (1+x⁻¹)**x }   #= Bernoulli's limit

无论哪种方式,Raku 的高级内省设施意味着我们可以在程序执行过程中检索该文档,只需调用块的 .WHY 方法(因为注释应该是告诉你"为什么"(WHY),而不是 “如何”(HOW))。

所以我们可以通过询问代码块的行号和文档来生成我们的描述,如果没有提供文档注释,就用一个空字符串来解决("[$block.line()] {$block.WHY // ''}")。

现在我们已经准备好运行测试并生成报告了。我们首先在报告中插入一个空行……通过声明一个 Spacer 对象。然后,我们在值列表中迭代,以及它们相关的标签,用 “zip” 操作符(Z)将两者交错在一起。标签是我们之前建立的初始标签,然后是一个空字符串,然后是一个 “whatever”(*),它告诉 zip 操作符重复使用前面的空字符串,次数不限,以匹配 @values 的剩余元素。这样一来,我们在试验的第一行就能得到完整的标签,但之后就不会重复了。

zip 产生了一系列的双元素数组:一个值,一个标签。然后我们在这些数组中迭代,为每个测试建立一个合适的 Result 对象:

Result.new:
    desc  => "$label ($param.name() = $value)",
    value => timed { $block($value) },
    check => { last if .timing > TIMEOUT }

报告的描述是标签,后面是块的参数名($param.name())和本次迭代传递给它的当前试用值($value)。报表的值只是调用当前传递给它的试用值的块的定时结果(timed { $block($value) })。

但是我们也希望在测试速度太慢的时候停止测试,所以我们把 Result 构造函数也传给一个 check 参数,一旦结果被添加到报告列表中,就会导致其 tweak 子方法执行 check 块。然后我们安排 check 块在新对象的定时超过我们选择的超时时长时终止周围的 for 循环(即 last)。对 .timing 的调用将在 Result 对象上被调用,而 Result 对象(本身没有 timing 方法)将把调用委托给它的 $.value 属性,我们指定它来处理对象本身无法处理的"任何"方法。

唯一要做的就是覆盖用户完全不提供测试值的边缘情况(即没有 @values 可以迭代的时候)。这可能发生在传递给 assert 的显式空默认值列表时,或者是传递进来的块根本没有显式声明参数时(因此默认为隐式主题参数: $_)。无论哪种情况,我们都需要准确地调用该块一次,不需要任何参数。因此,如果测试值数组没有元素(if @values.elems == 0),我们就会实例化一个单一的 Result 对象,向它传递完整的标签,以及用(字面上)空参数列表调用块的定时结果:

if !@values {
    Result.new:
        desc  => $label,
        value => timed { $block(Empty) }
}

然后我们就完成了。我们现在有了一种内省和紧凑的方法,可以在一个块上进行多次试验,跨越一系列合适的值,无论是显式的还是推断的,每次试验都能自动计时和超时。

所以,让我们回到计算 𝑒 的值……

𝑒 is for eversion

常数 𝑒 与雅各布·伯努利有另一种完全无关的联系。在他 1712 年的遗作 Ars Conjectandi 中,伯努利探讨了二项式试验的数学:随机试验的结果是严格的二进制:成功或失败,真或假,是或不是,头或尾。

二项式理论的一个结果与极度倒霉的概率有关。如果我们进行一个随机的二项式试验,成功的概率是 ¹⁄ₖ,那么失败的概率一定是 1 - ¹⁄ₖ。也就是说,如果我们重复进行同样的随机试验 k 次,那么每次失败的概率就是 (1 - ¹⁄ₖ)ᵏ。随着 k 的增大,这个概率从 0(k=1),到 0.25(k=2),到 0.296296296…(k=3),再到 0.296296…(对于 k=3),到 0.31640625(对于 k=4),逐渐收敛于一个渐进值 0.36787944117144233…(对于 k=∞)。而 0.36787944117144233 的值正好是 1∕𝑒。

也就是说,当 k 趋向于 ∞ 时,(1 - ¹⁄ₖ)⁻ᵏ 一定趋向于 𝑒。所以我们可以试一试:

#| Bernoulli trials
assess -> \k=2..∞  { (1 - k⁻¹) ** -k }

…它告诉我们:

    [20] Bernoulli trials (k = 2): 4                     [0.000𝑠]
                          (k = 4): 3.160494              [0.000𝑠]
                          (k = 8): 2.910285368           [0.000𝑠]
                         (k = 16): 2.808403965576447     [0.000𝑠]
                         (k = 32): 2.762009089976452     [0.000𝑠]
                         (k = 64): 2.7398271814872146    [0.000𝑠]
                        (k = 128): 2.7289767306942787    [0.000𝑠]
                        (k = 256): 2.7236100544173554    [0.000𝑠]
                        (k = 512): 2.720941162086445     [0.000𝑠]
                       (k = 1024): 2.7196103037796897    [0.000𝑠]
                       (k = 2048): 2.7189457686628256    [0.001𝑠]
                       (k = 4096): 2.7186137242488035    [0.002𝑠]
                       (k = 8192): 2.7184477577823865    [0.011𝑠]
                      (k = 16384): 2.718364788478643     [0.059𝑠]
                      (k = 32768): 2.7183233073084274    [0.265𝑠]
                      (k = 65536): 2.718302567593645     [1.232𝑠]
                     (k = 131072): 2.7182921979538235    [6.135𝑠]

遗憾的是,这比伯努利最初的高利贷计划收敛得再快。

𝑒 is for exclamation

尽管伯努利试验的表现令人失望,但伯努利试验的方法突出了一个有用的想法。常数 𝑒 出现在许多其他数学方程中,所有这些方程我们都可以重新排列,产生一个公式,开始:𝑒 = …

例如,在1733年,亚伯拉罕-德-莫伊夫尔首先发表了一个阶乘函数的渐近式,(仅仅几天后!)詹姆斯-斯特林(James Stirling)对这个渐近式进行了改进,现在这个渐近式以他的名字命名。这个近似值是:n! ≅ √2πn × (n∕𝑒)ⁿ。随着 n 越大,近似值越精确。

将该方程重新排列,我们得到: 𝑒 ≅ ⁿ⁄ₙ√n!

对于这个方程,我们需要一个 n 次方根运算符和一个阶乘运算符。这两个都是标准 Raku 中所缺少的。这两个都是很容易添加到它:

# n-th root of x...
sub infix:<√> (Int $n, Int $x --> Numeric) is tighter( &[**] ) {
    $x ** $n⁻¹
}

# Factorial of x...
sub postfix:<!> (Int $x --> Int)  {
    [×] 1..$x
}

新的中缀 √ 运算符被指定为比现有的 ** 指数运算符(is tighter( &[**] ))更高的优先级。它只是将 √ 之后的数值提高到 √ 之前数值的倒数。

新的后缀 ! 运算符将所有小于或等于它的单个参数(1..$x)的数字乘在一起,使用乘法([×])化简

有了这两个新的运算符,我们现在可以写:

#| de Moivre/Stirling approximation
assess -> \n { n / n√n! }

遗憾的是,结果却不尽如人意:

  [30] de Moivre/Stirling approximation (n = 1): 1                     [0.000𝑠]
                                        (n = 2): 1.414213562373095     [0.000𝑠]
                                        (n = 4): 1.80720400721969      [0.000𝑠]
                                        (n = 8): 2.1252005594420327    [0.000𝑠]
                                       (n = 16): 2.352779665665871     [0.000𝑠]
                                       (n = 32): 2.501898064970579     [0.000𝑠]
                                       (n = 64): 2.5938155227045367    [0.001𝑠]
                                      (n = 128): 2.6481531254600723    [0.001𝑠]
                                      (n = 256): 0                     [0.001𝑠]
                                      (n = 512): 0                     [0.006𝑠]
                                     (n = 1024): 0                     [0.011𝑠]

这种方法收敛在 𝑒 上的速度甚至比之前的还要慢,这也是 Raku 第一次在数值计算上真正让我们失望。这些零是说明编译器已经无法取 256 的阶乘的 256 次方根。或者更高数字的深根。

它计算的是阶乘本身(也就是 8578177753428426541190822716812326251577815202 79485619859655650377269452553147589377440291360451408450375885342336584306 15719683469369647532228928849742602567963733256336878644267520762679456018 79688679715211433077020775266464514647091873261008328763257028189807736717 81454170250523018608495319068138257481070252817559459476987034665712738139 28620523475680821886070120361108315209350194743710910172696826286160626366) 但是取这个巨大数字的 256 方根(将其升到 ¹⁄₂₅₆ 的幂)却失败了,它错误地产生了值 ∞,于是 n 除以那个错误的结果产生了零。失败点似乎在 170! 或 10³⁰⁸ 左右,这看起来很可疑,像是内部 64 位浮点表示极限。

当然,我们可以通过改变计算n 次方根的方式来解决这个限制。例如,一个数字的 n 次方根也是由以下公式给出的: 10ˡᵒᵍ₁₀ˣ⁄ₙ

所以我们可以试试:

# n-th root of x...
sub infix:<√> (Int $n, Int $x --> Numeric) is tighter(&[**]) {
    10 ** (log10($x) / $n)
}

…但我们得到了完全相同的问题:内置的 log10 函数在 10³⁰⁸ 左右就失效了,这意味着我们不能将其应用于大于 170 的数字的阶乘。

但是当 X 是一个整数时,有一个出人意料的 log₁₀X 近似方法,它完全没有 10³⁰⁸ 的限制:我们只需计算 X 中的数字数,然后减去 ½。在最坏的情况下,这个结果的偏差不超过 ±0.5,但这意味着它在足够大的数值范围内的平均误差…是零。

用这种近似方法计算 n 次方根:

# n-th root of x...
sub infix:<√> (Int $n, Int $x --> Numeric) is tighter(&[**]) {
    10 ** (($x.chars - ½) / $n)
}

…我们现在可以把对德莫弗尔/斯特林近似的评估扩大到 n=170 以外。在这种情况下,我们发现:

  [30] de Moivre/Stirling approximation (n = 1): 0.31622776601683794   [0.000𝑠]
                                        (n = 2): 1.1246826503806981    [0.000𝑠]
                                        (n = 4): 1.6867860137143293    [0.000𝑠]
                                        (n = 8): 2.190735707411489     [0.000𝑠]
                                       (n = 16): 2.2928201123791405    [0.000𝑠]
                                       (n = 32): 2.4875680967640825    [0.000𝑠]
                                       (n = 64): 2.5570691577279274    [0.001𝑠]
                                      (n = 128): 2.6522607579574027    [0.001𝑠]
                                      (n = 256): 2.6898269484112647    [0.005𝑠]
                                      (n = 512): 2.69742667825642      [0.006𝑠]
                                     (n = 1024): 2.7080909004243128    [0.012𝑠]
                                     (n = 2048): 2.711166091674111     [0.040𝑠]
                                     (n = 4096): 2.7150077947508775    [0.159𝑠]
                                     (n = 8192): 2.716181527473568     [0.678𝑠]
                                    (n = 16384): 2.7171648275814224    [3.049𝑠]

令我们厌恶的是,这种方法的收敛性显然在 n 值越高的情况下没有加快。

事实上,即使 n 超过100万,结果仍然只能精确到小数点后三位;甚至比经典的回文分数近似法更差。

#| Palindromic fraction
assess {   878
         / 323
       }

…这给了我们:

[40] Palindromic fraction: 2.718266              [0.001𝑠]

唉,搜索还在继续。

𝑒 is for estimation

我们再来看看概率论。想象一系列均匀分布在 0..1 范围内的随机数,例如:

0.162532, 0.682623, 0.4052625, 0.42167261, 0.99765261, ...

如果我们开始将这些数字相加,那么我们需要将序列的多少项相加才会使总数大于1呢?在上面的例子中,我们需要将前三个值相加才能超过1。但在其他随机序列中,我们只需要加两个项(如 0.76217621,0.55326178,…),偶尔我们还需要加不少项(如 0.1282827,0.00262671,0.39838722,0.386272617,0.77282873,…)。然而,在大量的试验中,和值超过1所需的平均随机值数是(你猜对了): 𝑒。

因此,如果我们有一个范围为 0..1 的均匀随机值的来源,我们可以通过反复添加足够的随机值来获得一个近似值 𝑒 以超过1,然后平均每次所需的值的数量,在多个试验中。这在 Raku 看来是这样的:

#| Stochastic approximation
assess -> \trials=(10², 10³, 10⁴, 10⁵, 10⁶) {
    sum( (1..trials).map:
            { ([\+] rand xx ∞).first(* > 1):k + 1 }
    )
    / trials
}

即:我们进行指定数量的试验(1..trials),对于每一次试验(.map:),我们都会产生无限个均匀的随机值(rand xx ∞),然后我们将其转换为一个渐进式部分和的列表([\+])。然后我们寻找这些部分和中第一个超过1的(.first(* > 1)),找出它在列表中的索引(:k)并加1(因为索引从0开始,而计数从1开始)。

结果是我们每次试验中需要多少随机值超过1的计数列表。由于 𝑒 是这些计数的平均值,我们将它们相加,然后除以试验次数。然后发现。

   [50] Stochastic approximation (trials = 100): 2.62                 [0.020𝑠]
                                (trials = 1000): 2.7                  [0.045𝑠]
                               (trials = 10000): 2.7179               [0.357𝑠]
                              (trials = 100000): 2.7256               [3.651𝑠]
                             (trials = 1000000): 2.718698            [43.185𝑠]

这很好…对于随机猜测。如果我们让它运行更多的试验,我们最终会得到合理的准确性。但它不可靠:有时随着试验次数的增加而失去准确性。而且它的速度很慢:在一百万次试验之后,只有三位数的准确率…而且计算时间为40秒。为了得到一个有用的正确数字,我们需要数十亿次,可能是数万亿次的试验…这就需要数万或数千小时的计算。

所以我们继续…

𝑒 is for efficiency

或者说,我们回到了过去。回到1669年 数学天才中的艾萨克-牛顿身边 艾萨克-牛顿 在一份题为De analysi per aequationes numero terminorum infinitas的手稿中,牛顿提出了一种解无限级数方程的一般方法,这种方法意味着一种高效的方法来确定 𝑒。

具体地说,就是 𝑒 = Σ ¹⁄ₖ!,其中 k 从 0 到 ∞。

那我们就试试吧:

#| Newton's series
assess -> \k=0..∞ { sum (0..k)»!»⁻¹ }

在这里,我们通过从无限项系列(\k=0..∞)中取长度为 k 的越来越长的子集来计算 𝑒。我们取所选项的指数((0..k)),并对每一个项(»)取其阶乘(!)。然后对于每个阶乘(»)我们取其倒数()。结果是一个连续项 ¹⁄ₖ! 的列表,最后我们把这些项加在一起(sum),得到 𝑒。

    [60] Newton's series (k = 0): 1                       [0.001𝑠]
                         (k = 1): 2                       [0.000𝑠]
                         (k = 2): 2.5                     [0.000𝑠]
                         (k = 3): 2.666667                [0.000𝑠]
                         (k = 4): 2.708333                [0.001𝑠]
                         (k = 5): 2.716667                [0.001𝑠]
                         (k = 6): 2.718056                [0.002𝑠]
                         (k = 7): 2.718254                [0.001𝑠]
                         (k = 8): 2.718279                [0.005𝑠]
                         (k = 9): 2.718282                [0.001𝑠]
                        (k = 10): 2.718281801             [0.001𝑠]
                        (k = 11): 2.718281826             [0.001𝑠]
                        (k = 12): 2.7182818283            [0.002𝑠]
                        (k = 13): 2.718281828447          [0.002𝑠]
                        (k = 14): 2.7182818284582         [0.002𝑠]
                        (k = 15): 2.71828182845899        [0.002𝑠]
                        (k = 16): 2.7182818284590423      [0.002𝑠]
                        (k = 17): 2.718281828459045       [0.002𝑠]
                        (k = 18): 2.718281828459045227    [0.002𝑠]
                        (k = 19): 2.7182818284590452      [0.003𝑠]
                        (k = 20): 2.71828182845904523534  [0.007𝑠]
                        ⋮

终于…取得了一些真正的进展。在只对无限的 ¹⁄ₖ! 数列的20个项相加后,我们有了 19 个正确的小数点,(终于!)有了一个相当准确的值 𝑒。

我们还可以做得更好。𝑒 可以分解成许多其他的无限和。比如说:哈兰兄弟和约翰-诺克斯发现 𝑒 = Σ ⁽²ᵏ⁺¹⁾⁄₍₂ₖ₎! 我们可以用以下方法评估:

#| Brothers/Knox series
assess -> \k=0..∞ {
    sum (2 «×« (0..k) »+» 1) »/« (2 «×« (0..k))»!
}

通过将每个算术运算符都变成超运算符,我们可以在一个向量表达式中计算出整个 0..k 项系列,然后将它们相加得到 𝑒。

哈兰兄弟公式的代码可能比牛顿的原始公式要丑陋一些,但它的收敛速度是牛顿公式的两倍,仅从系列的前十个项中就能得到 19个 正确的小数点。

    [70] Brothers/Knox series (k = 0): 1                       [0.003𝑠]
                              (k = 1): 2.5                     [0.000𝑠]
                              (k = 2): 2.708333                [0.001𝑠]
                              (k = 3): 2.718056                [0.001𝑠]
                              (k = 4): 2.718279                [0.001𝑠]
                              (k = 5): 2.718281801             [0.001𝑠]
                              (k = 6): 2.7182818283            [0.001𝑠]
                              (k = 7): 2.7182818284582         [0.001𝑠]
                              (k = 8): 2.7182818284590423      [0.001𝑠]
                              (k = 9): 2.718281828459045227    [0.001𝑠]
                             (k = 10): 2.71828182845904523534  [0.001𝑠]
                             ⋮

𝑒 is for epsilon

我们在用这两个牛顿数列精确计算 𝑒 的问题上取得了扎实的进展,但那个丑陋的(可能是可怕的)超运算符计算代码还是挺让人生气的。而且太容易出错了。

鉴于这些高效方法的工作方式都是一样的-对无限系列项的(初始子集)求和–如果我们有一个函数来帮我们做这件事,也许会更好。而且,如果函数能够自行计算出它究竟需要包含多少系列的初始子集才能产生一个准确的答案….而不是要求我们手动梳理多次试验的结果来发现这一点,那肯定会更好。

而且,就像 Raku 中经常出现的那样,构建我们所需要的东西是出奇的容易:

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

我们称这个子程序为 Σ,因为这是这个函数的常用数学符号(好吧,是的,只是因为我们可以)。我们指定它接受一个参数块,并返回一个数值(Unary $block --> Numeric),返回值代表块指定的系列的"足够精确"的求和。

为了达到这个目的,它首先建立无限的项索引列表(0..∞),并将每个项索引依次传递给块(.map($block))。然后,它对所产生的项建立连续的部分和(.produce(&[+]))。也就是说,如果块是 { 1 / k! },那么 .map 创建的列表将是: (1, 1, 1∕2, 1∕6, 1∕24, 1∕120, ...) 然后 .produce(&[+]) 方法的调用将逐步汇总这些项,产生列表: 1, 1+1, 1+1+1∕2, 1+1+1∕2+1∕6, 1+1+1∕2+1∕6+1∕24, 1+1+1∕2+1∕6+1∕24+1∕120, ...)也就是:(1, 2, 2.5, 2.666, 2.708, 2.717, ...)。

我们之所以这样设置计算,是因为我们需要能够计算出什么时候停止添加项,也就是当 .produce 列表的连续元素收敛到一个一致的值时。换句话说,就是当两个连续的值相差不大的时候。

而且,方便的是,Raku 有一个运算符可以准确地告诉我们:近似相等运算符(或其 ASCII 等价物:=~=)。

所以我们可以写一个子程序,接受一个列表,并从中返回第一个 “收敛"的值,就像这样:

sub converge(Iterable $list --> Numeric) {
    $list.rotor(2=>-1).first({ .head ≅ .tail }).tail
}

.rotor 方法从其列表中提取 N 个元素的子列表。在这里,我们告诉它一次提取两个元素,但也让它们重叠,在每次提取之间退后一位(-1)。结果是,从一个列表中,如(1, 2, 2.5, 2.666, 2.708, 2.717 …),我们得到一个由每一对相邻值组成的列表:((1, 2), (2, 2.5), (2.5, 2.667), (2.667, 2.708), (2.708, 2.717) …)。

然后,我们简单地在列表中步步为营,寻找头元素和尾元素近似相等的第一个子列表(.first({ .head ≅ .tail }))。这样我们就得到了一个单一的子列表,我们只从中提取比较准确的尾元素(.tail)。

然后我们就可以直接将这个函数应用到 Σ 中产生的部分和的列表中:

(0..∞).map($block).produce(&[+]).&converge

…这只是将整个列表从左到右传递给 converge 的一种方便的方式。当然,我们也可以把它写成一个普通的 name-arglist 风格的子例程调用:

converge  (0..∞).map($block).produce(&[+])

…如果这能让我们更舒服的话。

关于这一切的超酷之处–也是它工作的唯一原因–是 Σ 和 converge 这两个方法调用链中的每一个组件都被懒惰地计算了。这意味着 .map 并不会在任何一个(0..∞)值上实际执行块,.produce 也不会将其中任何一个值逐步加在一起,.rotor 也不会提取其中任何一对重叠的值,.first 也不会搜索其中任何一个值…除非最终的结果要求它们这样做。.map 只会根据需要多次调用块,以获得足够的值让 .produce 添加,获得足够的值让 .rotor 提取,获得足够的值让 .first 找到第一个大致相等的对。

我真的很喜欢 Raku 的这一点:它不仅能让你写出简洁、表现力强、优雅的代码,你简洁、表现力强、优雅的代码自然也是高效的。

同时,我们现在有了一个更好的(即更简洁、更有表现力、更优雅)的方式来生成一个高度准确的 𝑒 值:

#| Newton's series
assess { Σ -> \k { 1 / k! } }

#| Brothers/Knox series
assess { Σ -> \k { (2×k + 1) / (2×k)! } }

…这给了我们直接的答案,我们正在寻找:

         [80] Newton's series: 2.718281828459045227    [0.007𝑠]

    [83] Brothers/Knox series: 2.71828182845904523534  [0.002𝑠]

𝑒 is for epigram

这将是故事的结束…除了我们仍然忽略了一半的可能性。到目前为止,我们所尝试的每一种技术都是数学性质的。但 Raku 不只是为代数学家、统计学家或概率理论家准备的。它也适用于语言学家、作家和诗人以及所有其他自然语言的爱好者。

那么,我们如何利用自然语言计算出一个准确的 𝑒 值呢?

好吧,事实证明,代数学家、统计学家和概率理论家已经做了几个世纪。因为一旦你花了几个小时(或几天,或几周)手动计算了 𝑒 的前八位数字,你就再也不想这样做了! 所以你想出了一个记忆法:一句话,帮助你记住这些来之不易的数字。

大多数这类记忆法的工作原理是将每个数字编码在连续单词的长度中。例如,要记住常数 1.23456789,你可以把它编码为:“我是当地唯一一个公开的澳大利亚算盘手”(I am the only local Aussie abacist publicly available)。然后你只需计算每个单词的字母来提取数字。

向往常一样,这在 Raku 中实现起来很琐碎:

sub mnemonic(Str $text --> Str) {
    with $text.words».trans(/<punct>+/ => '')».chars {
        return "{.head}.{.tail(*-1).join}"
    }
}

我们首先从文本中提取每一个字($text.words),对于其中的每一个字(»),我们通过将它们翻译成空字符串(.trans(/<punct>+/ => ''))来删除任何标点符号,最后计算每一个字中剩余的字符数(».chars)。然后,我们从该字长列表中取出第一个元素(.head),加上一个点,然后将其余的 N-1 个字符数的连接(.tail(*-1).join)追加,并将其作为结果数的字符串表示返回。

然后我们对其进行测试:

#| Mnemonic test
assess {
    mnemonic "I am the only local Aussie abacist
              publicly available"
}

…它打印出:

    [90] Mnemonic test: 1.23456789            [0.002𝑠]

…这显然是一个非常糟糕的近似于 𝑒。

但三个世纪以来,人们一直在编造更好的。其中使用最广泛的是略带自嘲的:

#| Modern mnemonic
assess {
    mnemonic "It enables a numskull to remember
              a sequence of numerals"
}

…这样我们就有了9个正确的小数位:

    [100] Modern mnemonic: 2.718281828           [0.003𝑠]

如果我们需要更准确,我们可以直接组成一个更长的句子。例如:

#| Titular mnemonic
assess {
    mnemonic "To compute a constant of calculus:
              (A treatise on multiple ways)"
}

…它产生一个额外的数字:

    [110] Titular mnemonic: 2.7182818284          [0.002𝑠]

或者,为了比目前的任何一种数学方法更精确,我们可以使用 Zeev Barel 的自述:

#| Extended mnemonic
assess {
    mnemonic "We present a mnemonic to memorize a constant
              so exciting that Euler exclaimed: '!'
              when first it was found. Yes, loudly: '!'.
              My students perhaps will compute e, via power
              or Taylor series, an easy summation formula.
              Obvious, clear, elegant!"
}

…这给我们提供了远比我们实际需要的更多的准确性:

 [120] Extended mnemonic: 2.718281828459045235360287471352662497757    [0.010s]

𝑒 is for eclectic

甚至这还不是故事的结尾。就像 π(它的主要对手是世界上最可怕的数学常数)一样,𝑒 对数学上的奇思妙想施加了一种奇怪的魅力。

例如,1988年,达里奥-卡斯特拉诺斯公布了以下公式:𝑒=(π⁴ + π⁵)¹⁄₆。其中,翻译成 Raku 就是:

#| Castellano's coincidence
assess { (π⁴ + π⁵) ** ⅙ }

…并给出了小数点后七位的精度:

[130] Castellano's coincidence: 2.718281808611915     [0.001𝑠]

而我们这些能记住数字1到9的人(或者认识澳大利亚算盘高手)可以用理查德-萨贝的公式: 𝑒 = (1 + 2-3×(4+5)).6×.7+89:

#| Sabey's digits
assess { (1+2**(-3×(4+5)))**(.6×.7+8⁹) }

…(恰如其分地)我们总共有9个正确的数字:

          [140] Sabey's digits: 2.718281826838823     [0.000𝑠]

但与真正的探索高地 𝑒 相比,这些好奇心如同无物。

比如,马克西米利安-皮斯科罗夫斯基发现,如果你恰好有8个空闲的9,就可以计算出 𝑒=(9∕9+9-99)999,精确到小数点后 3.69 亿位多一点。

我们可以很容易地把它翻译成 Raku。

#| Piskorowski's eight 9s
assess { (9/9 + 9**-9**9) ** 9**9**9 }

但是,可惜的是,我们的评估失败了(最终),产生了令人失望的结果:

    [150] Piskorowski's eight 9s: 1             [332242.506𝑠]

…因为 9⁻⁹ 的值是如此的微不足道,以至于在 Raku 中把它加到 9∕9 上只产生1,然后 1999 还是1。

Piskorowski 的不可估量的 ⅓ 亿个小数位的准确度似乎是这个探索的(il)逻辑终点。但直到前述理查德-萨比破坏了大家的游戏,他将上述的泛数字公式重新表述为 𝑒 = (1 + 9-47×6)3285

这精确到惊人的 18457734525360901453873570 个小数位(即18.4八十亿位)…但悲剧的是,当试图计算初始的 9-442时,它立刻溢出了 Raku 的数字表示法;Raku 无法准确表示 ¹⁄₉ 的 19342813113834066795298816 次方。

同时,如果你想进一步探索这些和其他𝑒相关的奇闻异事,可以去 Erich Friedman 的令人着迷的数学魔法网站看看。

𝑒 is for exquisite

在我们追求 𝑒 的过程中,还有一个更高的数学高峰需要我们去攀登。数学优雅的巅峰,算术的巅峰,所有数学中最美的一个方程: 欧拉方程

1748年,莱昂哈德-欧拉发表了他的代表作: Introductio in analysin infinitorum,其中他首次提到了一般的欧拉公式: 𝑒ⁱˣ = cos x + i sin x. 虽然在论文中没有明确提到,这个公式暗示了一个显著的特殊情况(在 x=π 时): 𝑒iπ = -1 + 0i。

这种特殊情况下的平等现在更多的是写成: 𝑒iπ+1=0 这样就把数学中最重要的五个常数统一起来了。

除了这个格外可爱的数学基础的统一之外,对于我们的目的来说,这里重要的一点是,五个成分之一是 𝑒。而且,如果我们重新排列公式来分离这个常数,我们得到: 𝑒 = (-1)1∕πi 我们可以很容易地在 Raku 中评估。

#| From Euler's Identity
assess { (-1) ** (π×i)⁻¹ }

不幸的是,Raku 还不够聪明,无法从虚指数中推断出它需要使用复数版的 ** 运算符(嗯,是的,当然复数已经内置在 Raku 中了)。所以我们回来了:

    [160] From Euler's Identity: NaN                   [0.000𝑠]

…因为当给定一个复数指数时,非复数指数就会失效,产生 NaN+NaN/i 的值。

但是如果我们从-1的适当的二维版本(即-1+0i)开始,就像这样:

#| From Euler's Identity
assess { (-1+0i) ** (π×i)⁻¹ }

…我们得到正确的复数运算,并从中得到高度准确的答案:

    [160] From Euler's Identity: 2.718281828459045    [0.000𝑠]

𝑒 is for effortless

我们现在已经绕了一圈:从欧拉公式的特殊情况中找到欧拉数,这很合适,这就是欧拉等式。可以说,欧拉的欧罗伯罗斯。

所有的努力… 最终多少回到了我们开始的地方。这是非常合适的,因为我们开始的地方已经有了答案。首先,以标准 exp 函数的形式,它返回的值是 𝑒ˣ。

所以我们可以尝试一下:

#| Built-in exp() function
assess -> \x=1 { exp(x) }

…这将给我们:

    [170] Built-in exp() function (x = 1): 2.718281828459045    [0.000𝑠]

但即使是这样,我们在 Raku 中也需要付出更多的努力。因为(当然)常量本身也是直接内置在语言中的:

#| Built-in e constant
assess { e }

…并产生同样准确的结果:

    [180] Built-in e constant: 2.718281828459045    [0.000𝑠]

所以原任务的最优解只有五个字符长。

say e

…此后,将被称为:“欧拉的单行程序”。

by Damian Conway