[BACK]Return to ts_straighten.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Schubert

Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_straighten.adb, Revision 1.1

1.1     ! maekawa     1: with text_io,integer_io;                 use text_io,integer_io;
        !             2: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
        !             3: with Brackets,Brackets_io;               use Brackets,Brackets_io;
        !             4: with Bracket_Monomials;                  use Bracket_Monomials;
        !             5: with Bracket_Polynomials;                use Bracket_Polynomials;
        !             6: with Bracket_Polynomials_io;             use Bracket_Polynomials_io;
        !             7: with Straightening_Syzygies;             use Straightening_Syzygies;
        !             8:
        !             9: procedure ts_straighten is
        !            10:
        !            11: -- DESCRIPTION :
        !            12: --   Enumerates all brackets in the Pluecker embedding and generates
        !            13: --   the van der Waerden syzygies.
        !            14:
        !            15:   function Number_of_Brackets ( n,d : natural ) return natural is
        !            16:
        !            17:   -- DESCIPTION :
        !            18:   --   Returns the number of brackets of d entries chosen from n numbers.
        !            19:
        !            20:     a,b : natural;
        !            21:
        !            22:   begin
        !            23:     a := 1;
        !            24:     for i in d+1..n loop
        !            25:       a := a*i;
        !            26:     end loop;
        !            27:     b := 1;
        !            28:     for i in 1..n-d loop
        !            29:       b := b*i;
        !            30:     end loop;
        !            31:     return a/b;
        !            32:   end Number_of_Brackets;
        !            33:
        !            34:   function Number_of_Linear_Equations ( n,d : natural ) return natural is
        !            35:
        !            36:   -- DESCRIPTION :
        !            37:   --   Returns the number of linear equations you need to determine a
        !            38:   --   d-plane in C^n completely.
        !            39:
        !            40:   begin
        !            41:     return (n-d)*d;
        !            42:   end Number_of_Linear_Equations;
        !            43:
        !            44:   function Number_of_Zero_Brackets ( n,d : natural ) return natural is
        !            45:
        !            46:   -- DESCRIPTION :
        !            47:   --   Returns the maximal number of brackets that can be set to zero.
        !            48:
        !            49:   begin
        !            50:     return (Number_of_Brackets(n,d) - Number_of_Linear_Equations(n,d) - 1);
        !            51:   end Number_of_Zero_Brackets;
        !            52:
        !            53:   function Create ( b1,b2 : Bracket ) return Bracket_Term is
        !            54:
        !            55:     bm : Bracket_Monomial := Create(b1);
        !            56:     bt : Bracket_Term;
        !            57:
        !            58:   begin
        !            59:     Multiply(bm,b2);
        !            60:     bt.coeff := Create(1.0);
        !            61:     bt.monom := bm;
        !            62:     return bt;
        !            63:   end Create;
        !            64:
        !            65:   procedure Enumerate ( n,d,k,start : in natural; accu : in out Bracket;
        !            66:                         cnt : in out natural ) is
        !            67:
        !            68:   -- DESCRIPTION :
        !            69:   --   Lexicographic enumeration of all brackets.
        !            70:
        !            71:   begin
        !            72:     if k > d
        !            73:      then -- put(accu); new_line;
        !            74:           cnt := cnt + 1;
        !            75:      else for l in start..n loop
        !            76:             accu(k) := l;
        !            77:             Enumerate(n,d,k+1,l+1,accu,cnt);
        !            78:           end loop;
        !            79:     end if;
        !            80:   end Enumerate;
        !            81:
        !            82:   procedure Enumerate2 ( n,d,k,start : in natural; b : in Bracket;
        !            83:                          accu : in out Bracket;
        !            84:                          cntstd,cntnonstd : in out natural;
        !            85:                          nonstd : in out Bracket_Polynomial ) is
        !            86:
        !            87:   -- DESCRIPTION :
        !            88:   --   Lexicographic enumeration of all brackets, with b < accu and with
        !            89:   --   a test whether the pair b*accu forms a Standard monomial.
        !            90:
        !            91:   -- ON ENTRY :
        !            92:   --   n             total number of indices to choose from;
        !            93:   --   d             number of indices in the brackets;
        !            94:   --   k             current entry in accu that is to be filled;
        !            95:   --   start         indices will be taken in between start..n;
        !            96:   --   b             previously enumerated bracket;
        !            97:   --   accu          accumulating parameter, filled in upto (k-1)th entry;
        !            98:   --   cntstd        current number of quadratic standard monomials;
        !            99:   --   cntnonstd     current number of quadratic nonstandard monomials;
        !           100:   --   nonstd        current polynomial of quadratic nonstandard monomials.
        !           101:
        !           102:   -- ON RETURN :
        !           103:   --   accu          updated accumulating parameter, accu(k) is filled in;
        !           104:   --   cnstd         updated number of quadratic standard monomials;
        !           105:   --   cntnonstd     updated number of quadratic nonstandard monomials;
        !           106:   --   nonstd        updated polynomial of quadratic nonstandard monomials.
        !           107:
        !           108:     s : natural;
        !           109:
        !           110:   begin
        !           111:     if k > d
        !           112:      then if b < accu
        !           113:            then -- put(b); put(accu);
        !           114:                 s := Brackets.Is_Standard(b,accu);
        !           115:                 if s = 0
        !           116:                  then -- put_line(" is a Standard monomial.");
        !           117:                       cntstd := cntstd + 1;
        !           118:                  else -- put(" is not a Standard monomial with s = ");
        !           119:                       -- put(s,1); new_line;
        !           120:                       cntnonstd := cntnonstd + 1;
        !           121:                       Add(nonstd,Create(b,accu));
        !           122:                 end if;
        !           123:           end if;
        !           124:      else for l in start..n loop
        !           125:             accu(k) := l;
        !           126:             Enumerate2(n,d,k+1,l+1,b,accu,cntstd,cntnonstd,nonstd);
        !           127:           end loop;
        !           128:     end if;
        !           129:   end Enumerate2;
        !           130:
        !           131:   procedure Enumerate1 ( n,d,k,start : natural; acc1,acc2 : in out Bracket;
        !           132:                          cntstd,cntnonstd : in out natural;
        !           133:                          nonstd : in out Bracket_Polynomial ) is
        !           134:
        !           135:   -- DESCRIPTION :
        !           136:   --   Lexicographic enumeration of all brackets, with acc1 < acc2 and with
        !           137:   --   a test whether the pair acc1*acc2 forms a Standard monomial.
        !           138:   --   Counts the standard and nonstandard monomials and adds every
        !           139:   --   nonStandard monomial to the polynomial nonstd.
        !           140:
        !           141:   begin
        !           142:     if k > d
        !           143:      then Enumerate2(n,d,1,acc1(1),acc1,acc2,cntstd,cntnonstd,nonstd);
        !           144:      else for l in start..n loop
        !           145:             acc1(k) := l;
        !           146:             Enumerate1(n,d,k+1,l+1,acc1,acc2,cntstd,cntnonstd,nonstd);
        !           147:           end loop;
        !           148:     end if;
        !           149:   end Enumerate1;
        !           150:
        !           151:   procedure Enumerate_Pairs ( n,d : in natural;
        !           152:                               nonstd : in out Bracket_Polynomial ) is
        !           153:
        !           154:   -- DESCRIPTION :
        !           155:   --   Enumerates all pairs (b1,b2), with b1 < b2 and checks whether
        !           156:   --   the corresponding bracket monomial b1*b2 is standard or not.
        !           157:
        !           158:     b1,b2 : Bracket(1..d);
        !           159:     cntstd,cntnonstd : natural := 0;
        !           160:
        !           161:   begin
        !           162:     Enumerate1(n,d,1,1,b1,b2,cntstd,cntnonstd,nonstd);
        !           163:     put("#Standard bi-monomials    : "); put(cntstd,1);    new_line;
        !           164:     put("#nonStandard bi-monomials : "); put(cntnonstd,1); new_line;
        !           165:   end Enumerate_Pairs;
        !           166:
        !           167:   procedure Enumerate_Brackets ( n,d : in natural ) is
        !           168:
        !           169:     acc : Bracket(1..d);
        !           170:     cnt : natural := 0;
        !           171:     nonstd : Bracket_Polynomial;
        !           172:
        !           173:   begin
        !           174:     Enumerate(n,d,1,1,acc,cnt);
        !           175:     put("#brackets : "); put(cnt,1); new_line;
        !           176:     Enumerate_Pairs(n,d,nonstd);
        !           177:     put_line("The polynomial of all quadratic nonstandard monomials :");
        !           178:     put(nonstd);
        !           179:   end Enumerate_Brackets;
        !           180:
        !           181:   procedure Read_Bracket ( b : in out Bracket ) is
        !           182:
        !           183:     d : constant natural := b'last;
        !           184:     sig : integer;
        !           185:
        !           186:   begin
        !           187:     put("Give "); put(d,1); put(" row indices : "); get(b,sig);
        !           188:     put("The sorted bracket : ");
        !           189:     if sig > 0
        !           190:      then put("+");
        !           191:      else put("-");
        !           192:     end if;
        !           193:     put(b);
        !           194:     if Is_Zero(b)
        !           195:      then put_line(" is zero.");
        !           196:      else put_line(" is different from zero.");
        !           197:     end if;
        !           198:   end Read_Bracket;
        !           199:
        !           200:   procedure Test_Straighten ( n,d : in natural ) is
        !           201:
        !           202:     b,bb : Bracket(1..d);
        !           203:     s : natural;
        !           204:     ans : character;
        !           205:     bp : Bracket_Polynomial;
        !           206:
        !           207:   begin
        !           208:     loop
        !           209:       Read_Bracket(b);
        !           210:       Read_Bracket(bb);
        !           211:       if Is_Equal(b,bb)
        !           212:        then put_line("Both brackets are equal.");
        !           213:       end if;
        !           214:       if b < bb
        !           215:        then put(b); put(" < "); put(bb);
        !           216:             s := Brackets.Is_Standard(b,bb);
        !           217:             put(" and "); put(b); put(bb);
        !           218:             bp := Straightening_Syzygy(b,bb);
        !           219:        else put(b); put(" >= "); put(bb);
        !           220:             s := Brackets.Is_Standard(bb,b);
        !           221:             put(" and "); put(bb); put(b);
        !           222:             bp := Straightening_Syzygy(bb,b);
        !           223:       end if;
        !           224:       if s = 0
        !           225:        then put_line(" is a Standard monomial.");
        !           226:        else put(" is not a Standard monomial, with s = "); put(s,1); new_line;
        !           227:       end if;
        !           228:       put_line("The straightening syzygy : "); put(bp);
        !           229:       put("Do you want to test more pairs of brackets ? (y/n) "); get(ans);
        !           230:       exit when ans /= 'y';
        !           231:     end loop;
        !           232:   end Test_Straighten;
        !           233:
        !           234:   procedure Test_Laplace ( n,d : natural ) is
        !           235:
        !           236:     p : Bracket_Polynomial := Laplace_Expansion(n,d);
        !           237:
        !           238:   begin
        !           239:     put("The Laplace expansion of "); put(n,1); put("*"); put(n,1);
        !           240:     put("-determinant as product of "); put(d,1); put("- and ");
        !           241:     put(n-d,1); put_line("-blocks : ");
        !           242:     put(p);
        !           243:   end Test_Laplace;
        !           244:
        !           245:   function Decompose ( n,d : natural; p : Bracket_Polynomial )
        !           246:                      return Bracket_Polynomial is
        !           247:
        !           248:   -- DESCRIPTION :
        !           249:   --   Decomposes the ideal of square-free nonStandard monomials.
        !           250:
        !           251:     res : Bracket_Polynomial;
        !           252:     np : constant natural := Number_of_Monomials(p);
        !           253:     brackmat : array(1..np,1..2) of Link_to_Bracket;
        !           254:
        !           255:     procedure Initialize is
        !           256:
        !           257:     -- DESCRIPTION :
        !           258:     --   Initializes the matrix of the brackets that occur in p.
        !           259:
        !           260:       cnt_term : natural := 0;
        !           261:
        !           262:       procedure List_Term ( t : in Bracket_Term; continue : out boolean ) is
        !           263:
        !           264:         cnt_mon : natural := 0;
        !           265:
        !           266:         procedure Store_Bracket ( b : in Bracket; continue : out boolean ) is
        !           267:         begin
        !           268:           cnt_mon := cnt_mon + 1;
        !           269:           brackmat(cnt_term,cnt_mon) := new Bracket'(b);
        !           270:           continue := true;
        !           271:         end Store_Bracket;
        !           272:         procedure Store_Brackets is
        !           273:           new Bracket_Monomials.Enumerate_Brackets(Store_Bracket);
        !           274:
        !           275:       begin
        !           276:         cnt_term := cnt_term+1;
        !           277:         Store_Brackets(t.monom);
        !           278:         continue := true;
        !           279:       end List_Term;
        !           280:       procedure List_Terms is new Enumerate_Terms(List_Term);
        !           281:
        !           282:     begin
        !           283:       List_Terms(p);
        !           284:     end Initialize;
        !           285:
        !           286:     procedure Write is
        !           287:
        !           288:     -- DESCRIPTION :
        !           289:     --   Writes the matrix of brackets.
        !           290:
        !           291:     begin
        !           292:       for i in 1..np loop
        !           293:         for j in 1..2 loop
        !           294:           put("  ");
        !           295:           put("b("); put(i,1); put(","); put(j,1); put(") :");
        !           296:           put(brackmat(i,j).all);
        !           297:         end loop;
        !           298:         new_line;
        !           299:       end loop;
        !           300:     end Write;
        !           301:
        !           302:     function Occured_Yet ( k : natural; b : Bracket ) return boolean is
        !           303:
        !           304:     -- DESCRIPTION :
        !           305:     --   Returns true if the bracket b occurs in one of the rows <= k.
        !           306:
        !           307:     begin
        !           308:       for i in 1..k loop
        !           309:         for j in 1..2 loop
        !           310:           if Is_Equal(brackmat(i,j).all,b)
        !           311:            then return true;
        !           312:           end if;
        !           313:         end loop;
        !           314:       end loop;
        !           315:       return false;
        !           316:     end Occured_Yet;
        !           317:
        !           318:     procedure Solve is
        !           319:
        !           320:     -- DESCRIPTION :
        !           321:     --   Solves the initial ideal of quadratic nonStandard monomials.
        !           322:
        !           323:       accu : Bracket_Monomial;
        !           324:       nzrb : constant natural := Number_of_Zero_Brackets(n,d);
        !           325:
        !           326:       procedure Solve ( k : in natural; acc : in Bracket_Monomial ) is
        !           327:       begin
        !           328:         if k > np
        !           329:          then Add(res,Create(acc));
        !           330:          else if Divisible(acc,brackmat(k,1).all)
        !           331:                 or else Divisible(acc,brackmat(k,2).all)
        !           332:                then Solve(k+1,acc);
        !           333:                else if Number_of_Brackets(acc) < nzrb
        !           334:                      then for i in 1..2 loop
        !           335:                             declare
        !           336:                               newacc : Bracket_Monomial;
        !           337:                             begin
        !           338:                               Copy(acc,newacc);
        !           339:                               Multiply(newacc,brackmat(k,i).all);
        !           340:                               Solve(k+1,newacc);
        !           341:                             end;
        !           342:                           end loop;
        !           343:                     end if;
        !           344:               end if;
        !           345:         end if;
        !           346:       end Solve;
        !           347:
        !           348:     begin
        !           349:       Solve(1,accu);
        !           350:     end Solve;
        !           351:
        !           352:   begin
        !           353:     Initialize;
        !           354:    -- Write;
        !           355:     Solve;
        !           356:     return res;
        !           357:   end Decompose;
        !           358:
        !           359:   procedure Enumerate_Straightening_Syzygies ( n,d : in natural ) is
        !           360:
        !           361:   -- DESCRIPTION :
        !           362:   --   Sets up of the straightening syzygies for all nonStandard
        !           363:   --   quadratic monomials.
        !           364:
        !           365:     nonstd,decom : Bracket_Polynomial;
        !           366:
        !           367:     procedure Write_Syzygy ( s : in Bracket_Polynomial; cont : out boolean ) is
        !           368:     begin
        !           369:       put(s); cont := true;
        !           370:     end Write_Syzygy;
        !           371:     procedure Write_Syzygies is new Enumerate_Syzygies(Write_Syzygy);
        !           372:
        !           373:   begin
        !           374:     put("(n,d) = "); put("("); put(n,1); put(","); put(d,1); put(")");
        !           375:     put("  #Brackets : "); put(Number_of_Brackets(n,d),1);
        !           376:     put("  #Linear Equations : ");
        !           377:     put(Number_of_LInear_Equations(n,d),1); new_line;
        !           378:     put_line("Type of linear equation :");
        !           379:     put(Laplace_Expansion(n,n-d));
        !           380:     nonstd := nonStandard_Monomials(n,d);
        !           381:     put_line("The polynomial with all quadratic nonStandard monomials :");
        !           382:     put(nonstd);
        !           383:     put_line("All quadratic straightening syzygies : ");
        !           384:     Write_Syzygies(nonstd);
        !           385:     put("#Zero-Brackets : "); put(Number_of_Zero_Brackets(n,d),1); new_line;
        !           386:     decom := Decompose(n,d,nonstd);
        !           387:     put_line("Decomposition of the ideal of nonStandard monomials :");
        !           388:     put(decom);
        !           389:     put("Number of solutions : "); put(Number_of_Monomials(decom),1);
        !           390:   end Enumerate_Straightening_Syzygies;
        !           391:
        !           392:   procedure Main is
        !           393:
        !           394:     n,d : natural;
        !           395:     ans : character;
        !           396:
        !           397:   begin
        !           398:     new_line;
        !           399:     put_line("Interactive testing of straightening algorithm.");
        !           400:     loop
        !           401:       new_line;
        !           402:       put_line("Choose one of the following : ");
        !           403:       put_line("  0. Exit this program.");
        !           404:       put_line("  1. Enumerate all brackets.");
        !           405:       put_line("  2. Test the straightening algorithm.");
        !           406:       put_line("  3. Test Laplace expansion.");
        !           407:       put_line("  4. Enumerate straightening syzygies.");
        !           408:       put("Make your choice (0,1,2,3,or 4) : "); get(ans);
        !           409:       exit when ans = '0';
        !           410:       new_line;
        !           411:       put("Give the number of entries in bracket : "); get(d);
        !           412:       put("Give the number of elements to choose from : "); get(n);
        !           413:       case ans is
        !           414:         when '1' => Enumerate_Brackets(n,d);
        !           415:         when '2' => Test_Straighten(n,d);
        !           416:         when '3' => Test_Laplace(n,d);
        !           417:         when '4' => Enumerate_Straightening_Syzygies(n,d);
        !           418:         when others => put_line("Bad answer.");
        !           419:       end case;
        !           420:     end loop;
        !           421:   end Main;
        !           422: begin
        !           423:   Main;
        !           424: end ts_straighten;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>