5
$\begingroup$

Given natural numbers $N, K, m, C$, with $3^{m/3}K>C$, I want to be able to write an algorithm to exactly compute the number

$$ \left\lceil \log_3 \left(\frac{N}{3^{m/3}K-C}\right) \right\rceil $$

Since I think I see how to exactly compare an integer against this number, there is the trivial method of "try larger numbers until one works" (maybe this could be sped up by doubling the number every time, and doing a binary search once an upper bound is found?).

However, I expect there are better methods. Can anyone point me to a reference on doing this sort of exact numeric computation, or at least tell me how to solve this particular problem in a non-terrible way?

EDIT: Since it seems a number of people have missed the point -- I want to compute this number exactly, and I want the output to be provably correct. Yes, in practice, doing it with doubles will get the right answer; I already have a version of the program that does that. But I want to write a version that is provably correct, and will, given sufficient memory, work even if the numbers input far exceed the range of a double (though I hope nobody ever does this). In short, any potential failure other than my computer running out of memory is unacceptable.

Since the program is written in Haskell, pointing me to an exact arithmetic library in Haskell that actually works and is provably correct, will suffice for an answer. But please do some sort of check that this is actually or at least likely to be true of the library you are proposing I use; for instance, the last one I tried, ERA, told me that $\lfloor 3-3^{-128}\rfloor=3$.

EDIT again: Some experimentation with what I thought would be a better method suggests maybe brute-force/binary search is the way to go after all. Oops? Well, I'm not certain yet, but my point is, if the answer really is "Yes, binary search really is the a good way to do this, no fancier way is needed," I will accept that as an answer, not demand some better one.

  • 0
    What language do you have in mind for the implementation?2012-04-27
  • 0
    Compute it inexactly using floating-point arithmetic, and use that as an initial guess to start the search?2012-04-27
  • 0
    @Rahul Narain: I want everything I'm doing to be provably correct, so I'd rather avoid floating-point. But using it as a guideline for searches may well work. Of course, there is the problem that in general the numbers may get larger than native floating point can handle.2012-04-27
  • 0
    @deoxygerbe: The program is written in Haskell. I suppose I could rewrite it in Mathematica or something but I'd rather not. And so far none of the exact arithmetic libraries I've found for Haskell have been quite satisfactory, which is why I'm asking this -- but perhaps you know one?2012-04-27
  • 0
    @HarryAltman: can you take natural logarithms and evaluate real^real easily in those libraries?2012-04-27
  • 0
    Also, I'd edit the title to mention Haskell -- which makes this much more interesting/nontrivial than doing it in C, which would look something like double H; int u; H=N/(pow(3,m/3)*K-C); u = ceil(log(H)/log(3));2012-04-27
  • 0
    @deoxygerbe: I think you've missed the point; I want the answer to be provably correct, and it's not clear to me that using doubles will be. After all, that C solution, rephrased, works perfectly well in Haskell; it's just not clear to me that it will always return the correct answer. While I expect that double precision will always be enough in practice, I don't see that I have any way of knowing this in advance.2012-04-27
  • 0
    @deoxygerbe: So far I've tried BigFloat, in which basic calculations sometimes fail to terminate, and ERA, where floor(3-3^-128) yielded 3 rather than 2. Again, I want exact computation; I already have a floating-point version. I know there are more out there, but that two already failed suggests maybe I should do this myself.2012-04-27
  • 0
    What is wrong with the method you suggested?2012-04-27
  • 0
    @copper.hat: It will work, certainly, but I expect there are better methods and I am hoping someone can point me to one.2012-04-27
  • 0
    How big are the numbers we're talking about here?2012-04-27
  • 0
    I presume you are doing something like finding the integer i that satisfies $3^{i-1} (3^{m/3}K-C) < N \leq 3^{i} (3^{m/3}K-C)$. That is, computing $\log_3$ is not 'part of the problem' as such.2012-04-27
  • 0
    Sorry, I should be clear -- the numbers will not actually in practice so big that doubles will fail. But do I have a guarantee of this beforehand? No. And I want the output to be *provably correct*. Therefore, I must cover the case when the numbers are as big as that, even if this case will not come up in practice. And if I need a separate, more general algorithm for that case, I may as well use that in general. The only acceptable failure in any hypothetical case is for my computer to run out of memory.2012-04-27
  • 0
    This might have something useful? http://stackoverflow.com/questions/2568446/the-best-cross-platform-portable-arbitrary-precision-math-library2012-04-27
  • 0
    Just as a minor note, 'provable correctness' for Haskell means something different than 'gives the right numerical answer'. What algorithmic criteria are you using for 'provably correct?2012-04-27
  • 0
    I have no idea what you're getting at. The choice of language is irrelevant except insofar as it restricts what is already done for me (since it's Haskell I have bigints) and what libraries I can use. It would be the same in any language, assuming I have bigints on hand (which here I do). By provably correct I mean it is guaranteed to give the correct answer on all possible inputs (or would be if my computer had infinite memory). I.e. there exists a proof of correctness, ideally one I could write down.2012-04-27
  • 0
    To put it another way: What I want is that if I were to write "Theorem: [statement] is true for numbers up to [bound]. Proof: Above I showed that [this algorithm] can check whether [statement] is true for a given number; I ran the algorithm for numbers up to [bound]; it always returned true." then people would accept that this is a proof. This requires that the algorithm provably can check the truth of the statement, for arbitrary input.2012-04-27
  • 0
    @HarryAltman: As an aside: "provable correctness" in the context of Haskell does have a very specific meaning, see: http://stackoverflow.com/questions/4077970/can-haskell-functions-be-proved-model-checked-verified-with-correctness-properti2012-04-27

2 Answers 2

1

Well, it turns out binary search was faster than high precision floating-point computations. Not a very interesting answer or mathematical answer, but that's what happened, so I'm going to mark this closed.

1

I have absolutely no experience using Haskell, but given a library supporting arbitrary precision floating point numbers and arbitrary size integers I'd do the following:

  • Define a function $\tt{check}$ that takes an integer as input and a bool as output with $$ {\tt check}(r) = {\tt true} ~~~ \Leftrightarrow ~~~ 3^{3r+m} K^3 \ge (N + 3^r)^3, $$ i.e. a function that, as you put it, exactly compares the integer to the desired output.
  • Compute the expression in floating point arithmetic with some starting precision, store the result as an arbitrarily large integer $r$.
  • Check whether ${\tt check}(r) = {\tt true}$ and ${\tt check}(r-1) = {\tt false}$. If that is the case, $r$ is proven to be the correct output. If not, double the precision and repeat.

I'm not sure if that counts as a non-terrible solution, but I think it should terminate relatively fast. Indeed, in most cases it should terminate on its first try, given a sensible starting precision. (I also assume that $3^{m/3}K-C > 0$, otherwise $\text{log}$ isn't well-defined anyway.)

If you want to avoid floating point arithmetic altogether, your idea of doubling integers and then using binary search should be fast enough. You only need something like $4 \text{log}( )$ checks to do this. I expect a floating-point-guess, integer-check method to be slightly faster, though.

  • 0
    What bugs me here is the "convert to integer" step. I mean, if there's not enough precision, it doesn't really represent a particular integer... but I guess that's why you double the precision and try again; whatever arbitrary-precision library I'd be using has some way of guessing one, and I don't really care what it is so long as it's eventually correct.2012-04-28