################################################# 11-14-97
## <SHAREFILE=algebra/crystal/crystal26.mpl >
## <DESCRIBE>
##                A Maple package which decomposes the tensor product of two
##                irreducible finite dimensional "fundamental" representations
##                of a simple Lie algebra into irreducible constituents using
##                "crystal" graphs.  This is an implementation of a theorem
##                of Kashiwara.  The package makes use of John Stembridge's
##                coxeter and weyl packages in the Maple share library.
##
##                AUTHOR: David Joyner, wdj@nadn.navy.mil
##                AUTHOR: Roland Martin, rem@sma.usna.navy.mil
##                AUTHOR: Michael Fourte, m982208@nadn.navy.mil
## <DESCRIBE>
## <UPDATE=R5>
#################################################

`crystal/reader` := proc(x) 
  eval(`crystal/`.x) 
end:

`crystal/weight_system` := proc(x,y)
local i,n,s,Cartan,HighestWeight,Nbase,Sbase,Weights;
  if y = `A1` then
    Weights:=[[1],[-1]];
    RETURN(Weights);
  else
    Sbase := coxeter[base](y);
    Nbase := nops(Sbase);
    for i to Nbase do  
      Cartan[i] := [seq(2*coxeter[iprod](s,op(i,Sbase))/coxeter[iprod](s,s),s = Sbase)];
    od;
    HighestWeight := [[seq(-2*coxeter[iprod](s,x)/coxeter[iprod](s,s),s = Sbase)]];
    crystal[addweight](HighestWeight,y)
  fi;
end:

`crystal/weight_coords` := proc(x)
local s,Sbase;
  Sbase := coxeter[base](R);
  [seq(-2*coxeter[iprod](s,x)/coxeter[iprod](s,s),s = Sbase)]
end:

`crystal/Stemweight_coords` := proc(x) 
local s,Sbase;
  Sbase := coxeter[base](R); 
  [seq(2*coxeter[iprod](s,x)/coxeter[iprod](s,s),s = Sbase)] 
end:

`crystal/addweight` := proc(x,y)
local a,b,i,j,k,s,Cartan,CartanVector,Nbase,Sbase,Weights,WeightVector;
  Sbase := coxeter[base](y);
  Nbase := nops(Sbase);
  for i to Nbase do  
    Cartan[i] := [seq(2*coxeter[iprod](s,op(i,Sbase))/coxeter[iprod](s,s),s = Sbase)];
  od;
  Weights := x;
  for i to nops(x) do
    WeightVector[i] := linalg[vector](op(i,x));
    for j to Nbase do
      if 0 < WeightVector[i][j] then
        for k to WeightVector[i][j] do
          CartanVector[j] := linalg[vector](Cartan[j]);
          a[j] := evalm(WeightVector[i]-k*CartanVector[j]);
          b[j] := convert(a[j],list);
          if not member(b[j],[op(Weights)]) then
            Weights := [op(Weights),b[j]];
          fi;
        od;
      fi;
    od;
  od;
  if nops(x) <> nops(Weights) then 
    crystal[addweight](Weights,y);
  else
    RETURN(Weights); 
  fi;
end:

`crystal/graphrep` := proc(x,y,z)
local a,i,j,k,s,Cartan,Nbase,Sbase;
  Sbase := coxeter[base](y);
  Nbase := nops(Sbase);
  for i to Nbase do  
    Cartan[i] := [seq(2*coxeter[iprod](s,op(i,Sbase))/coxeter[iprod](s,s),s = Sbase)];
  od;
  z := networks[new]();
  for i from 1 to nops(x) do
    networks[addvertex](v.i,weights = [op(i,x)],z);
  od;
  for i from 1 to nops(x) do
    for j from 1 to Nbase do
      for k from i to nops(x) do
        a := evalm(linalg[vector](op(i,x)) - linalg[vector](Cartan[j]));
        if linalg[equal](a,linalg[vector](op(k,x))) then         
          networks[addedge]([v.i,v.k],weights = j,z);
        fi;
      od;
    od;
  od;
  print(`Graph formed.`);
end:

`crystal/showgraph` := proc(x,y,z)
local i,j,k,label,labels,labelset,lines,lineset,points,pointset;
  labels := [];
  lines := [];
  points := [];
  for i from 1 to nops(networks[vertices](x)) do
    if member(v.i,networks[vertices](x)) then
      if i <= z and i >= y then
        points := [op(points),PLOT(POINTS([i,0],SYMBOL(CIRCLE)))];
      fi;
    fi;
  od;
  pointset := op(points);
  for i from 1 to nops(networks[vertices](x)) do
    for j from i to nops(networks[vertices](x)) do
      for k from 1 to nops(networks[edges](x)) do
        if i <= z then
          if networks[tail](e.k,x) = v.i then
            if y <= j then
              if networks[head](e.k,x) = v.j then
                if j-i=1 then
                  lines := [op(lines),PLOT(CURVES([[i,0],[j,0]]))];
                else
                  lines := [op(lines),PLOT(CURVES([[i,0],[(i+j)/2,1/(j-i)],[(i+j)/2,1/(j-i)],[j,0]]))];
                fi;
                label := convert(networks[eweight](e.k,x),string);
                if j-i=1 then
                  labels := [op(labels),PLOT(TEXT([(i+j)/2,0],`label`,ALIGNABOVE))];
                else
                  labels := [op(labels),PLOT(TEXT([(i+j)/2,1/(j-i)],`label`,ALIGNABOVE))];
                fi;
              fi;
            fi;
          fi;
        fi;
      od;
    od;
  od;
  lineset := op(lines);
  labelset := op(labels);
  plots[display]([pointset,lineset,labelset],axes = none,scaling = constrained);
end:

`crystal/graphprodrep` := proc(x,y,z)
local a,b,c,i,j,k,l,m,badxweights,badyweights,X;
  z := networks[new]();
  for i from 1 to nops(networks[vertices](x)) do
    a := linalg[vector](networks[vweight](v.i,x));
    for j from 1 to nops(networks[vertices](y)) do
      b := linalg[vector](networks[vweight](v.j,y));
      c := convert(evalm(a+b),list);
      networks[addvertex](v.i.X.j,weights = [c],z);
    od;
  od;
  for i from 1 to nops(networks[vertices](y)) do
    for j from 1 to nops(networks[vertices](x)) do
      for k from 1 to nops(networks[edges](x)) do
        if networks[head](e.k,x) = v.j then
          for m from 1 to j do
            badxweights := [];
            for l from 1 to nops(networks[edges](z)) do
              if networks[head](e.l,z) = v.j.X.i then
                badxweights := [op(badxweights),networks[eweight](e.l,z)];
              fi;
            od;
            if networks[tail](e.k,x) = v.m then
              if not member(networks[eweight](e.k,x),badxweights) then 
                a := networks[eweight](e.k,x);
                networks[addedge]([v.m.X.i,v.j.X.i],weights = [a],z);
              fi;
            fi;
          od;
        fi;
      od;
    od;    
    for j from i to nops(networks[vertices](y)) do
      for k from 1 to nops(networks[edges](y)) do
        if networks[tail](e.k,y) = v.i then
          for m from 1 to nops(networks[vertices](x)) do
            badyweights := [];
            for l from 1 to nops(networks[edges](z)) do
              if networks[tail](e.l,z) = v.m.X.i then
                badyweights := [op(badyweights),networks[eweight](e.l,z)];
              fi;
            od;
            if networks[head](e.k,y) = v.j then
              if not member(networks[eweight](e.k,y),badyweights) then
                a := networks[eweight](e.k,y);
                networks[addedge]([v.m.X.i,v.m.X.j],weights = [a],z);
              fi;
            fi;
          od;
        fi;
      od;
    od;        
  od;
  print(`Graph formed.`);
end:

`crystal/showgraphprod` := proc(u,w,x,y,z)
local a,b,i,j,k,l,m,n,cut,finish,index0,label,labels,labelset,
  lines,lineset,points,pointset,start,table0,X;
  a := 0;
  b := 0;
  labels := [];
  lines := [];
  points := [];
  for i from 1 to nops(networks[vertices](u)) do
    for j from 1 to nops(networks[vertices](u)) do
#print(`.8`,i,j,v.i.X.j,networks[vertices](u));
      if member(v.i.X.j,networks[vertices](u)) then
        if (i >= w and i <= x) then
          if (j >= y and j <= z) then
            points := [op(points),PLOT(POINTS([j,i],SYMBOL(CIRCLE)))];
#print(`.9`,i,j);
          fi;
        fi;
        if i > a then 
          a := i;
        fi;
        if j > b then
          b := j;
        fi;
      fi;
#print(`1`,i,j,a,b);
    od;
  od;
  pointset := op(points);
  for i from 1 to a do
    table0[i,1] := convert(i,string);
  od;
  for j from 1 to b do
    table0[1,j] := convert(j,string);
  od;    
  for m from 1 to nops(networks[edges](u)) do
    label := convert(networks[eweight](e.m,u),string);
    start := convert(networks[tail](e.m,u),string);
    finish := convert(networks[head](e.m,u),string);
    cut[1] := SearchText(X,start);
    cut[2] := SearchText(X,finish);
    index0[3] := substring(start,2..(cut[1]-1));
    index0[4] := substring(finish,2..(cut[2]-1));
    index0[1] := substring(start,(cut[1]+1)..2024);    
    index0[2] := substring(finish,(cut[2]+1)..2024);
    for n from 1 to b do
#print(`1.1`,index0[1], table0[1,n],n);
      if index0[1] = table0[1,n] then
        i := n;
      fi;
#print(`1.2`,index0[2], table0[1,n],n);
      if index0[2] = table0[1,n] then
        j := n;
      fi;
    od;
    for n from 1 to a do
#print(`1.3`,index0[3], table0[n,1],n);
      if index0[3] = table0[n,1] then
        k := n;
      fi;
#print(`1.4`,index0[4], table0[n,1],n);
      if index0[4] = table0[n,1] then
        l := n;
      fi;
    od;
    if k = l then
      if j-i = 1 then
        lines := [op(lines),PLOT(CURVES([[i,k],[j,k]]))];
        labels := [op(labels),PLOT(TEXT([i+1/2,k],`label`,ALIGNABOVE))];
#print(`2`,i,j,k);
      else
        lines := [op(lines),PLOT(CURVES([[i,k],[(i+j)/2,k+1/(j-i)],[(i+j)/2,k+1/(j-i)],[j,k]]))];
        labels := [op(labels),PLOT(TEXT([(i+j)/2,k+1/(j-i)],`label`,ALIGNABOVE))];
#print(`3`,i,j,k);
      fi;
    fi;
    if i = j then
      if l-k = 1 then
        lines := [op(lines),PLOT(CURVES([[i,k],[i,l]]))];
        labels := [op(labels),PLOT(TEXT([i,k+1/2],`label`,ALIGNLEFT))];
#print(`4`,i,k,l);
      else
        lines := [op(lines),PLOT(CURVES([[i,k],[i-1/(l-k),(k+l)/2],[i-1/(l-k),(k+l)/2],[i,l]]))];
        labels := [op(labels),PLOT(TEXT([i-1/(l-k),(k+l)/2],`label`,ALIGNLEFT))];
#print(`5`,i,k,l);
      fi;
    fi;
#print(`6`);
  od;
  lineset := op(lines);
  labelset := op(labels);
#print(`7`);
  plots[display]([pointset,lineset,labelset],axes = none,scaling = constrained);
end:

`crystal/linsubgraphrep` := proc(w,x,y,z)
local a,b,c,d,e,f,i,j,k,l,m,n,o,s,cut,index0,start,str,table0,tempgraph,
workgraph,Cartan,CartanMatrix,HighestWeight,IgusaWeight,IgusaWeights,Nbase,Order,Sbase,X;
  CartanMatrix := [];
  Sbase := coxeter[base](x);
  Nbase := nops(Sbase);
  for i to Nbase do  
    Cartan[i] := [seq(2*coxeter[iprod](s,op(i,Sbase))/coxeter[iprod](s,s),s = Sbase)];
    CartanMatrix := [op(CartanMatrix),Cartan[i]];
  od;
  str := convert(w,string);
  if searchtext(e,str) <> 0 then
    HighestWeight := [seq(-2*coxeter[iprod](s,w)/coxeter[iprod](s,s),s = Sbase)];
    for i from 1 to nops(networks[vertices](y)) do
      for j from 1 to nops(networks[vertices](y)) do
        if member(v.i.X.j,networks[vertices](y)) then        
          if networks[vweight](v.i.X.j,y) = HighestWeight then
            a := i;
            b := j;
          fi;
          break;          
        fi;
      od;
      break;
    od;
    tempgraph := networks[new]();
    for i from 1 to nops(networks[components](y)) do
      if member(v.a.X.b,op(i,networks[components](y))) then
        tempgraph := networks[induce](op(i,networks[components](y)),y);
      fi;
    od;
    start[1] := 1;
    start[2] := 1;
  else
    for i from 1 to trunc(sqrt(nops(networks[vertices](y)))) do
      table0[i] := convert(i,string);
    od;
    HighestWeight := networks[vweight](w,y);
    cut := SearchText(X,str);
    index0[1] := substring(str,2..(cut-1));
    index0[2] := substring(str,(cut+1)..2024);
    for i from 1 to trunc(sqrt(nops(networks[vertices](y)))) do
      if index0[1] = table0[i] then 
        a := i; 
      fi;
      if index0[2] = table0[i] then 
        b := i; 
      fi;
    od;   
    tempgraph := networks[new]();
    for i from 1 to nops(networks[components](y)) do
      if member(w,op(i,networks[components](y))) then
        tempgraph := networks[induce](op(i,networks[components](y)),y);
      fi;
    od;
    start[1] := a;
    start[2] := b;
  fi;
  d := 0;
  e := 0;
  f := [];
  for i from start[1] to nops(networks[vertices](tempgraph)) do
    for j from start[2] to nops(networks[vertices](tempgraph)) do
      if member(v.i.X.j,networks[vertices](tempgraph)) then
        a := linalg[vector](HighestWeight);
        b := linalg[vector](networks[vweight](v.i.X.j,y));
        c := convert(linalg[linsolve](linalg[matrix](CartanMatrix),evalm(a-b)),list);
        IgusaWeight[i][j] := sum(op(k,c),k=1..nops(c));  
        f := [op(f),IgusaWeight[i][j]];
        if i > d then
          d := i;
        fi;
        if j > e then
          e := j;
        fi;                
      fi;
    od;
  od;
  f:=sort(f);
  IgusaWeights := [op(1,f)];
  for i from 1 to nops(f) do
    if not member(op(i,f),IgusaWeights) then 
      IgusaWeights:=[op(IgusaWeights),op(i,f)];
    fi;
  od;
  f:=1;
  workgraph := networks[new]();
  z := networks[new]();
  for i from 0 to nops(IgusaWeights) do
    for j from 1 to d do
      for k from 1 to e do
        if member(v.j.X.k,networks[vertices](tempgraph)) then  
          if IgusaWeight[j][k]=op(i,IgusaWeights) then
            networks[addvertex](v.j.k.f,weights = [networks[vweight](v.j.X.k,y)],workgraph);
            networks[addvertex](v.f,weights = [networks[vweight](v.j.X.k,y)],z);
            f:=f+1;
          fi;
        fi;
      od;
    od;  
  od;
  for i from 1 to d do
    for j from i to d do
      for k from 1 to e do
        for l from k to e do
          for m from 1 to nops(networks[edges](y)) do
            if networks[tail](e.m,y) = v.i.X.k then
              if networks[head](e.m,y) = v.j.X.l then
                for n from 1 to f-1 do
                  for o from n to f-1 do
                    if member(v.i.k.n,networks[vertices](workgraph)) then
                      if member(v.j.l.o,networks[vertices](workgraph)) then
                        networks[addedge]([v.n,v.o],weights = [networks[eweight](e.m,y)],z);
                      fi;
                    fi;
                  od;
                od;
              fi;
            fi;
          od;
        od;
      od;
    od;
  od;
  print(`Graph formed.`);
end:

