###################################################################3


Help:=proc():

if nargs=0 then

print(` Version of May 3 , 2003 `):
print(`   `):
print(` Welcome to LLL, a Maple package`):
print(`  `):
print(` written by Mohamud Mohammed as a project for the course`):
print(`  `):
print(` Math 583, Combinatorics II, by Prof. Doron Zeilberger`):
print(`    `):
print(`Please send all bugs to mohamudm@math.rutgers.edu`):
print(`  `):
print(`This is an implementation of the Integer lattice `):
print(`   `):
print(` reduction algorithm,( LLL algorithm) by `):
print(`     `):
print(` A.K. Lenstra, H.W. Lenstra, L. Lovasz with some applications`):
print(`    `):
print(` From the original paper: Factoring polynomials`):
print(`    `):
print(` with Rational Coefficients,`):
print(`    `):
print(` Math. Ann., 261, 515-534(1982)`): 
print(`    `):
print(` The default inner product is that of Euclid: If one wants`):
print(`     `):
print(` can change the procedure VP and use with different spaces`):
print(`    `):
print( ` One instance is on the function space with`):
print((f,g)=Int(f(x)*g(x),x=a..b)); 

print(`  `):

print(`Contains the following procedures:`):
print(`   `):
print(`VP(inner product), LLL, SP:(mult. by scalar.),`):
print(`  `):
print(`  VD(vector difference), VS (vector sum)`):
print(`  `):
print(` FirstCond, SecondCond, GrScOr, IntRel, MinPoly, LLLApp`):

elif nargs=1 then

if args[1]=VP then
print(`VP(v,u) : given two vectors of the same dimensions`):
print(`u and v, returns their Euclidean inner product`):
print(`Try InnPr([1,2],[5,6]): it should return: 17`):
print(`  `):
fi:

if args[1]=InnProduct then
print(`InnProduct(u,v) : Finds the inner product of u and v`):
print(`according to the rule specified.The default here is`):
print(`Euclidean but one can also change in a trivial`):
print(`way to suit ones need.`):
fi:

if args[1]=GrScOr then
print(`GrScOr(basis) : given independent vectors, basis, it `):
print(`computes the orthogonal basis using Gram Schmidt`):
print(` orthogonalization and parameters associated to it.`):
print(` If you want to out put the orthogonal basis do : `):
print(`GrScOr(basis)[3]`):
print(`  `):
fi:

if args[1]=LLLApp then
print(`LLLApp(basis,C) : The same as LLL(basis,C) but we don't use`):
print(` Gauss elimination for independence here, since in most `):
print(` applications we have real numbers and Gauss ellim. `):
print(` is over rational field`):
print(`   `):
fi:

if args[1]=IntRel then
print(`IntRel(basis,C) : given a vector of real `):
print(`numbers and a constant C, IntRel finds a vector`):
print(`of integers which annihilates basis: a^T*b=0 `):
print(`Try IntRel([1,Pi^2,Pi^4,Pi^6,Pi^8,V],10^15)`):
print(`where`):
print( V = Int(sqrt(x)*ln(x)^5/(1-x)^5,x=0..infinity)):
print(`should output :[0, 120, 140, -15, 0, 24,0.1*10^(-11)]  `):
print(`  `): 
print(` more interesting example , try:`):
print(`IntRel([a1,a2],10^15] , where`):

print(a1=Sum(1/n^3,n=1..infinity));
print(a2= Sum((-1)^(k+1)/k^3/binomial(2*k,k) ,k=1..infinity)); 
print(` Output :[-2, 5, -.10e-213], [Zeta(3), 1/2*hypergeom([1, 1, 1, 1],`):
print(`[2, 2, 3/2],-1/4)]`):
print(`  `): 
print(` cool !  Apery's formula.`):

print(`  `):
print(`Even More ! Try :  IntRel([c,d],10^15): where `):

print(c= Sum((-1)^n*n!^10*(205*n^2+250*n+77)/(2*n+1)!^5,n=0..infinity));

