2
$\begingroup$

Setup. Suppose we are given a function $$f:\mathbb N\to\{\text{False},\text{True}\}$$ such that $$f(n)=\text{True}\implies f(n+1)=\text{True}$$ and such that $f(n)=\text{True}$ for some $n$ large enough.

In natural language. The function $f$ imposes a condition on the natural numbers which is fulfilled once $n$ is large enough.

My question 1. How can I, for any given $f$, find the smallest $n$ such that $f(n)=\text{True}$?

A first idea would be to start with $n=1$ and to increment $n$ by one until $f(n)$ is True. However, this is fairly slow. The next step is a "binary search" algorithm (see below).

My question 2. How can I implement a good algorithm in Mathematica?

Here is an example of the first algorithm that I thought of, implemented in Python:

def bin_search(cond):
    n = 1
    while not cond(n):
        n *= 2
    lower_bound = n//2
    upper_bound = n
    middle = (lower_bound + upper_bound)//2
    while upper_bound - lower_bound > 1:
        if cond(middle):
            upper_bound = middle
        else:
            lower_bound = middle
        middle = (lower_bound + upper_bound)//2
    return upper_bound

For example, one such condition would be $$f(n)=[H_n\geq 10],$$

where $$H_n=\sum_{i=1}^n \frac 1i$$ is the $n$th harmonic number.

EDIT: A remark about the harmonic numbers: In fact I wanted to compare the runtime of a good brute force method against an intelligent method for this problem. For solving $H_n\geq 10$, it suffices to find the smallest integer $n$ such that $$\ln(n)+\gamma\geq 10,$$ where $\gamma$ denotes the Euler-Mascheroni constant. So we can find $$n=[ \exp(10-\gamma)],$$ where $[\cdot]$ denotes the nearest integer function.

$\endgroup$
8
  • 2
    $\begingroup$ 1. You are asking two different questions: (Q1) an algorithm derivation one, and (Q2) a code translation one. Which one do you want? 2. If Q1, then there is no need for ugly Python code... :) $\endgroup$ Dec 3, 2019 at 13:26
  • $\begingroup$ @AntonAntonov Since I am looking at general $f$, my question is: What is a good implementation in Mathematica of a good algorithm which is able to find the smallest $n$ such that ...? I don't want to use any specific properties of $f$ except the two mentioned above... My Python code was just an example of the simplest algorithm I could come up with, namely, a sort of binary search in order to find the minimal $n$ $\endgroup$ Dec 3, 2019 at 13:32
  • $\begingroup$ @AntonAntonov By the way, what is wrong with the Python implementation? 😯 $\endgroup$ Dec 3, 2019 at 18:59
  • $\begingroup$ You are aware, of course, that the problem has no solution in general, right? One can always construct a sequence growing slower than the $2^n$ built in the algorithm. $\endgroup$
    – yarchik
    Dec 3, 2019 at 19:15
  • 2
    $\begingroup$ @yarchik Actually, based on your Input, I will consider asking a part of this question on the StackExchange for theoretical Computer Science. That might reduce the confusion $\endgroup$ Dec 3, 2019 at 19:34

2 Answers 2

4
$\begingroup$

Here is one implementation of your algorithm:

firstTrue[f_] := next[f, 1, Infinity]

next[f_, last_, cur_] := If[last + 1 == cur,
    If[f[last], last, cur],
    If[f[last], next[f, dec[last], last], next[f, inc[last], cur]]
]

inc[n_] := With[{e = IntegerExponent[n, 2]},
    If[n == 2^e, 2^(e+1), n + 2^(e-1)]
]

dec[n_] := With[{e = IntegerExponent[n, 2]},
    If[n == 2^e, 2^(e-1) + 2^(e-2), n - 2^(e-1)]
]

Your example:

firstTrue[HarmonicNumber[#] >= 10&]

12367

Of course, for analytic functions like HarmonicNumber there are much better algorithms to do this.

$\endgroup$
2
  • $\begingroup$ Thank you! About HarmonicNumber: I will add some motivation; of course you are right $\endgroup$ Dec 3, 2019 at 18:20
  • $\begingroup$ Ok I have added some motivation to the Harmonic Number problem. In my opinion it is a very nice solution using Euler-Mascheroni $\endgroup$ Dec 3, 2019 at 18:23
2
$\begingroup$

As Carl Woll says:

Of course, for analytic functions like HarmonicNumber there are much better algorithms to do this.

Here is some code related to that statement -- it uses one of the universal WL functions:

Clear[H];
H[n_] := Sum[1/i, {i, n}];

FindArgMin[{H[n], H[n] >= 10, n > 0}, n]
(* {12366.5} *)

N@H[Ceiling[%[[1]]]]
(* 10. *)

N@H[Floor[%%[[1]]]]
(* 9.99996 *)
$\endgroup$
1
  • $\begingroup$ This is interesting because it simply uses numerical methods. A caveat is the innacuracy: For example, FindArgMin[{H[n], H[n] >= 30, n > 0}, n, WorkingPrecision -> 1000, AccuracyGoal -> 1000, PrecisionGoal -> 1000] returns $$6.000022445096458984375\cdot 10^{12}$$ and the correct result is $$6000022499693$$ But it is indeed a very quick way to get close to the correct solution $\endgroup$ Dec 3, 2019 at 18:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.