[Exercise] Arithmetical computation

Joe

This one is a bit tricky. Write a procedure that will solve the following equation:

x^x = 1000

The goal is to show a code that will compute the above value of x, as well as an approximation of this value to the 7th decimal digit.

Ra8

nice exercise

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

int main()
{
    double upper=5,lower=1,middle;
    cout<<setprecision(10);
    for(int i=0;i<100;i++)
    {
        middle=(lower+upper)/2.0;
        if(pow(middle,middle)>1000) upper=middle;
        else lower=middle;
    }
    cout << middle << " - " << pow(middle,middle);
    return 0;
}

output: 4.555535705

arithma

Touche Ra8. Beautiful solution.

Joe

Ra8, I like your solution, very elegant. You could maybe improve it by having the algorithm define initial upper and lower automatically.

My solution is a bit more complex (heavily inspired by my work on SICP):

x^x = 1000
x * log(x) = log(1000)
x = log(1000)/log(x)
x = f(x) where f: x -> log(1000)/log(x)

So the solution consists of finding the fixed point of the function f that maps x to log(1000)/log(x)

Here's a function for finding the fixed point of a procedure f in Scheme:

;; the tolerance can be changed
(define tolerance 0.00001)


(define (close-enough? v1 v2)
  (< (abs (- v1 v2)) tolerance))

;; first guess is arbitrary
(define (fixed-point f guess)
  (if (close-enough? guess (f guess))
      (f guess)
      (fixed-point f (f guess))))

So calculating the value of x consists of:

;; first-guess cannot be 1.0, so we take it to 2.0
(fixed-point (lambda (x) (/ (log 1000) (log x))
                  2.0)

arithma

@rahmu: Can you please explain the scheme code. I couldn't get how the fixed point finder works. Some pseudo code could really help :)

mesa177

I like the simplicity of Ra8's answer, but there is one thing I don't agree upon, and that is your ceiling limits. How do you know that your upper ceiling should be 5 and not 4? Also, why did you choose your lower ceiling to be 1? Why not 0.5?

When I read the problem, I thought of this (consider it a pseudo-code if you will):

x^x = ct where ct is any constant value different from 0 and x >0
xln(x) = ln(ct)
let y(x) = xln(x) - ln(ct) => 
according to Newton's method for derivation, x[i+1] = x[i] - y(x[i]) / y'(x[i])

=> x[i+1] = x[i] - x[i]*ln(x[i])/[ln(x[i]+1]

take x[i] = ct (since x^x = ct, x must be less than ct => safe high value to start truncation)
keep evaluating x[i+1] until y(x[i+1]) <= 0 

let z = x[i] for that value of x[i+1] and start process of zero-crossing detection 
for y(x) (i.e. y(x[i]) = 0)
start by decreasing value of z by 0.1 
(y(x[i]) is known to be positive => value of y(x[i]) must be decreased to reach zero-crossing)
keep evaluating y(x[i]) = x[i]ln(x[i]) - ln(ct) and decreasing value of z by 0.1 until y(x[i]) <=0

now start by increasing value of z by 0.01 and evaluating y(x[i]) until y(x[i]) >=0

keep repeating consecutive processes of decreasing and increasing values of 10^-n 
where n ranges from 3 to 7 (approximation must be accurate to the 7th decimal digit) 
as in prior step

When ct is set to 1000, x reaches a value of 4.555535705 as well.

Ra8

Since x^x will increase at a very high rate I think the step 1 is enough starting with 0 to find the lower and upper limits

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

int main()
{
    double upper=0,lower=0,middle,target=1000;
    while(true)
    {
        if(pow(lower,lower)<target)lower++;
        else
        {
            upper=lower--;
            break;
        }
    }
    cout<<setprecision(10);
    for(int i=0;i<100;i++)
    {
        middle=(lower+upper)/2.0;
        if(pow(middle,middle)>target) upper=middle;
        else lower=middle;
    }
    cout << middle << " - " << pow(middle,middle);
    return 0;
}

arithma

@mesa: Newton raphson for the win. I thought about that too. Engineers, yeah?

Joe

arithma wrote@rahmu: Can you please explain the scheme code. I couldn't get how the fixed point finder works. Some pseudo code could really help :)

I edited the code to make it more readable.

In short, the function will call itself as many times until f(x) ~ x.

Basically, consider delta the value:

delta= abs( f(x) - x )

the function keep trying to improve its result until delta is smaller than the defined tolerance.

How do I improve?

It happens that in this case (as in many other), you can find the fixed point of f by repeatedly calling the function on its own result:

let y = f(y)

y = f( f( f( f( f( f( f(... f(initial-guess)))))))))

You call f enough times for delta to be smaller than tolerance.

Initial-guess is arbitrary. The easiest case would take 1 (or 1.0 in Scheme). but in this case you can't because:
log(1000)/log(1) would be dividing by zero. So I took initial-guess=2.0.

Hope I was clear enough.

arithma

Crystal dude, I got stumped by searching where the fixed point was being iterated and forgot it could be recursed.