## ⍟⍵

APL-related discussions - a stream of APL consciousness.
Not sure where to start a discussion ? Here's the place to be
Forum rules
This forum is for discussing APL-related issues. If you think that the subject is off-topic, then the Chat forum is probably a better place for your thoughts !

### ⍟⍵

The recent APL Chat Forum post * versusPerformance [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 ○10 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 ...      ⍴b64      8 8⍴b0 1 0 0 0 0 0 00 0 0 0 1 0 0 10 0 1 0 0 0 0 11 1 1 1 1 0 1 10 1 0 1 0 1 0 00 1 0 0 0 1 0 00 0 1 0 1 1 0 10 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*523.14159      pie - ○10`

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 1e991 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      e3      x1.125      x × 2*e9      t0.0588235 0.0000678472 1.40859E¯7 3.48144E¯10 9.36952E¯13                                         2.65258E¯15 7.76642E¯18 2.32903E¯20      (2×+/⌽t) + e×⍟22.19722      log 92.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

Roger|Dyalog

Posts: 236
Joined: Thu Jul 28, 2011 10:53 am