File: [local] / OpenXM_contrib / PHC / Ada / Schubert / ts_straighten.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:32 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with text_io,integer_io; use text_io,integer_io;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Brackets,Brackets_io; use Brackets,Brackets_io;
with Bracket_Monomials; use Bracket_Monomials;
with Bracket_Polynomials; use Bracket_Polynomials;
with Bracket_Polynomials_io; use Bracket_Polynomials_io;
with Straightening_Syzygies; use Straightening_Syzygies;
procedure ts_straighten is
-- DESCRIPTION :
-- Enumerates all brackets in the Pluecker embedding and generates
-- the van der Waerden syzygies.
function Number_of_Brackets ( n,d : natural ) return natural is
-- DESCIPTION :
-- Returns the number of brackets of d entries chosen from n numbers.
a,b : natural;
begin
a := 1;
for i in d+1..n loop
a := a*i;
end loop;
b := 1;
for i in 1..n-d loop
b := b*i;
end loop;
return a/b;
end Number_of_Brackets;
function Number_of_Linear_Equations ( n,d : natural ) return natural is
-- DESCRIPTION :
-- Returns the number of linear equations you need to determine a
-- d-plane in C^n completely.
begin
return (n-d)*d;
end Number_of_Linear_Equations;
function Number_of_Zero_Brackets ( n,d : natural ) return natural is
-- DESCRIPTION :
-- Returns the maximal number of brackets that can be set to zero.
begin
return (Number_of_Brackets(n,d) - Number_of_Linear_Equations(n,d) - 1);
end Number_of_Zero_Brackets;
function Create ( b1,b2 : Bracket ) return Bracket_Term is
bm : Bracket_Monomial := Create(b1);
bt : Bracket_Term;
begin
Multiply(bm,b2);
bt.coeff := Create(1.0);
bt.monom := bm;
return bt;
end Create;
procedure Enumerate ( n,d,k,start : in natural; accu : in out Bracket;
cnt : in out natural ) is
-- DESCRIPTION :
-- Lexicographic enumeration of all brackets.
begin
if k > d
then -- put(accu); new_line;
cnt := cnt + 1;
else for l in start..n loop
accu(k) := l;
Enumerate(n,d,k+1,l+1,accu,cnt);
end loop;
end if;
end Enumerate;
procedure Enumerate2 ( n,d,k,start : in natural; b : in Bracket;
accu : in out Bracket;
cntstd,cntnonstd : in out natural;
nonstd : in out Bracket_Polynomial ) is
-- DESCRIPTION :
-- Lexicographic enumeration of all brackets, with b < accu and with
-- a test whether the pair b*accu forms a Standard monomial.
-- ON ENTRY :
-- n total number of indices to choose from;
-- d number of indices in the brackets;
-- k current entry in accu that is to be filled;
-- start indices will be taken in between start..n;
-- b previously enumerated bracket;
-- accu accumulating parameter, filled in upto (k-1)th entry;
-- cntstd current number of quadratic standard monomials;
-- cntnonstd current number of quadratic nonstandard monomials;
-- nonstd current polynomial of quadratic nonstandard monomials.
-- ON RETURN :
-- accu updated accumulating parameter, accu(k) is filled in;
-- cnstd updated number of quadratic standard monomials;
-- cntnonstd updated number of quadratic nonstandard monomials;
-- nonstd updated polynomial of quadratic nonstandard monomials.
s : natural;
begin
if k > d
then if b < accu
then -- put(b); put(accu);
s := Brackets.Is_Standard(b,accu);
if s = 0
then -- put_line(" is a Standard monomial.");
cntstd := cntstd + 1;
else -- put(" is not a Standard monomial with s = ");
-- put(s,1); new_line;
cntnonstd := cntnonstd + 1;
Add(nonstd,Create(b,accu));
end if;
end if;
else for l in start..n loop
accu(k) := l;
Enumerate2(n,d,k+1,l+1,b,accu,cntstd,cntnonstd,nonstd);
end loop;
end if;
end Enumerate2;
procedure Enumerate1 ( n,d,k,start : natural; acc1,acc2 : in out Bracket;
cntstd,cntnonstd : in out natural;
nonstd : in out Bracket_Polynomial ) is
-- DESCRIPTION :
-- Lexicographic enumeration of all brackets, with acc1 < acc2 and with
-- a test whether the pair acc1*acc2 forms a Standard monomial.
-- Counts the standard and nonstandard monomials and adds every
-- nonStandard monomial to the polynomial nonstd.
begin
if k > d
then Enumerate2(n,d,1,acc1(1),acc1,acc2,cntstd,cntnonstd,nonstd);
else for l in start..n loop
acc1(k) := l;
Enumerate1(n,d,k+1,l+1,acc1,acc2,cntstd,cntnonstd,nonstd);
end loop;
end if;
end Enumerate1;
procedure Enumerate_Pairs ( n,d : in natural;
nonstd : in out Bracket_Polynomial ) is
-- DESCRIPTION :
-- Enumerates all pairs (b1,b2), with b1 < b2 and checks whether
-- the corresponding bracket monomial b1*b2 is standard or not.
b1,b2 : Bracket(1..d);
cntstd,cntnonstd : natural := 0;
begin
Enumerate1(n,d,1,1,b1,b2,cntstd,cntnonstd,nonstd);
put("#Standard bi-monomials : "); put(cntstd,1); new_line;
put("#nonStandard bi-monomials : "); put(cntnonstd,1); new_line;
end Enumerate_Pairs;
procedure Enumerate_Brackets ( n,d : in natural ) is
acc : Bracket(1..d);
cnt : natural := 0;
nonstd : Bracket_Polynomial;
begin
Enumerate(n,d,1,1,acc,cnt);
put("#brackets : "); put(cnt,1); new_line;
Enumerate_Pairs(n,d,nonstd);
put_line("The polynomial of all quadratic nonstandard monomials :");
put(nonstd);
end Enumerate_Brackets;
procedure Read_Bracket ( b : in out Bracket ) is
d : constant natural := b'last;
sig : integer;
begin
put("Give "); put(d,1); put(" row indices : "); get(b,sig);
put("The sorted bracket : ");
if sig > 0
then put("+");
else put("-");
end if;
put(b);
if Is_Zero(b)
then put_line(" is zero.");
else put_line(" is different from zero.");
end if;
end Read_Bracket;
procedure Test_Straighten ( n,d : in natural ) is
b,bb : Bracket(1..d);
s : natural;
ans : character;
bp : Bracket_Polynomial;
begin
loop
Read_Bracket(b);
Read_Bracket(bb);
if Is_Equal(b,bb)
then put_line("Both brackets are equal.");
end if;
if b < bb
then put(b); put(" < "); put(bb);
s := Brackets.Is_Standard(b,bb);
put(" and "); put(b); put(bb);
bp := Straightening_Syzygy(b,bb);
else put(b); put(" >= "); put(bb);
s := Brackets.Is_Standard(bb,b);
put(" and "); put(bb); put(b);
bp := Straightening_Syzygy(bb,b);
end if;
if s = 0
then put_line(" is a Standard monomial.");
else put(" is not a Standard monomial, with s = "); put(s,1); new_line;
end if;
put_line("The straightening syzygy : "); put(bp);
put("Do you want to test more pairs of brackets ? (y/n) "); get(ans);
exit when ans /= 'y';
end loop;
end Test_Straighten;
procedure Test_Laplace ( n,d : natural ) is
p : Bracket_Polynomial := Laplace_Expansion(n,d);
begin
put("The Laplace expansion of "); put(n,1); put("*"); put(n,1);
put("-determinant as product of "); put(d,1); put("- and ");
put(n-d,1); put_line("-blocks : ");
put(p);
end Test_Laplace;
function Decompose ( n,d : natural; p : Bracket_Polynomial )
return Bracket_Polynomial is
-- DESCRIPTION :
-- Decomposes the ideal of square-free nonStandard monomials.
res : Bracket_Polynomial;
np : constant natural := Number_of_Monomials(p);
brackmat : array(1..np,1..2) of Link_to_Bracket;
procedure Initialize is
-- DESCRIPTION :
-- Initializes the matrix of the brackets that occur in p.
cnt_term : natural := 0;
procedure List_Term ( t : in Bracket_Term; continue : out boolean ) is
cnt_mon : natural := 0;
procedure Store_Bracket ( b : in Bracket; continue : out boolean ) is
begin
cnt_mon := cnt_mon + 1;
brackmat(cnt_term,cnt_mon) := new Bracket'(b);
continue := true;
end Store_Bracket;
procedure Store_Brackets is
new Bracket_Monomials.Enumerate_Brackets(Store_Bracket);
begin
cnt_term := cnt_term+1;
Store_Brackets(t.monom);
continue := true;
end List_Term;
procedure List_Terms is new Enumerate_Terms(List_Term);
begin
List_Terms(p);
end Initialize;
procedure Write is
-- DESCRIPTION :
-- Writes the matrix of brackets.
begin
for i in 1..np loop
for j in 1..2 loop
put(" ");
put("b("); put(i,1); put(","); put(j,1); put(") :");
put(brackmat(i,j).all);
end loop;
new_line;
end loop;
end Write;
function Occured_Yet ( k : natural; b : Bracket ) return boolean is
-- DESCRIPTION :
-- Returns true if the bracket b occurs in one of the rows <= k.
begin
for i in 1..k loop
for j in 1..2 loop
if Is_Equal(brackmat(i,j).all,b)
then return true;
end if;
end loop;
end loop;
return false;
end Occured_Yet;
procedure Solve is
-- DESCRIPTION :
-- Solves the initial ideal of quadratic nonStandard monomials.
accu : Bracket_Monomial;
nzrb : constant natural := Number_of_Zero_Brackets(n,d);
procedure Solve ( k : in natural; acc : in Bracket_Monomial ) is
begin
if k > np
then Add(res,Create(acc));
else if Divisible(acc,brackmat(k,1).all)
or else Divisible(acc,brackmat(k,2).all)
then Solve(k+1,acc);
else if Number_of_Brackets(acc) < nzrb
then for i in 1..2 loop
declare
newacc : Bracket_Monomial;
begin
Copy(acc,newacc);
Multiply(newacc,brackmat(k,i).all);
Solve(k+1,newacc);
end;
end loop;
end if;
end if;
end if;
end Solve;
begin
Solve(1,accu);
end Solve;
begin
Initialize;
-- Write;
Solve;
return res;
end Decompose;
procedure Enumerate_Straightening_Syzygies ( n,d : in natural ) is
-- DESCRIPTION :
-- Sets up of the straightening syzygies for all nonStandard
-- quadratic monomials.
nonstd,decom : Bracket_Polynomial;
procedure Write_Syzygy ( s : in Bracket_Polynomial; cont : out boolean ) is
begin
put(s); cont := true;
end Write_Syzygy;
procedure Write_Syzygies is new Enumerate_Syzygies(Write_Syzygy);
begin
put("(n,d) = "); put("("); put(n,1); put(","); put(d,1); put(")");
put(" #Brackets : "); put(Number_of_Brackets(n,d),1);
put(" #Linear Equations : ");
put(Number_of_LInear_Equations(n,d),1); new_line;
put_line("Type of linear equation :");
put(Laplace_Expansion(n,n-d));
nonstd := nonStandard_Monomials(n,d);
put_line("The polynomial with all quadratic nonStandard monomials :");
put(nonstd);
put_line("All quadratic straightening syzygies : ");
Write_Syzygies(nonstd);
put("#Zero-Brackets : "); put(Number_of_Zero_Brackets(n,d),1); new_line;
decom := Decompose(n,d,nonstd);
put_line("Decomposition of the ideal of nonStandard monomials :");
put(decom);
put("Number of solutions : "); put(Number_of_Monomials(decom),1);
end Enumerate_Straightening_Syzygies;
procedure Main is
n,d : natural;
ans : character;
begin
new_line;
put_line("Interactive testing of straightening algorithm.");
loop
new_line;
put_line("Choose one of the following : ");
put_line(" 0. Exit this program.");
put_line(" 1. Enumerate all brackets.");
put_line(" 2. Test the straightening algorithm.");
put_line(" 3. Test Laplace expansion.");
put_line(" 4. Enumerate straightening syzygies.");
put("Make your choice (0,1,2,3,or 4) : "); get(ans);
exit when ans = '0';
new_line;
put("Give the number of entries in bracket : "); get(d);
put("Give the number of elements to choose from : "); get(n);
case ans is
when '1' => Enumerate_Brackets(n,d);
when '2' => Test_Straighten(n,d);
when '3' => Test_Laplace(n,d);
when '4' => Enumerate_Straightening_Syzygies(n,d);
when others => put_line("Bad answer.");
end case;
end loop;
end Main;
begin
Main;
end ts_straighten;