luc(n) := block([L: 4], for i: 1 thru n - 1 do L: L^2 - 2, L);

merpp(n) := block([m: 2^n - 1, x: 4],
	for i: 1 thru n - 2 do x: remainder(x^2 - 2, m),
	if x = 0
		then return("prime")
		else return("composite")
);

prup(a, n) := block([c: 0],
	for i: 1 thru n do
		if primep(i) then
			if primep(ru(a, i)) then c: c + 1, c
);

rmmpt(f) := block([fx, fxx, s, t, Lz: [], Lp: [], Ln : [], realonly: true],
	fx: diff(f, x),
	fxx: diff(f, x, 2),
	s: algsys([fx], [x]),
	if length(s) = 0 then
		return("Stationary Point Not Found")
	else (
		for i: 1 thru length(s) do (
			t: subst(s[i], fxx),
			if t = 0 then
				Lz: cons(s[i][1], Lz)
			elseif t > 0 then
				Lp: cons(subst(s[i], [x, f]), Lp)
			else
				Ln: cons(subst(s[i], [x, f]), Ln)
		),
		if length(Lz) # 0 then (
			print("UNKNOWN:"),
			for i in Lz do disp(i)
		),
		if length(Lp) # 0 then (
			print("Relative Minimum Point:"),
			for i in Lp do disp(i)
		),
		if length(Ln) # 0 then (
			print("Relative Maximum Point:"),
			for i in Ln do disp(i)
		)
	)
);

rmmpt_two(f) := block([fx, fy, fxx, fxy, fyy, J, s, t, u, Lz: [], Lp: [], Ln : [], realonly: true],
	fx: diff(f, x), fy: diff(f, y),
	fxx: diff(f, x, 2), fxy: diff(f, x, 1, y, 1), fyy: diff(f, y, 2),
	J: fxx * fyy - fxy^2,
	s: algsys([fx, fy], [x, y]),
	if length(s) = 0 then
		return("Critical Point Not Found")
	else (
		for i: 1 thru length(s) do (
			t: subst(s[i], J),
			if t = 0 then
				Lz: cons(subst(s[i], [x, y, f]), Lz)
			elseif t > 0 then (
				u: subst(s[i], fxx),
				if u > 0 then
					Lp: cons(subst(s[i], [x, y, f]), Lp)
				else
					Ln: cons(subst(s[i], [x, y, f]), Ln)
			)
		),
		if length(Lz) # 0 then (
			print("UNKNOWN:"),
			for i in Lz do disp(i)
		),
		if length(Lp) # 0 then (
			print("Relative Minimum Point:"),
			for i in Lp do disp(i)
		),
		if length(Ln) # 0 then (
			print("Relative Maximum Point:"),
			for i in Ln do disp(i)
		)
	)
);

eratosthenes(N) := block([p: 2, L],
	L: makelist(1, i, 1, N),
	while(p <= sqrt(N)) do (
		if(L[p] = 1) then (
			k: p,
			while(k + p <= N) do (
				k: k + p,
				L[k]: 0
			)
		),
		p: p + 1
	),
	for i: 2 thru N do
		if L[i] = 1 then
			print(i)
);

factor_e(M) := block([p: 2, N: floor(sqrt(M)), L, S: []],
	load("numericalio"),
	L: makelist(1, i, 1, N),
	while(p <= sqrt(N)) do (
		if(L[p] = 1) then (
			k: p,
			while(k + p <= N) do (
				k: k + p,
				L[k]: 0
			)
		),
		p: p + 1
	),
	with_stdout("eratosthenes_output.txt",
		for i: 2 thru N do
			if L[i] = 1 then
				print(i)
	),
	for i in read_list("eratosthenes_output.txt") do
		while(remainder(M, i) = 0) do (
			S: cons(i, S),
			M: M / i
		),
	if M # 1 then S: cons(M, S),
	S
);

fib_s(n) := block([k],
	if n = 0 or n = 1 then
		return(n)
	elseif evenp(n) then (
		k: n/2,
		return(fib_2(k) * (fib_2(k) + 2 * fib_2(k-1)))
	) else (
		k: (n + 1) / 2,
		return(fib_2(k)^2 + fib_2(k-1)^2)
	)
);

fib_t(n) := block([c, a: 0, b: 1],
	for i: 1 thru n - 1 do (
		c: a + b,
		a: b,
		b: c
	),
	return(c)
);

binary(a, n) := block([x: 1],
	while ( n # 0 ) do (
		if oddp(n) then (
			x: x * a,
			n: n - 1
		),
		a: a * a,
		n: n / 2
	),
	return(x)
);

binary_fib(n) := block([X: ident(2), F: matrix([1, 1], [1, 0])],
	if (n = 0 or n = 1) then
		retuen(n)
	else (
		n: n - 1,
		while ( n # 0 ) do (
			if oddp(n) then (
				X: X . F,
				n: n - 1
			),
			F: F . F,
			n: n / 2
		),
		return(X[1, 1])
	)
);

count(file) := block([s, L: makelist(0, i, 1, 26)],
	for i in read_list(file) do (
		s: sdowncase(printf(false, "~a", i)),
		for j in charlist(s) do
			if alphacharp(j) then L[cint(j) - 96]: L[cint(j) - 96] + 1
	),
	return(L)
);

count_g(file, n) := block([a: [], s: [], m, L: makelist(0, i, 1, 26^n)],
	for i in read_list(file) do (
		a: charlist(sdowncase(printf(false, "~a", i))),
		for j in a do (
			if alphacharp(j) then (
				if length(s) = n-1 then (
					s: endcons(cint(j)-97, s*26),
					m: apply("+", s),
					L[m+1]: L[m+1] + 1,
					s: rest(s, 1)
				) else (
					s: endcons(cint(j)-97, s*26)
				)
			)
		)
	),
	return(L)
);
to_alphabet(N, n) := block([q: N, r, L: []],
	while(q # 0) do (
		[q, r]: divide(q, 26),
		L: cons(r, L)
	),
	L: append(makelist(0, i, 1, n), L),
	return(map(ascii, rest(L, length(L)-n) + 97))
);
next_char(L, c) := block([L1: [], L2: []],
	L1: makelist([L[i], i-1], i, 1, length(L)),
	L2: sort(L1, ordergreatp),
	for i in L2 do
		if quotient(i[2], 26) = cint(c)-97 then 
			print([to_alphabet(i[2], 2), i[1]])
);

correl(L1, L2) := block([x1, x2],
	x1: apply("+", L1)/length(L1),
	x2: apply("+", L2)/length(L2),
	sum((L1[i] - x1)*(L2[i] - x2), i, 1, length(L1))
		/sqrt(sum((L1[i] - x1)^2, i, 1, length(L1)))
		/sqrt(sum((L2[i] - x2)^2, i, 1, length(L2)))
);

solve_identity(EQ, X) := block([e, v, w, linsolvewarn: false],
	local(L),
	e: fullratsimp(lhs(EQ) - rhs(EQ)),
	e: expand(e * denom(e)),
	v: listofvars(e),
	if listp(X) then w: X else w: [X],
	if length(v) >= length(w) + 1 and length(v) = length(unique(append(v, w))) then (
		L[0]: [e],
		for j: 1 thru length(w) do (
			v: delete(w[j], v),
			L[j]: [],
			for k: 1 thru length(L[j-1]) do
				L[j]: append(L[j], makelist(coeff(L[j-1][k], w[j], i), i, 0, hipow(L[j-1][k], w[j])))
		),
		L[0]: unique(L[length(w)]),
		solve(L[0], v)
	)
);

elim_2nd([L]) := block([f, F, x, p, q, h, X, V, list: []],
	f: expand(L[1]),
	if length(L) = 2 then (
		x: L[2],
		h: hipow(f, x),
		if h >= 2 then (
			p: factor(coeff(f, x, h - 1)/(h * coeff(f, x, h))),
			F: collectterms(expand(subst(X - p, x, f)), X),
			return(subst(x + p, X, F))
		)
	) else (
		V: listofvars(f),
		for x in V do (
			h: hipow(f, x),
			if h >= 2 then (
				p: factor(coeff(f, x, h - 1)/(h * coeff(f, x, h))),
				F: collectterms(expand(subst(X - p, x, f)), X),
				list: cons(subst(x + p, X, F), list)
			)
		),
		return(list)
	)
);

f_pickup(F, V) := block([partswitch: true],
	if not(freeof(V, F)) then (
		for i: 1 while (true) do (
			G: part(F, i),
			if not(freeof(V, G)) then return(f_pickup(G, V))
			elseif G = end then return(F)
		)
	)
);

pickup(f, v) := block([f_pickup, L: []],
	f_pickup(F, V) := block([partswitch: true, d],
		if not(freeof(V, F)) then (
			d: 0,
			for i: 1 while (true) do (
				G: part(F, i),
				if not(freeof(V, G)) then (
					d: 1,
					f_pickup(G, V)
				) elseif G = end then (
					if (d = 0) then (
						L: cons(F, L)
					),
					return()
				)
			)
		)
	),
	f_pickup(f, v),
	L
);

unit_frac_decomp_g(r) := block([R: r, L: []],
	if (ratnump(R) and R > 0 and R < 1 ) then (
		for j while(R # 0) do (
			n: ceiling(1/R),
			R: R - 1/n,
			L: endcons(1/n, L)
		),
		return(L)
	) else (
		return(false)
	)
);