print(d= Sum(1/n^3,n=1..infinity));

print(`to get :  [1,-64,0], [ c,d]`):

print(`c and d as above.  `):

print(` Amdberhan and Zeilberger's accelerating formula for zeta(3) !`)

fi:

if args[1]=MinPoly then
print(`If alpha is a real number, then `):
print(`alpha is an algebraic number if for some r, the vector `):
print(`[1,alpha,alph^2,...,alpha^r] has an integer relation`):
print(`MinPoly(vect): given vect in the above form finds the`):
print(`minimal polynomial for alpha.`): 
print(`To write down your min. poly., set the vector product`):
print(`of the first vec. less the last component with the second to`):
print(` zero and replace alpha by x and write everything interms of x`):
print(`Try MinPoly([1,sqrt(2),2]): it should return: `):
print(`[-2,0,1,0],[1,sqrt(2),2]`):
print(`  `):
fi:


if args[1]=LLL then
 print(`LLL(list) : given a list of independent vectors,(b[1],..,b[n])`):
 print(`returns the reduced linearly independent base vectors : `):
 print(`(b[1]*,..,b[n]*)`):
 print(`Try:  LLL([1,2,3],[2,1,6],[1,5,7]), it should output :`): 
 print(`[[-1, 1, 1], [2, 1, 2], [-1, -2, 1]]`):
print(`  `):
fi:

fi:

end:


with(linalg):

VP:= proc(a,b)
local i:
if nops(a)<>nops(b) then
  RETURN(`Error: the two vectors should have the same dimension`):
fi:

convert([seq(a[i]*b[i],i=1..nops(a))],`+`):

end: 

SP:=proc(c,w)
local i:
[seq(c*w[i],i=1..nops(w))]:
end:


VD:=proc(v,u)
local i:
[seq(v[i]-u[i],i=1..nops(v))]:
end:

VS:=proc(v,u)
local i:
[seq(v[i]+u[i],i=1..nops(v))]:
end:


GrScOr:= proc(basis)
 local i,j,k,mu,B,new_basis,dim:

dim:=nops(basis):

for i from 1 to  dim do  
  new_basis[i]:=basis[i]:    
 
  for j from 1 to i-1 do
    mu[i,j]:= VP(basis[i],new_basis[j])/B[j]: 
    new_basis[i]:=VD(new_basis[i], SP(mu[i,j],new_basis[j])):
  od:
       
  B[i]:=VP(new_basis[i],new_basis[i]):
od:


for j from 1 to dim do
 new_basis[j]:=basis[j]:
od:

[seq([seq(mu[i,j],j=1..dim)],i=1..dim)],
[seq(B[i],i=1..dim)],
[seq(new_basis[j],j=1..dim)]:

end:


FirstCond:= proc(k,l,A,basis,dim)
 local mu, new_basis,q,i,j: 
 
 mu:=A:
 new_basis:=basis:
  
if abs(mu[k,l]) > 1/2 then
    q:= floor(mu[k,l]+1/2):
    new_basis[k]:= VD(new_basis[k],SP(q,new_basis[l])):

    for j from 1 to l-1 do
      mu[k,j]:=mu[k,j]-q*mu[l,j]:
    od:

    mu[k,l]:=mu[k,l]-q:
 fi:

[seq([seq(mu[i,j],j=1..dim)],i=1..dim)],
[seq(new_basis[j],j=1..dim)]:

end:


SecondCond:= proc(k,A,b,basis,dim)
    local new_basis,mu,B,u_temp,u_perm,u_prime,
          B_temp,b_temp,i,j:
   
    B:=b:
    new_basis:=basis:
    mu:=A:

    u_temp:=mu[k,k-1]:
    B_temp:=B[k]+u_temp^2*B[k-1]:
    mu[k,k-1]:=u_temp*B[k-1]/B_temp:
    B[k]:=B[k-1]*B[k]/B_temp:
    B[k-1]:=B_temp:


    b_temp:=new_basis[k]:
    new_basis[k]:= new_basis[k-1]:     
    new_basis[k-1]:=b_temp:       


   for j from 1 to  k-2  do
    u_perm:=mu[k,j]:
    mu[k,j]:=mu[k-1,j]:
    mu[k-1,j]:=u_perm: 
   od:

    