`crystal/dynkin` := proc(R)
  local i,n,arrows,arrowset,index0,label,labels,labelset,letter,lines,
   lineset,number,points,pointset;
  points := [];
  lines := [];
  labels := [];
  arrows := [];
  letter := substring(R,1..1);
  number := substring(R,2..8);
  for i from 1 to 2048 do
    index0 := convert(i,string);
    if index0 = number then
      n := i;
      break;
    fi;
  od;
  if letter = `A` then
    for i from 1 to n do
      points := [op(points),PLOT(POINTS([i/n,0],SYMBOL(CIRCLE)))];
      label := convert(i,string);
      labels := [op(labels),PLOT(TEXT([i/n,1/50],`label`,ALIGNABOVE))];
    od;
    pointset := op(points);
    labelset := op(labels);
    lines := [op(lines),PLOT(CURVES([[1/n,0],[1,0]]))];
    lineset := op(lines);
    plots[display]([pointset,lineset,labelset],title = `Dynkin diagram for A`.n,axes = none,scaling = constrained);
  elif letter = `B` then
    for i from 1 to n do
      points := [op(points),PLOT(POINTS([i/n,0],SYMBOL(CIRCLE)))];
      label := convert(i,string);
      labels := [op(labels),PLOT(TEXT([i/n,1/50],`label`,ALIGNABOVE))];
    od;
    pointset := op(points);
    labelset := op(labels);
    lines := [op(lines),PLOT(CURVES([[1/n,0],[1-1/n,0]]))];
    lines := [op(lines),PLOT(CURVES([[1-1/n,1/150],[1,1/150]]))];
    lines := [op(lines),PLOT(CURVES([[1-1/n,-1/150],[1,-1/150]]))];
    lineset := op(lines);
    arrows := [op(arrows),PLOT(CURVES([[1-1/(1.5*n),0],[1-1/(4*n),(1.25)/(8*n)]]))];
    arrows := [op(arrows),PLOT(CURVES([[1-1/(1.5*n),0],[1-1/(4*n),-(1.25)/(8*n)]]))];
    arrowset := op(arrows);
    plots[display]([pointset,labelset,lineset,arrowset],title = `Dynkin diagram for B`.n,axes = none,scaling = constrained);
  elif letter = `C` then
    for i from 1 to n do
      points := [op(points),PLOT(POINTS([i/n,0],SYMBOL(CIRCLE)))];
      label := convert(i,string);
      labels := [op(labels),PLOT(TEXT([i/n,1/50],`label`,ALIGNABOVE))];
    od;
    pointset := op(points);
    labelset := op(labels);
    lines := [op(lines),PLOT(CURVES([[1/n,0],[1-1/n,0]]))];
    lines := [op(lines),PLOT(CURVES([[1-1/n,1/150],[1,1/150]]))];
    lines := [op(lines),PLOT(CURVES([[1-1/n,-1/150],[1,-1/150]]))];
    lineset := op(lines);
    arrows := [op(arrows),PLOT(CURVES([[1-1/(4*n),0],[1-1/(1.5*n),(1.25)/(8*n)]]))];
    arrows := [op(arrows),PLOT(CURVES([[1-1/(4*n),0],[1-1/(1.5*n),-(1.25)/(8*n)]]))];
    arrowset := op(arrows);
    plots[display]([pointset,labelset,lineset,arrowset],title = `Dynkin diagram for C`.n,axes = none,scaling = constrained);
  elif letter = `D` then
    for i from 1 to n-2 do
      points := [op(points),PLOT(POINTS([i/n,0],SYMBOL(CIRCLE)))];
      label := convert(i,string);
      labels := [op(labels),PLOT(TEXT([i/n,1/50],`label`,ALIGNABOVE))];
    od;
    points := [op(points),PLOT(POINTS([1-1/n,1/(2*n)],SYMBOL(CIRCLE)))];
    label := convert(n-1,string);
    labels := [op(labels),PLOT(TEXT([1-1/n,1/(2*n)+1/50],`label`,ALIGNABOVE))];
    points := [op(points),PLOT(POINTS([1-1/n,-1/(2*n)],SYMBOL(CIRCLE)))];
    label := convert(n,string);
    labels := [op(labels),PLOT(TEXT([1-1/n,-1/(2*n)+1/50],`label`,ALIGNABOVE))];
    pointset := op(points);
    labelset := op(labels);
    lines := [op(lines),PLOT(CURVES([[1/n,0],[1-2/n,0]]))];
    lines := [op(lines),PLOT(CURVES([[1-2/n,0],[1-1/n,1/(2*n)]]))];
    lines := [op(lines),PLOT(CURVES([[1-2/n,0],[1-1/n,-1/(2*n)]]))];
    lineset := op(lines);
    plots[display]([pointset,labelset,lineset],title = `Dynkin diagram for D`.n,axes = none,scaling = constrained);
  elif letter = `E` then
    for i from 1 to n-1 do
      points := [op(points),PLOT(POINTS([i/n,0],SYMBOL(CIRCLE)))];
      label := convert(i,string);
      labels := [op(labels),PLOT(TEXT([i/n,-1/20],`label`,ALIGNABOVE))];
    od;
    points := [op(points),PLOT(POINTS([3/n,1/n],SYMBOL(CIRCLE)))];
    label := convert(n,string);
    labels := [op(labels),PLOT(TEXT([3/n,1/n+1/50],`label`,ALIGNABOVE))];
    pointset := op(points);
    labelset := op(labels);
    lines := [op(lines),PLOT(CURVES([[1/n,0],[1-1/n,0]]))];
    lines := [op(lines),PLOT(CURVES([[3/n,0],[3/n,1/n]]))];
    lineset := op(lines);
    plots[display]([pointset,labelset,lineset],title = `Dynkin diagram for E`.n,axes = none,scaling = constrained);
  elif letter = `F` then
    n := 4;
    for i from 1 to n do
      points := [op(points),PLOT(POINTS([i/n,0],SYMBOL(CIRCLE)))];
      label := convert(i,string);
      labels := [op(labels),PLOT(TEXT([i/n,-1/20],`label`,ALIGNABOVE))];
    od;
    pointset := op(points);
    labelset := op(labels);
    lines := [op(lines),PLOT(CURVES([[1/n,0],[2/n,0]]))];
    lines := [op(lines),PLOT(CURVES([[2/n,1/150],[3/n,1/150]]))];
    lines := [op(lines),PLOT(CURVES([[2/n,-1/150],[3/n,-1/150]]))];
    lines := [op(lines),PLOT(CURVES([[3/n,0],[1,0]]))];
    lineset := op(lines);
    arrows := [op(arrows),PLOT(CURVES([[2/n+1/(4*n),0],[2/n+1/(1.5*n),(1.25)/(8*n)]]))];
    arrows := [op(arrows),PLOT(CURVES([[2/n+1/(4*n),0],[2/n+1/(1.5*n),-(1.25)/(8*n)]]))];
    arrowset := op(arrows);
    plots[display]([pointset,labelset,lineset,arrowset],title = `Dynkin diagram for F`.n,axes = none,scaling = constrained);
  else
    n := 2;
    for i from 1 to n do
      points := [op(points),PLOT(POINTS([i/n,0],SYMBOL(CIRCLE)))];
      label := convert(i,string);
      labels := [op(labels),PLOT(TEXT([i/n,-1/20],`label`,ALIGNABOVE))];
    od;
    pointset := op(points);
    labelset := op(labels);
    lines := [op(lines),PLOT(CURVES([[1/n,1/175],[1,1/175]]))];
    lines := [op(lines),PLOT(CURVES([[1/n,0],[1,0]]))];
    lines := [op(lines),PLOT(CURVES([[1/n,-1/150],[1,-1/150]]))];
    lineset := op(lines);
    arrows := [op(arrows),PLOT(CURVES([[1/n+1/(4*n),0],[1/n+1/(1.5*n),(1.25)/(8*n)]]))];
    arrows := [op(arrows),PLOT(CURVES([[1/n+1/(4*n),0],[1/n+1/(1.5*n),-(1.25)/(8*n)]]))];
    arrowset := op(arrows);
    plots[display]([pointset,labelset,lineset,arrowset],title = `Dynkin diagram for G`.n,axes = none,scaling = constrained);
  fi;
end:

`crystal/branch` := proc(x,y,z)
local a,i,j;
  z := networks[new]();
  for i from 1 to nops(networks[vertices](y)) do
    networks[addvertex](v.i,z);
  od;
  for i from 1 to nops(networks[vertices](y)) do
    for j from 1 to nops(networks[edges](y)) do
      if networks[tail](e.j,y) = v.i then
        a := networks[eweight](e.j,y);
        if not member(a,x) then
          networks[addedge]([v.i,networks[head](e.j,y)],weights = a,z);
        fi;
      fi;
    od;
  od;
  print(`Branch graph formed.`);
end:

init_crystal := proc()
global crystal;
  crystal[weight_system] := `crystal/reader`(weight_system);
  crystal[weight_coords] := `crystal/reader`(weight_coords);
  crystal[Stemweight_coords] := `crystal/reader`(Stemweight_coords);
  crystal[addweight] := `crystal/reader`(addweight);
  crystal[graphrep] := `crystal/reader`(graphrep);
  crystal[graphprodrep] := `crystal/reader`(graphprodrep);
  crystal[showgraph] := `crystal/reader`(showgraph);
  crystal[showgraphprod] := `crystal/reader`(showgraphprod);
  crystal[linsubgraphrep] := `crystal/reader`(linsubgraphrep);
  crystal[dynkin] := `crystal/reader`(dynkin);
  crystal[branch] := `crystal/reader`(branch);
  print(`Crystal initialized.`);
end:

`help/crystal/text/graphrep` := TEXT(
`    `,
`FUNCTION : graphrep(L,R,G);`,
`    `,
`CALLING SEQUENCE : crystal[graphrep](L,R,G);`,
`    `,
`PARAMETERS :`,
` L is a list outputed by crystal[weight_system].`,
` R is a string denoting a simple Lie algebra from the list`,
` A1,A2,...,B2,B3,...,C3,C4,...,D4,D5,....,E6,E7,E8,F4,G2.`,
` G is a name denoting the output of the program.`,
`    `,
`SYNOPSIS :`,
` This program creates a crystal graph, denoted G, with`,
` weighted vertices [v1,v2,...,vn] and weighted, directed edges`,
` [e1,e2,...,em]. The crystal graph G may be viewed with the graph`,
` plotting program crystal[showgraph].`,
`    `,
`EXAMPLES :`,
` L:=crystal[weight_system](e1,A3);`,
` crystal[graphrep](L,A3,G);`,
`    `,
`SEE ALSO : crystal[weight_system], crystal[showgraph]`):

`help/crystal/text/graphprodrep` := TEXT(
`    `,
`FUNCTION : graphprodrep(G1,G2,G);`,
`    `,
`CALLING SEQUENCE : crystal[graphprodrep](G1,G2,G);`,
`    `,
`PARAMETERS :`,
` G1 and G2 are crystal graphs constructed by crystal[graphrep].`,
` G is a name denoting the output of the program.`, 
`    `,
`SYNOPSIS :`,
` This program computes the crystal graph product of the`,
` crystal graphs G1 and G2 in the sense of M. Kashiwara as`,
` detailed in his paper in Comm. Math. Physics, 1990. The crystal`,
` graph G may be viewed with the graph plotting program`,
` crystal[showgraphprod].`,
`    `,
`EXAMPLES :`,
` crystal[graphprodrep](G1,G2,G);`,
`    `,
`SEE  ALSO : crystal[graphrep], crystal[showgraphprod]`):

`help/crystal/text/showgraphprod` := TEXT(
`    `,
`FUNCTION : showgraphprod(G,x1,x2,y1,y2);`,
`    `,
`CALLING SEQUENCE : crystal[showgraphprod](G,x1,x2,y1,y2);`,
`    `,
`PARAMETERS :`,
` G is a graph created by crystal[graphprodrep]. (x1,y1) are the`,
` coordinates of the lower left hand vertex which the program`,
` will plot. (x2,y2) are the coordinates of the upper right`,
` hand vertex which the program will plot.`,
`    `,
`SYNOPSIS :`,
` This program will display in two dimensions a crystal graph`,
` created by crystal[graphprodrep].`,
`    `,
`EXAMPLES :`,
` L:=crystal[weight_system](e1,A3);`,
` crystal[graphrep](L,A3,G1);`,
` crystal[graphrep](L,A3,G2);`,
` crystal[graphprodrep](G1,G2,G);`,
` crystal[showgraphprod](G,1,4,1,4);`,
`    `,
`SEE ALSO : crystal[graphprodrep]`):

`help/crystal/text/weight_system` := TEXT(
`    `,
`FUNCTION : weight_system(v,R)`,
`    `,
`CALLING SEQUENCE : crystal[weight_system](v,R);`,
`    `,
`PARAMETERS :`,
` R is a string denoting a simple Lie algebra from the list of`,
` strings A1,A2,...,B2,B3,...,C3,C4,...,D4,D5,...,E6,E7,E8,F4 and`,
` G2. v is a positive rational linear combination of e1,...,en,`,
` where n denotes the rank of R and the ei denote the standard`,
` basis for R.`,
`    `,
`SYNOPSIS :`,
` This program will return the list of weights corresponding`,
` to the representation of R with highest weight v. However,`,
` crystal[weight_system] does not compute this list of weights with`,
` multiplicity. Note that one can use the crystal package to`,
` construct the crystal graph of most representations by building up`,
` from the appropriate fundamental representations which generally`,
` have multiplicity free weight systems.`,
`    `,
`EXAMPLES :`,
` crystal[weight_system](e1,A3);`,
`    `,
`SEE ALSO : the TeX documentation crystal.tex`):

`help/crystal/text/showgraph` :=TEXT(
`    `,
`FUNCTION : showgraph(G,x1,x2);`,
`    `,
`CALLING SEQUENCE : crystal[showgraph](G,x1,x2);`,
`    `,
`PARAMETERS :`,
` G is a graph created by crystal[graphrep]. x1 is the`,
` coordinate of the left hand vertex from which the `,
` program will begin plotting. x2 is the coordinate of the`,
` right hand vertex where the program will terminate plotting.`,
`    `,
`SYNOPSIS :`,
` This program will display in a single dimension a crystal`,
` graph created by crystal[graphrep].`,
`    `,
`EXAMPLES :`,
` L:=crystal[weight_system](e1,A3);`,
` crystal[graphrep](L,A3,G);`,
` crystal[showgraph](G,1,4);`,
`    `,
`SEE ALSO : crystal[graphrep]`):

