for large x, log (x::Integer) :: Double
i'm looking for an accurate way to take the log of a very large integer, for example: let x :: Integer = 1301427272151881160612765560226881966218101403436917787184856303672382623256898455416763978959067300249652773943715743032733292602624834984761739233232794619193611954735720284761058146899246611113236700853600891798968920775344491685185906922596026543915321367577774522912315930144523472702386240645993859368230855941019371447058664115974032571881072431604651385520393674840678811793554266595013773947434115579588912967969680150473258236727830867832149867100437142705476716669039640252677955201589378051836112800268367331455296715904387732836350613539218249950829555418397197909288345303407194983545308212828662999623279222913080214196287140117582811769188486693208227570257136851945943408206281672555558289460256867016896063334140640075708083581866297494610834545554864846306383014549439540479675828018496049574066533167553894586573246931377586176000 this # is approximately 10^850. not surprisingly, if i do: log (fromIntegral x) :: Double i get "Infinity". i can think of several ways to try to combat this: (1) divide x by something large (10^850) and then take the log of the integer part (requires finding that x ~ 10^850 by some search method); (2) multiply x by something large (10^10) and then take the log of this *as an Integer* and then convert that to a double. has anyone encountered this problem before? surely there must be a "good way" to do this... - hal -- Hal Daume III | hdaume@isi.edu "Arrest this man, he talks in maths." | www.isi.edu/~hdaume
On 30/6/04 14:47, Hal Daume III wrote:
i'm looking for an accurate way to take the log of a very large integer, for example:
let x :: Integer = 130142727215188116061276556022688196621810140343691778718485630367238262325689 845541676397895906730024965277394371574303273329260262483498476173923323279461 919361195473572028476105814689924661111323670085360089179896892077534449168518 590692259602654391532136757777452291231593014452347270238624064599385936823085 594101937144705866411597403257188107243160465138552039367484067881179355426659 501377394743411557958891296796968015047325823672783086783214986710043714270547 671666903964025267795520158937805183611280026836733145529671590438773283635061 353921824995082955541839719790928834530340719498354530821282866299962327922291 308021419628714011758281176918848669320822757025713685194594340820628167255555 828946025686701689606333414064007570808358186629749461083454555486484630638301 4549439540479675828018496049574066533167553894586573246931377586176000
this # is approximately 10^850.
Hal, Based on the old manual way of doing this, I would first obtain the log10 value. This is easiest if you can convert the number to a string: - if the original number is an integer, the string length gives you the exponent, which in log10 forms the integer part of the result - then take the first few digits of the string, concat '0.' in front, convert to a double, and obtain the log10 of that, which forms the decimal part of the result. And given the log10, if you want loge or ln then all you have to do is multiply by loge(10), if I remember correctly. Howard. Dr Howard Oakley The Works columnist for MacUser magazine (UK) http://www.macuser.co.uk/
On 30/6/04 19:56, Howard Oakley wrote:
- if the original number is an integer, the string length gives you the exponent, which in log10 forms the integer part of the result
Well, almost. The string length is actually 1 greater than the integer part of the result, as log10(1) = 0.0, and log10(10) = 1.0, if the dusty room at the back of my mind is not deceiving me :-) Howard. Dr Howard Oakley The Works columnist for MacUser magazine (UK) http://www.macuser.co.uk/
At 06:47 30/06/04 -0700, Hal Daume III wrote:
i'm looking for an accurate way to take the log of a very large integer,
This may not help, but just in case... there is an effective (approximate) algorithm for computing a log-gamma function, documented in "Numerical Recipies" by Press/Flannery/Teukolsky/Vetterling. If you want more info, I have an implementation in Mathematica somewhere, and could dig out a more specific reference. The alternative way I'd consider involves building a power series of (2^(2^i)) and performing a compare-and-divide-and-sum so that you get n and x in 2^n*x, for x in the range (0.0..1.0) (or any other range that you may choose). #g -- At 06:47 30/06/04 -0700, Hal Daume III wrote:
i'm looking for an accurate way to take the log of a very large integer, for example:
let x :: Integer = 1301427272151881160612765560226881966218101403436917787184856303672382623256898455416763978959067300249652773943715743032733292602624834984761739233232794619193611954735720284761058146899246611113236700853600891798968920775344491685185906922596026543915321367577774522912315930144523472702386240645993859368230855941019371447058664115974032571881072431604651385520393674840678811793554266595013773947434115579588912967969680150473258236727830867832149867100437142705476716669039640252677955201589378051836112800268367331455296715904387732836350613539218249950829555418397197909288345303407194983545308212828662999623279222913080214196287140117582811769188486693208227570257136851945943408206281672555558289460256867016896063334140640075708083581866297494610834545554864846306383014549439540479675828018496049574066533167553894586573246931377586176000
this # is approximately 10^850.
not surprisingly, if i do:
log (fromIntegral x) :: Double
i get "Infinity".
i can think of several ways to try to combat this: (1) divide x by something large (10^850) and then take the log of the integer part (requires finding that x ~ 10^850 by some search method); (2) multiply x by something large (10^10) and then take the log of this *as an Integer* and then convert that to a double.
has anyone encountered this problem before? surely there must be a "good way" to do this...
- hal
-- Hal Daume III | hdaume@isi.edu "Arrest this man, he talks in maths." | www.isi.edu/~hdaume
_______________________________________________ Haskell mailing list Haskell@haskell.org http://www.haskell.org/mailman/listinfo/haskell
------------ Graham Klyne For email: http://www.ninebynine.org/#Contact
On Wed, 30 Jun 2004 06:47:49 -0700 (PDT), Hal Daume III <hdaume@ISI.EDU> wrote:
i'm looking for an accurate way to take the log of a very large integer, for example:
The standard algorithm works along these lines: (Assume for the purposes of discussion that your numbers are represented in some kind of binary notation.) 0) Call your integer i. 1) Compute d, the smallest power of 2 that is greater than i. 2) Compute the floating-point quotient q = i / d. This will be a value between 0.5 and 1.0. 3) Compute n, the base-2 log of d. (This is trivial, since d is a power of 2.) 4) Use your CPU's built-in floating-point unit to compute g, the base-2 log of q. (For an Intel CPU, see the FYL2X instruction. If q is very close to 1, you may want to use the FYL2XP1 instruction instead to preserve more significant bits; in that case, pass it q - 1 rather than q.) 5) The base-2 log of i then simply n + g, and you can convert that to another base b by multiplying it by the base-2 log of b (or, equivalently, by dividing it by the base-b log of 2). If your integer is stored in some other fashion (i.e., not binary), then you'll want to adjust the divisor and the logarithm base appropriately. The idea is that you want to choose your initial divisor so that it is a power of the base, where the value of that base is chosen to make the process of computing the quotient q efficent and accurate. Steve Schafer Fenestra Technologies Corp http://www.fenestra.com/
Steve Schafer wrote:
4) Use your CPU's built-in floating-point unit to compute g, the base-2 log of q. (For an Intel CPU, see the FYL2X instruction. If q is very close to 1, you may want to use the FYL2XP1 instruction instead to preserve more significant bits; in that case, pass it q - 1 rather than q.)
My colleague from the Java-for-x86 world reminds me that the floating-point transcendentals on x86 are chock full of little infelicities which limit the conditions under which they may be used. My feeling is you want to understand these infelicities before you start using the instructions in your code. (For example, you should be able to explain, and correct, Mr. Shafer's explanation of FYL2X versus FYL2XP1 above when 0.5 <= q <= 1.0. I haven't quite figured it out myself...) Here be dragons. For the rest of us, there are very nice floating-point libraries. If you're lucky, their authors actually did the above work... :-) -Jan
participants (5)
-
Graham Klyne -
Hal Daume III -
Howard Oakley -
Jan-Willem Maessen - Sun Labs East -
Steve Schafer