#Nathan Fox
#Homework 8
#I give permission for this work to be posted online

#Read procedures from class
read(`C9.txt`):
read(`C8.txt`):

Help:=proc():
    print(` ApplyNCg(f,x,n,a,b) , ApplyNCgg(f,x,n,a,b,K) `):
    print(` InitCompCounter(digits, start:=1) `):
    print(` IncCompCounter(counter) , BestNC(f,x,M) `):
    print(` ApplyGauss(f,x,n) , GaussVsNC(f,x,n) `):
end:

##PROBLEM 1##

#ApplyNCg(f,x,n,a,b): uses Newton-Cotes for approximating 
#int(f,x=a..b)
ApplyNCg:=proc(f, x, n, a, b) local u:
    return (b-a)*ApplyNC(subs(x=(b-a)*u+a, f), u, n):
end:

##PROBLEM 2##

#ApplyNCgg(f,x,n,a,b,K): approximates the integral of an
#expression f from x=a to x=b by chopping it into K equal
#intervals, applying ApplyNCg to each subinterval,
#and adding them up.
ApplyNCgg:=proc(f, x, n, a, b, K) local i:
    return add(ApplyNCg(f, x, n, a+(b-a)*i/K,
        a+(b-a)*(i+1)/K), i=0..K-1):
end:

##PROBLEM 3##

#Compound counter, each digit has a different cap,
#start from start
InitCompCounter:=proc(digs, start:=1) local i:
    return [[seq(start,i=1..nops(digs))],digs,start]:
end:

#Increment compound counter; return the new counter
#Return false if overflowed the counter
IncCompCounter:=proc(counter) local ctr,i:
    ctr:=counter:
    i:=1:
    while i<=nops(ctr[1]) and ctr[1][i]=ctr[2][i] do
        ctr[1][i]:=ctr[3]:
        i:=i+1:
    od:
    if i<=nops(ctr[1]) then
        ctr[1][i]:=ctr[1][i] + 1:
    else
        return false:
    fi:
    return ctr:
end:

#BestNC(f,x,M): inputs an expression f of x, a positive integer M,
#and outputs the pair of integers [n,K] such that M=n*K and 
#ApplyNCgg(f,x,n,0,1,K) gets the closest to the true int(f,x=0..1)
BestNC:=proc(f, x, M)
local pf, ctr, i, intval, n, K, champ, rec, nc:
    pf:=ifactors(M)[2]:
    intval:=evalf(Int(f, x=0..1)):
    champ:=[0,0]:
    rec:=infinity:
    ctr:=InitCompCounter([seq(pf[i][2], i=1..nops(pf))], 0):
    while ctr <> false do
        n:=mul(pf[i][1]^ctr[1][i], i=1..nops(pf)):
        K:=M/n:
        nc:=evalf(abs(intval-ApplyNCgg(f, x, n, 0, 1, K))):
        print(nc, [n, K]):
        if nc < rec then
            rec:=nc:
            champ:=[n, K]
        fi:
        ctr:=IncCompCounter(ctr):
    od:
    return champ:print(` IncCompCounter(counter) , BestNC(f,x,M) `):
end:

#I did Digits:=50 for this
#BestNC(1/(1+x^2), x, 32); returned [32, 1]
#BestNC(1/(1+x^4),x, 32); returned [32, 1]
#BestNC(exp(x), x,32); returned [32, 1]
#BestNC(log(1+x),x, 32); returned [32, 1]
#I'm noticing a theme here...

##PROBLEM 4##

#ApplyGauss(f,x,n): inputs an expression f
#in the variable x, and uses the Gaussian Quadrature
#rule with n points
#to find an APPROXIMATION to int(f,x=0..1):
ApplyGauss:=proc(f,x,n) local L, X, H, i:
    L:=FindGaussSmart(n):
    X:=L[1]:
    H:=L[2]:
    return add(H[i]*subs(x=X[i], f), i=1..nops(X)):
end:

##PROBLEM 5##

#Compare Gauss to Newton-Cotes
#Fight to the death!
GaussVsNC:=proc(f, x, n) local intval, gauss, nc:
    intval:=evalf(Int(f, x=0..1)):
    gauss:=evalf(ApplyGauss(f, x, n)):
    nc:=evalf(ApplyNC(f, x, n)):
    printf("#f=%a,\tn=%d\n", f, n):
    printf("#Integral:\t%f\n#Gauss:\t%f\n#Newton-Cotes:\t%f\n",
        intval, gauss, nc):
    if gauss < intval then
        printf("#Gauss underestimated by %f\n", intval - gauss):
    else
        printf("#Gauss overestimated by %f\n", gauss - intval):
    fi:
    if nc < intval then
        printf("#Newton-Cotes underestimated by %f\n",
            intval - nc):
    else
        printf("#Newton-Cotes overestimated by %f\n", nc - intval):
    fi:
    if abs(intval - gauss) < abs(intval - nc) then
        printf("#Gauss wins!\n"):
    elif abs(intval - gauss) > abs(intval - nc) then
        printf("#Newton-Cotes wins!\n"):
    else #Happens with probability 0
        printf("#It's a tie!\n"):
    fi:
end:

#Did Digits:=50

#f=1/(1+x^2),	n=1
#Integral:	0.785398
#Gauss:	0.800000
#Newton-Cotes:	0.750000
#Gauss overestimated by 0.014602
#Newton-Cotes underestimated by 0.035398
#Gauss wins!
#f=1/(1+x^2),	n=2
#Integral:	0.785398
#Gauss:	0.786885
#Newton-Cotes:	0.783333
#Gauss overestimated by 0.001487
#Newton-Cotes underestimated by 0.002065
#Gauss wins!
#f=1/(1+x^2),	n=3
#Integral:	0.785398
#Gauss:	0.785267
#Newton-Cotes:	0.784615
#Gauss underestimated by 0.000131
#Newton-Cotes underestimated by 0.000783
#Gauss wins!
#f=1/(1+x^2),	n=4
#Integral:	0.785398
#Gauss:	0.785403
#Newton-Cotes:	0.785529
#Gauss overestimated by 0.000005
#Newton-Cotes overestimated by 0.000131
#Gauss wins!
#f=1/(1+x^2),	n=5
#Integral:	0.785398
#Gauss:	0.785398
#Newton-Cotes:	0.785470
#Gauss underestimated by 0.000000
#Newton-Cotes overestimated by 0.000071
#Gauss wins!


#f=1/(1+x^4),	n=1
#Integral:	0.866973
#Gauss:	0.941176
#Newton-Cotes:	0.750000
#Gauss overestimated by 0.074203
#Newton-Cotes underestimated by 0.116973
#Gauss wins!
#f=1/(1+x^4),	n=2
#Integral:	0.866973
#Gauss:	0.859522
#Newton-Cotes:	0.877451
#Gauss underestimated by 0.007450
#Newton-Cotes overestimated by 0.010478
#Gauss wins!
#f=1/(1+x^4),	n=3
#Integral:	0.866973
#Gauss:	0.867518
#Newton-Cotes:	0.871071
#Gauss overestimated by 0.000545
#Newton-Cotes overestimated by 0.004098
#Gauss wins!
#f=1/(1+x^4),	n=4
#Integral:	0.866973
#Gauss:	0.866956
#Newton-Cotes:	0.866425
#Gauss underestimated by 0.000017
#Newton-Cotes underestimated by 0.000548
#Gauss wins!
#f=1/(1+x^4),	n=5
#Integral:	0.866973
#Gauss:	0.866971
#Newton-Cotes:	0.866674
#Gauss underestimated by 0.000002
#Newton-Cotes underestimated by 0.000299
#Gauss wins!


#f=exp(x),	n=1
#Integral:	1.718282
#Gauss:	1.648721
#Newton-Cotes:	1.859141
#Gauss underestimated by 0.069561
#Newton-Cotes overestimated by 0.140859
#Gauss wins!
#f=exp(x),	n=2
#Integral:	1.718282
#Gauss:	1.717896
#Newton-Cotes:	1.718861
#Gauss underestimated by 0.000385
#Newton-Cotes overestimated by 0.000579
#Gauss wins!
#f=exp(x),	n=3
#Integral:	1.718282
#Gauss:	1.718281
#Newton-Cotes:	1.718540
#Gauss underestimated by 0.000001
#Newton-Cotes overestimated by 0.000258
#Gauss wins!
#f=exp(x),	n=4
#Integral:	1.718282
#Gauss:	1.718282
#Newton-Cotes:	1.718283
#Gauss underestimated by 0.000000
#Newton-Cotes overestimated by 0.000001
#Gauss wins!
#f=exp(x),	n=5
#Integral:	1.718282
#Gauss:	1.718282
#Newton-Cotes:	1.718282
#Gauss underestimated by 0.000000
#Newton-Cotes overestimated by 0.000000
#Gauss wins!


#f=ln(x+1),	n=1
#Integral:	0.386294
#Gauss:	0.405465
#Newton-Cotes:	0.346574
#Gauss overestimated by 0.019171
#Newton-Cotes underestimated by 0.039721
#Gauss wins!
#f=ln(x+1),	n=2
#Integral:	0.386294
#Gauss:	0.386595
#Newton-Cotes:	0.385835
#Gauss overestimated by 0.000301
#Newton-Cotes underestimated by 0.000460
#Gauss wins!
#f=ln(x+1),	n=3
#Integral:	0.386294
#Gauss:	0.386300
#Newton-Cotes:	0.386084
#Gauss overestimated by 0.000006
#Newton-Cotes underestimated by 0.000211
#Gauss wins!
#f=ln(x+1),	n=4
#Integral:	0.386294
#Gauss:	0.386294
#Newton-Cotes:	0.386288
#Gauss overestimated by 0.000000
#Newton-Cotes underestimated by 0.000006
#Gauss wins!
#f=ln(x+1),	n=5
#Integral:	0.386294
#Gauss:	0.386294
#Newton-Cotes:	0.386291
#Gauss overestimated by 0.000000
#Newton-Cotes underestimated by 0.000004
#Gauss wins!

#Verdict: Gauss beats Newton into submission!