`help/crystal/text/linsubgraphrep` := TEXT(
`    `,
`FUNCTION : linsubgraphrep(v,R,G1,G);`,
`    `,
`CALLING SEQUENCE : crystal[linsubgraphrep](v,R,G1,G);`,
`    `,
`PARAMETERS :`,
` R is a string denoting a simple Lie algebra from the list`,
` A1,A2,...,B2,B3,...,C3,C4,...,D4,D5,....,E6,E7,E8,F4,G2.`,
` v is a positive rational linear combination of e1,...,en`,
` indicating the highest weight of G1, n denoting the rank`,
` of R and the ei denoting the standard basis for R, or an`,
` explicit vertex of G1. G1 is a graph created by`,
` crystal[graphprodrep]. G is a name denoting the output of`,
` the program.`,
`    `,
`SYNOPSIS :`,
` This program creates the crystal graph G from the connected`,
` component of the crystal graph G1 with vertex of weight v`,
` if v is the highest weight or creates a crystal graph G from`,
` the connected component of the crystal graph G1 with vertex v`,
` of weight w and those vertices with weights less than w.`,
` The crystal graph G may be viewed with the graph plotting`,
` program crystal[showgraph].`,
`    `,
`EXAMPLES :`,
` L:=crystal[weight_system](e1,A3);`,
` crystal[graphrep](L,A3,G1);`,
` crystal[graphrep](L,A3,G2);`,
` crystal[graphprodrep](G1,G2,G3);`,
` crystal[linsubgraphrep](2*e1,A3,G3,G);`,
`    `,
`SEE ALSO : crystal[graphprodrep], crystal[showgraph]`):

`help/text/crystal`:=TEXT(
`    `,
`Help file for the crystal package.`,
`    `,
`CALLING SEQUENCE : crystal[<function>](<parameters>);`,
`    `,
`SYNOPSIS :`,
` This package will help decompose tensor products of `,
` irreducible finite dimensional representations of simple Lie`,
` algebras. The irreducible constituents of the decomposition `,
` correspond, by a theorem of Kashiwara, to the connected`,
` components of a crystal graph. This package uses the`,
` MAPLE shareware packages weyl and coxeter.`,
`    `,
` In MAPLEV3, type:`,
`  read 'crystal25.mpl';`,
`  init_crystal();`,
`  with(share);`,
`  readshare(coxeter,coxeter);`,
`  readshare(weyl,coxeter);  .`,
`    `,
` Help is available for the following functions:`,
`  graphprod, graphrep, weight_system, showgraph,`,             
`  linsubgraphrep and showgraphprod.`,
`    `,
`EXAMPLE :`,
` ?crystal[weight_system];`,
`    `,
`SEE ALSO : the TeX documentation crystal.tex.`):

`help/crystal/text/dynkin` := TEXT(
`    `,
`FUNCTION : dynkin(R);`,
`    `,
`CALLING SEQUENCE : crystal[dynkin](R);`,
`    `,
`PARAMETERS :`,
` R is a string denoting a simple Lie algebra from the list`,
` A1,A2,...,B2,B3,...,C3,C4,...,D4,D5,....,E6,E7,E8,F4,G2.`,
`    `,
`SYNOPSIS :`,
` This program plots the Dynkin diagram associated to R.`,
`    `,
`EXAMPLES :`,
` crystal[dynkin](A3);`):

`help/crystal/text/branch` := TEXT(
`    `,
`FUNCTION : branch(L,G,H);`,
`    `,
`CALLING SEQUENCE : crystal[branch](L,G,H);`,
`    `,
`PARAMETERS :`,
` L is a list of weights. G is a crystal graph.`,
` H is a name denoting the output of the program.`,
`    `,
`SYNOPSIS :`,
` This program creates a crystal graph, denoted H,`,
` associated to a branching rule for the simple Lie`,
` algebra associated to G.`,
`    `,
`EXAMPLES :`,
` L:=crystal[weight_system](e1,A2);`,
` crystal[graphrep](L,A2,G);`,
` crystal[branch]([2],G,H);`,
` crystal[showgraph](H,1,3);`,
`    `,
`SEE ALSO : crystal[weight_system], crystal[graphrep], crystal[showgraph]`):