# file "Rillinois2"  (related to Q3 of HW#2)

x = c(3,4,8,10,11)

y = c(5,10,21,35,40)

f = function(beta) { sum((y-x^beta)^2) }

g = fp = function(beta) { sum(-2*(y-x^beta)*x^beta*log(x)) }

tolerance = 0.000001

avec = c(0)
bvec = c(2)
ga = g(avec[1])
gb = g(bvec[1])
lastupdated = 3  # 1 for a, 2 for b

for (i in 1:100) {
    a = avec[length(avec)]
    b = bvec[length(bvec)]
    c = a - ga * (a-b) / (ga - gb)
    if (g(c)==0) {
	cat("answer equals ", c, "\n")
	break
    } else if ( abs(b-a) < tolerance ) {
	cat("answer is within", tolerance, "of", a, "and", b, "\n")
	break
    } else if (g(c) * g(b) < 0) {
	avec = c(avec, c)
	ga = g(c)
	if (lastupdated==1)
	    gb = gb / 2
	lastupdated = 1
    } else {
	bvec = c(bvec, c)
	gb = g(c)
	if (lastupdated==2)
	    ga = ga / 2
	lastupdated = 2
    }
}

print(avec)
print(bvec)