for i from k+1 to dim do
    u_prime:=mu[i,k]:
    mu[i,k]:=mu[i,k-1] - u_temp*u_prime:
    mu[i,k-1]:=u_prime + mu[k,k-1]*mu[i,k]:
od:


[seq([seq(mu[i,j],j=1..dim)],i=1..dim)],
[seq(B[i],i=1..dim)],
[seq(new_basis[j],j=1..dim)]:

end:


RevMM:= proc(basis)
 local i,j:
[seq(basis[nops(basis)-i],i=0..nops(basis)-1)]:

end: 


LLL:= proc(basis)
 local aa,l,j,new_basis,mu,n,dim,q,a,B,k,b: 

dim:=nops(basis):
a:=GrScOr(basis):

ffgausselim(matrix(basis),'r'):

if r < dim then
 RETURN(`ERROR : the vectors are linearly dependent`):
fi:


mu:=a[1]:
B:=a[2]:
new_basis:=a[3]:

k:=2:

while k <= dim do

a:=FirstCond(k,k-1,mu,new_basis,dim):     # First condition

mu:=a[1]:
new_basis:=a[2]:

if B[k] < (3/4-mu[k,k-1]^2)*B[k-1] then   # Second condition

b:=SecondCond(k,mu,B,new_basis,dim):

mu:=b[1]:
B:=b[2]:
new_basis:=b[3]:

if k > 2 then
  k:=k-1:
fi:

next:

fi:


for l from k-2 by -1 to 1 do

   aa:= FirstCond(k,l,mu,new_basis,dim):
   mu:= aa[1]:
   new_basis:=aa[2]:

od:

k:=k+1:

od: #outer while loop.

new_basis:

end:      


LLLApp:= proc(basis)
 local aa,l,j,new_basis,mu,n,dim,q,a,B,k,b: 

dim:=nops(basis):
a:=GrScOr(basis):


mu:=a[1]:
B:=a[2]:
new_basis:=a[3]:

k:=2:

while k <= dim do

a:=FirstCond(k,k-1,mu,new_basis,dim):     # First condition

mu:=a[1]:
new_basis:=a[2]:

if B[k] < (3/4-mu[k,k-1]^2)*B[k-1] then   # Second condition

b:=SecondCond(k,mu,B,new_basis,dim):

mu:=b[1]:
B:=b[2]:
new_basis:=b[3]:

if k > 2 then
  k:=k-1:
fi:

next:

fi:


for l from k-2 by -1 to 1 do

   aa:= FirstCond(k,l,mu,new_basis,dim):
   mu:= aa[1]:
   new_basis:=aa[2]:

od:

k:=k+1:

od: #outer while loop.

new_basis:

end:      


IntRel:=proc(basis1,C)
  local i,index,j,dim,new_basis,basis,a:

dim:=nops(basis1):
Digits:=50:
basis:= basis1:

for i from 1 to dim do
    if not(type(basis[i], rational)) then
     basis[i] := evalf(basis[i]):
    fi:
od:

new_basis:= []:

for i from 1 to dim do
  new_basis:=[op(new_basis),[seq(0,j=1..i-1),1,seq(0, j=1..dim-i),
             C*basis[i]]]:
od:

a:=LLLApp(new_basis):

index:=1:

for i from 2 to nops(dim) do
 if  a[i][nops(a[1])] < a[index][nops(a[1])] then
  index:=i:
 fi:
od:

a[index],basis1:
a;
end:

MinPoly:=proc(basis): IntRel(basis,10^(15)): end: