aka.nice
aka.nice

Reputation: 9572

Which guard is necessary for log2

The question can be formulated like this:

can you find or construct a pair of positive integers i and p such that:

(i raisedTo: p) highBit / i log2 < p

or

((i raisedTo: p) highBit - 1) / i log2 > p

with log2 as defined in Squeak trunk, using some IEEE754 double precision approximation (see below).

We can exclude case when i isPowerOfTwo which is trivial since log2 is exact in this case.

The context is for determining if an Integer is an integral power of another integer.

Integer>>isPowerOf: anInteger
    "Return true if the receiver is an integral power of anInteger."

    | n iLog2 guard |
    self < 0 ifTrue: [^anInteger < 0 and: [self abs isOddPowerOf: anInteger abs]].
    anInteger < 0 ifTrue: [^self isEvenPowerOf: anInteger abs].
    anInteger < 2 ifTrue: [^self = anInteger].
    anInteger = 2 ifTrue: [^self isPowerOfTwo].
    anInteger isPowerOfTwo ifTrue: [^self isPowerOfTwo and: [self highBit - 1 \\ (anInteger highBit - 1) = 0]].
    self = 1 ifTrue: [^true].
    self < anInteger ifTrue: [^false].
    iLog2 := anInteger log2.
    guard := 2.0 * iLog2 ulp.
    n := (self highBit / (iLog2 - guard)) truncated.
    ^self highBit - 1 / (iLog2 + guard) < n and: [(anInteger raisedToInteger: n) = self]

The question occurs because I wonder if I can write the last lines without a guard, but using lesser or equal:

    iLog2 := anInteger log2.
    n := (self highBit / iLog2) truncated.
    ^self highBit - 1 / iLog2 <= n and: [(anInteger raisedToInteger: n) = self]

Assume anInteger**n=self
Then n*anInteger log2=self log2
But we know that self highBit - 1 < self log2 < self highBit
The inequalities are strict because we remove trivial case when anInteger isPowerOfTwo
That means that (self highBit - 1) / anInteger log2 < n < (self highBit / anInteger log2)
Since anInteger > 2, anInteger log2 > 1, and 1/anInteger log2 < 1,
there is at most one integer n in this interval.
We know that (self highBit - 1) isAnExactFloat, otherwise memory size would be > 2**54 which is not possible in immediate future.

With current implementation of log2 in Squeak, I can find some integer for which the error on log2 is a bit more than 0.5 ulp, and for which either:

(i raisedTo: p) highBit / i log2 = p

or:

((i raisedTo: p) highBit - 1) / i log2 = p

For example:

i := 1<<54-45.
i log2 - (i asArbitraryPrecisionFloatNumBits: 106) log2 / i log2 ulp.
    "=> 0.507..."
(1 to: 10) collect: [:p | (i**p) highBit / i log2].
    "=> #(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0)"

or:

i := 1<<55+89.
i log2 - (i asArbitraryPrecisionFloatNumBits: 106) log2 / i log2 ulp.
    "=> -0.501..."
(1 to: 10) collect: [:p | (i**p) highBit - 1 / i log2].
    "=> #(1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0)"

But so far, I failed to construct a counter example...

The strategy of log2 is to convert asFloat, unless overflow:

Number>>log2
    "Answer the base-2 log of the receiver."
    ^self asFloat log2

Float>>log2
    "Answer the base 2 logarithm of the receiver.
    Arrange to answer exact result in case of exact power of 2."
    |  s  |
     s := self significand.
    ^s > 1.3333333333333333
        ifTrue: [(0.5 * s) ln / Ln2 + (1 + self exponent)]
        ifFalse: [s ln / Ln2 + self exponent]

LargePositiveInteger>>log2
    "This function is defined because super log2 might overflow."
    | res h |
    res := super log2.
    res isFinite ifTrue: [^res].
    h := self highBit.
    ^h + (self / (1 << h)) asFloat log2

Fraction>>log2
    "This function is defined because super log2 might overflow."
    | res |
    self <= 0 ifTrue: [DomainError signal: 'log2 is only defined for x > 0'].
    "Test self < 1 before converting to float in order to avoid precision loss due to gradual underflow."
    numerator < denominator ifTrue: [^self reciprocal log2 negated].
    res := super log2.
    res isFinite ifTrue: [^res].
    ^numerator log2 - denominator log2

Float is IEEE754 double precision format (64 bits).
You can assume that Integer and Fraction asFloat are correctly rounded to nearest Float, tie to even
ln is neperian logarithm and will invoke log function from libm.
I only tried on mac os libm which seems accurate enough.

I also tried with i having more than Float emax bits in order to go through Fraction algorithm, but since denominator isPowerOfTwo, log2 is exact, and the subtraction too, so I don't believe that it can much decrease accuracy.

Upvotes: 1

Views: 63

Answers (0)

Related Questions