Reputation: 9572
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