⍟⍵

[Roger|Dyalog]

The recent APL Chat Forum post * versus ⍟ Performance [Hui 2021] briefly mentioned how one might compute ⍟⍵. Computing ⍟⍵ on 64-bit floating-point numbers turns out to be easier than described.

64-Bit Floating-Point Numbers

We aim to provide a model of ⍟⍵ where ⍵ is a 64-bit floating-point number ("double"). Such numbers are specified by the IEEE 754-1985 standard [Wikipedia 2021]:

• 00-00 ___ a sign bit (0 means positive, 1 means negative)
• 1-110 ___ 11 bits of the base-2 exponent + 1023
• 12-63 ___ 52 bits of the mantissa, with an assumed leading bit of 1 when the number is nonzero

That is, if s is the sign bit, e the exponent bits, and x the mantissa bits, then the value of a nonzero number is

      (¯1*s) × (2*¯1023+2⊥e) × 1+(2⊥x)÷2*52

Dyadic ⎕dr provides a way to get at the bits of a floating-point number:

      ⍝
drep←{                                                         
  ⊃83 ⎕dr 1 256:,⊖8 8⍴11 ⎕dr ⍵,0.1    ⍝ reverse byte-order
  64⍴11 ⎕dr ⍵,0.1                     ⍝ normal  byte-order
}

      ⊢ b←drep ○1
0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 1 1 1 1 1 0 1 1 0 1 0 ...
      ⍴b
64
      8 8⍴b
0 1 0 0 0 0 0 0
0 0 0 0 1 0 0 1
0 0 1 0 0 0 0 1
1 1 1 1 1 0 1 1
0 1 0 1 0 1 0 0
0 1 0 0 0 1 0 0
0 0 1 0 1 1 0 1
0 0 0 1 1 0 0 0

      sb ← 1↑b      ⍝ sign bit
      eb ← 11↑1↓b   ⍝ exponent bits
      xb ← 12↓b     ⍝ mantissa bits

      ⊢ pie ← (¯1*sb) × (2*¯1023+2⊥eb) × 1+(2⊥xb)÷2*52
3.14159
      pie - ○1
0

Note that in the representation, the mantissa x for a non-zero number satisfies (1≤x)∧x<2.

Computing ⍟⍵

To compute ⍟⍵, we use equation 4.1.27 of Abramowitz and Stegun's Handbook of Mathematical Functions [Abramowitz and Stegun 1964]:

      ⍟⍵ ←→ 2 × +/ i ÷⍨ t*i ⊣ t←(⍵-1)÷(⍵+1) ⊣ i←1+2×⍳n

The following properties of ⍟⍵ are used:

      ⍟ a×b  ←→  (⍟a) + (⍟b)
      ⍟ a*e  ←→  e × ⍟a

Therefore, ⍟ x × 2*e ←→ (⍟x) + e×⍟2. Since the mantissa x is between 1 and 2, the quotient (x-1)÷(x+1) in equation 4.1.27 is between 0 and ÷3; moreover, since the equation uses odd powers of (x-1)÷(x-1), the terms in equation 4.27 decline by a factor of (÷3)*2 (=0.111111) or less. That is, a fast converging series.

      ⍝
assert←{⍺←'assertion failure' ⋄ 0∊⍵:⍺ ⎕signal 8 ⋄ shy←0}

log←{    ⍝ Model of ⍟⍵.  See Abramowitz and Stegun, equation 4.1.27
  ⍺←0
  ⎕ct←0
  assert 0=≢⍴⍵:        ⍝ scalars only please
  assert 0<⍵:          ⍝ we only do positive ⍵
  d←drep ⍵
  e←¯1023+2⊥12↑d       ⍝ base-2 exponent
  x←1+(2⊥12↓d)÷2*52    ⍝ mantissa
  r←(x-1)÷x+1
  t←(r×r){⍺ ∇⍣((⊃⍵)≠(⊃⍵)+t÷2+≢⍵) ⊢⍵,t←⍺×⊃⌽⍵}r
  0=⍺:(2×+/⌽t÷1+2×⍳≢t)+e×0.6931471805599453    ⍝ e×⍟2
  1=⍺:e,x,⊂t÷1+2×⍳≢t
}

For example:

      log¨ 1e¯99 0.9 9 90 1e99
¯227.956 ¯0.105361 2.19722 4.49981 227.956

      ⍟    1e¯99 0.9 9 90 1e99
¯227.956 ¯0.105361 2.19722 4.49981 227.956

      (⍟ = log¨) 1e¯99 0.9 9 90 1e99
1 1 1 1 1

log has an optional left argument. 1 log ⍵ returns a 3-item result of the base-2 exponent e, the mantissa x, and the terms t of the power series.

      (e x t) ← 1 log 9

      e
3
      x
1.125
      x × 2*e
9

      t
0.0588235 0.0000678472 1.40859E¯7 3.48144E¯10 9.36952E¯13 
                                        2.65258E¯15 7.76642E¯18 2.32903E¯20
      (2×+/⌽t) + e×⍟2
2.19722
      log 9
2.19722

The sum of t is computed by +/⌽t instead of +/t because it is slightly more accurate to accumulate the smaller terms first. For historical reasons, in Dyalog APL +/⍵ traverses the vector from left to right rather than the specified right to left.

References

• Abramowitz, M., and I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964.
• Hui, R.K.W., * versus ⍟ Performance, Dyalog APL Chat Forum, 2021-03-02 08:47.
• Wikipedia, Double-Precision Floating-Point Format, 2021.