⍟⍵
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 !
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 !
1 post
• Page 1 of 1
⍟⍵
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]:
Dyadic ⎕dr provides a way to get at the bits of a floating-point number:
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]:
The following properties of ⍟⍵ are used:
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.
For example:
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.
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
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
(¯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.
- Roger|Dyalog
- Posts: 238
- Joined: Thu Jul 28, 2011 10:53 am
1 post
• Page 1 of 1
Who is online
Users browsing this forum: No registered users and 1 guest
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group