VicTT <victhetraitor at yahoo.com> writes:
> I'll keep this short in the interest of preserving your time: I'm
> trying to find min(k) such that, given integers: x,r and rational
> epsilon, abs(x^(p/(2^k)-x^(1/r))<epsilon.
>
> Then, given integers: x,r,min(k) and rational epsilon, find p in the
> above inequality. How would I do this in Maxima? (first step is, I
> assume, rewriting this as an equality *somehow*) Also, if anyone knows
> any 3rd party add-ons that could solve this, I'd appreciate the
> contribution. I'm missing something, I just can't tell what. (This is
> neither homework, nor a task I've devoted less than 1 hour towards
> achieving.)
>
> Thank you in advance,
> Pirvu Paul Daniel
Ok, so this is a maths answer more than a Maxima answer, but
whatever. Firstly, as Laurent pointed out, I think you're missing a
closing parenthesis:
(%i1) ineq: abs(x^(p/(2^k))-x^(1/r))<epsilon;
! p !
! --!
! k!
! 1/r 2 !
(%o1) !x - x ! < epsilon
If I understand correctly, what you're trying to do is approximate 1/r
by p/2^k sufficiently well such that the above inequality holds. I'm
also going to assume that x is positive throughout. If not, you'll have
to tell me what you mean by x^(1/r).
So let's suppose for a second that we've chosen p/2^k such that
p/2^k = 1/r*(1+delta)
for some small delta. Substituting in gives
(%i2) ineq, p/2^k = 1/r*(1+delta);
! delta + 1 !
! --------- !
! r 1/r!
(%o2) !x - x ! < epsilon
(%i3) factor(%);
! 1/r! ! delta/r !
(%o3) !x ! !x - 1! < epsilon
And we should probably tell Maxima that x is positive, to get rid of the
pesky absolute value signs:
(%i5) %o3, ev;
1/r ! delta/r !
(%o5) x !x - 1! < epsilon
(The "ev" tells Maxima to re-evaluate the expression, which causes it to
notice the sign of x). This is starting to look a bit easier! Now we
need to find an allowable range for delta. Divide through by the first
bit:
(%i6) ineq2: map(lambda([u], u/x^(1/r)), %);
! delta/r ! epsilon
(%o6) !x - 1! < -------
1/r
x
(The weird "map" syntax is because something like (a < b) + 1 doesn't
expand to a+1 < b+1. I think I'll send a separate message to the mailing
list asking whether we should maybe expand them.)
Note that x is a positive integer, so x^(1/r) is positive and the sign
of the inequality doesn't change.
First, let's suppose that delta > 0. Unfortunately, telling Maxima that
delta is positive doesn't seem to be enough to get it to throw away the
absolute value signs, so we'll have to do the following:
(%i7) assume(x^(delta/r) > 1);
delta/r
(%o7) [x > 1]
(%i8) ineq2, ev;
delta/r epsilon
(%o8) x - 1 < -------
1/r
x
(%i9) map(lambda([u], u+1), %);
delta/r epsilon
(%o9) x < ------- + 1
1/r
x
(%i10) map(log, %);
delta log(x) epsilon
(%o10) ------------ < log(------- + 1)
r 1/r
x
x was a positive integer, so is greater than 1 and so log(x) >= 0. This
inequality is true for all delta if x = 1, so let's assume x > 1. Then
log(x) is positive. You also didn't say that r is positive, but I think
it should be(!), so let's assume that too so we can divide through to
give:
(%i11) map(lambda([u], u*r/log(x)), %);
epsilon
r log(------- + 1)
1/r
x
(%o11) delta < ------------------
log(x)
Right, that's our upper bound on delta:
(%i12) delta_ub: rhs(%);
epsilon
r log(------- + 1)
1/r
x
(%o12) ------------------
log(x)
Now we look at delta < 0, so we'll have to forget what we assumed:
(%i13) forget(x^(delta/r) > 1);
delta/r
(%o13) [x > 1]
(%i14) assume(x^(delta/r) < 1);
delta/r
(%o14) [x < 1]
(%i15) ineq2, ev;
delta/r epsilon
(%o15) 1 - x < -------
1/r
x
Do the same trick to rearrange and take logs, to get:
(%i16) map(lambda([u], u+x^(delta/r)-rhs(%)), %);
epsilon delta/r
(%o16) 1 - ------- < x
1/r
x
(%i17) map(log, %);
epsilon delta log(x)
(%o17) log(1 - -------) < ------------
1/r r
x
(%i18) map(lambda([u], u*r/log(x)), %);
epsilon
r log(1 - -------)
1/r
x
(%o18) ------------------ < delta
log(x)
(%i19) delta_lb: lhs(%);
epsilon
r log(1 - -------)
1/r
x
(%o19) ------------------
log(x)
(Notice that this is negative since the argument to the log is less than
1)
Now we've got some rather unpleasant bounds for delta. The next problem
is to translate them into a minimum k. Solving our original equation
that defined delta you get:
(%i20) solve(p/2^k = 1/r*(1+delta), delta);
k
p r - 2
(%o20) [delta = --------]
k
2
Now, remember that you can choose p as required, so we can guarantee
that the top of the fraction is at most r, giving that abs(delta) <
2^(-k). For delta_ub, you get:
(%i21) log(2^(-k))/log(2) < log(delta_ub)/log(2);
epsilon
r log(------- + 1)
1/r
x
log(------------------)
log(x)
(%o21) - k < -----------------------
log(2)
Let's tidy up the ridiculous tower of logarithms and fix the sign:
(%i22) radcan(%);
1/r
log(r log(x + epsilon) - log(x)) - log(log(x))
(%o22) - k < -------------------------------------------------
log(2)
(%i23) map(lambda([u], u+k-rhs(%)), %);
1/r
log(r log(x + epsilon) - log(x)) - log(log(x))
(%o23) - ------------------------------------------------- < k
log(2)
(%i24) k_lb_pos: first(%);
1/r
log(r log(x + epsilon) - log(x)) - log(log(x))
(%o24) - -------------------------------------------------
log(2)
Now do it all again, but for delta_lb (but use -delta_lb, so we get
positive stuff to take logs of):
(%i25) log(2^(-k))/log(2) < log(-delta_lb)/log(2);
epsilon
r log(1 - -------)
1/r
x
log(- ------------------)
log(x)
(%o25) - k < -------------------------
log(2)
(%i26) radcan(%);
1/r
log(log(x) - r log(x - epsilon)) - log(log(x))
(%o26) - k < -------------------------------------------------
log(2)
(%i27) map(lambda([u], u+k-rhs(%)), %);
1/r
log(log(x) - r log(x - epsilon)) - log(log(x))
(%o27) - ------------------------------------------------- < k
log(2)
(%i28) k_lb_neg: first(%);
1/r
log(log(x) - r log(x - epsilon)) - log(log(x))
(%o28) - -------------------------------------------------
log(2)
Now we have a (somewhat!) explicit formula for k:
(%i29) min_k: min(k_lb_pos, k_lb_neg);
1/r
log(log(x) - r log(x - epsilon)) - log(log(x))
(%o29) min(- -------------------------------------------------,
log(2)
1/r
log(r log(x + epsilon) - log(x)) - log(log(x))
- -------------------------------------------------)
log(2)
Now I confess that I cheated and plotted some values:
plot2d([k_lb_pos/k_lb_neg], [x, 2, 100]), r=30, epsilon=0.001;
changing r and epsilon and noticed that k_lb_pos > k_lb_neg. To prove
that this wasn't just a fluke, you can take their difference, then
Taylor expand in epsilon:
(%i30) taylor(k_lb_pos - k_lb_neg, epsilon, 0, 3);
3
epsilon epsilon
(%o30)/T/ ----------- + ---------------- + . . .
1/r 1/r 3
x log(2) 4 (x ) log(2)
This hopefully convinces you that k_lb_pos is always larger for at least
small epsilon. Right, so:
(%i31) min_k: k_lb_pos;
1/r
log(r log(x + epsilon) - log(x)) - log(log(x))
(%o31) - -------------------------------------------------
log(2)
All that remains is to find the best value of p given k,x and r. We need
to minimise p*r-2^k, so the answer will be one of
floor(2^k/r), floor(2^k/r)+1
Here's code to choose:
(%i32) find_p(k,x,r) :=
block([a: floor(2^k/r), b: floor(2^k/r)+1],
if is(abs(r*a - 2^k) < abs(r*b - 2^k)) then a else b)$
Finally, wrap everything up:
find_kp(x, r, epsilon) :=
block([k: ceiling(subst(['x = x, 'r = r, 'epsilon=epsilon], min_k)),
p],
p: find_p(k, x, r),
[k, p, float(abs(x^(p/2^k)-x^(1/r)))])$
Some examples:
(%i34) find_kp(5, 3, 10^(-3));
(%o34) [10, 341, 8.956312902090868e-4]
(%i35) find_kp(5, 4, 10^(-1));
(%o35) [3, 2, 0.0]
(%i36) find_kp(5, 4, 10^(-10));
(%o36) [33, 2147483648, 0.0]
Note that sometimes we're overly pessimistic. For the last example, I
should probably have returned [3,2,0] again. Fixing this is an exercise
for the reader...
Rupert
PS: Yep, this is a ridiculously long answer, but I was enjoying myself
working it out!
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 315 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20130131/a434e897/attachment.pgp>;