# Anthony Zaleski. Math 640 HW1. 24 Jan 2014.

#### 1.  Some examples from the Maple booklet:
(* 
#I turned all the command line stuff into a comment since otherwise I got an error when #reading this file into Maple.  Is there a better way to write this up?
p:=x^2+5*x+6;
p := x^2+5*x+6:
factor (p);
(x+3)*(x+2)

q:= x/abs(x);
q := x/abs(x):
assume (x>0):
simplify(q);
1

x='x':
readlib('rationalize'):
q:= 1/(1+sqrt(x)):
rationalize(%);
(-1+x^(1/2))/(-1+x)

subs(x=0, q);
1


solve(sin(z)=0, z);
0
solve(sin(z)= cos(z), z);
(1/4)*Pi


isprime(42);
false
ifactor(42);
``(2)*``(3)*``(7)
gcd(42, floor(Pi*10^10))*6;
42
*)


#### 2. Extending the M and K functions to work in base B:

Mbase:=proc(a,b,B) local d,a1,a2,b1,b2:

if a<B and b<B then

  RETURN(a*b):

fi:

#If a has 2d digits, let's write as a=B^d*a1+a2 and ditto for b
#where a1 and a2 have d digits

d:=log[B](max(a,b)):

if type(d,integer) then
d:=d+1:
else
 d:=ceil(d):
fi:

d:= ceil(d/2):

a2:=a mod B^d:  b2:=b mod B^d:

a1:=(a-a2)/B^d: b1:= (b-b2)/B^d:

# a*b=(B^d*a1 +a2)* (B^d*b1+b2)=
#(B^d)^2*a1*b1+B^d*a1*b2+a2*B^d*b1+a2*b2

Mbase(a1,b1,B)*(B^d)^2+ ( Mbase(a1,b2,B)+Mbase(a2,b1,B))*B^d + Mbase(a2,b2,B):

end:


#Karatsuba’s Method:

Kbase:=proc(a,b,B) local d,a1,a2,b1,b2:
option remember:

if a<B and b<B then

  RETURN(a*b):

fi:

#If a has 2d digits, let's write as a=B^d*a1+a2 and ditto for b
#where a1 and a2 have d digits

d:=log[B](max(a,b)):

if type(d,integer) then
d:=d+1:
else
 d:=ceil(d):
fi:

d:= ceil(d/2):

a2:=a mod B^d:  b2:=b mod B^d:

a1:=(a-a2)/B^d: b1:= (b-b2)/B^d:

# a*b=(B^d*a1 +a2)* (B^d*b1+b2)=
#(B^d)^2*a1*b1+B^d*a1*b2+a2*B^d*b1+a2*b2
#(B^d)^2*a1*b1+B^d*(a1*b2+a2*b1)+a2*b2
#a1*b2+a2*b2=(a1+a2)*(b1+b2)-a1*b1-a2*b2:

Kbase(a1,b1,B)*(B^2)^d+ (Kbase(a1+a2,b1+b2,B)- Kbase(a1,b1,B)-Kbase(a2,b2,B))*B^d + Kbase(a2,b2,B):

end:



#### 4. Newton’s Method.  The implementation here stops when successive approximations 
#become close enough.  It works with decimals.

NR := proc(f,x,x0,eps) local fprime,x1:

option remember:

fprime := subs(x=x0, diff(f,x)):

if fprime = 0 then
	printf("NR stopped due to a zero derivative.\n");
	RETURN(0):
else
	x1:= evalf(x0 - subs(x=x0,f)/fprime):
fi:


if evalf(abs(x1 - x0)) < evalf(eps) then
	RETURN(x1):
else
	NR(f,x,x1,eps):
fi:
	

end: