⍟⍵

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 !

⍟⍵

Postby Roger|Dyalog on Sat Mar 06, 2021 1:05 am

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 ○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

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

Return to APL Chat

Who is online

Users browsing this forum: No registered users and 1 guest