Annotation of OpenXM_contrib/PHC/Ada/Schubert/curves_into_grassmannian.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 2: with Standard_Natural_Vectors;
! 3: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 4:
! 5: package body Curves_into_Grassmannian is
! 6:
! 7: -- CREATOR :
! 8:
! 9: function Symbolic_Create ( m,p,q : natural; top,bottom : Bracket )
! 10: return Standard_Complex_Poly_Matrices.Matrix is
! 11:
! 12: res : Standard_Complex_Poly_Matrices.Matrix(1..(m+p),1..p);
! 13: rws : constant natural := (m+p)*(q+1);
! 14: n : constant natural := Number_of_Variables(top,bottom) + 2;
! 15: row,ind,s_deg,t_deg : natural;
! 16: t : Term;
! 17:
! 18: begin
! 19: for i in res'range(1) loop -- initialization
! 20: for j in res'range(2) loop
! 21: res(i,j) := Null_Poly;
! 22: end loop;
! 23: end loop;
! 24: t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
! 25: t.cf := Create(1.0);
! 26: ind := 0; -- ind counts #variables
! 27: for j in 1..p loop -- assign columnwise
! 28: t_deg := (bottom(j)-1)/(m+p); -- degree in t to homogenize
! 29: row := 0; s_deg := 0;
! 30: for i in 1..rws loop
! 31: row := row + 1;
! 32: if i >= top(j) and i <= bottom(j)
! 33: then ind := ind+1;
! 34: t.dg(n-1) := s_deg; t.dg(n) := t_deg;
! 35: t.dg(ind) := 1;
! 36: Add(res(row,j),t);
! 37: t.dg(n-1) := 0; t.dg(n) := 0;
! 38: t.dg(ind) := 0;
! 39: end if;
! 40: if i mod (m+p) = 0
! 41: then row := 0; s_deg := s_deg+1;
! 42: if t_deg > 0
! 43: then t_deg := t_deg-1;
! 44: end if;
! 45: end if;
! 46: end loop;
! 47: end loop;
! 48: Clear(t);
! 49: return res;
! 50: end Symbolic_Create;
! 51:
! 52: -- SELECTORS :
! 53:
! 54: function Number_of_Variables ( top,bottom : Bracket ) return natural is
! 55:
! 56: cnt : natural := 0;
! 57:
! 58: begin
! 59: for j in top'range loop
! 60: cnt := cnt + (bottom(j) - top(j) + 1);
! 61: end loop;
! 62: return cnt;
! 63: end Number_of_Variables;
! 64:
! 65: function Standard_Coordinate_Frame
! 66: ( m,p,q : natural; top,bottom : Bracket;
! 67: coeff : Standard_Complex_Matrices.Matrix )
! 68: return Standard_Natural_Matrices.Matrix is
! 69:
! 70: rws : constant natural := (m+p)*(q+1);
! 71: res : Standard_Natural_Matrices.Matrix(1..rws,1..p);
! 72: tol : constant double_float := 10.0**(-10);
! 73: first : boolean;
! 74:
! 75: begin
! 76: for j in 1..p loop
! 77: first := true;
! 78: for i in 1..rws loop
! 79: if i < top(j) or i > bottom(j)
! 80: then res(i,j) := 0;
! 81: elsif (first and (AbsVal(coeff(i,j)) > tol))
! 82: then res(i,j) := 1; first := false;
! 83: else res(i,j) := 2;
! 84: end if;
! 85: end loop;
! 86: end loop;
! 87: return res;
! 88: end Standard_Coordinate_Frame;
! 89:
! 90: function Eval ( c : Term; s,t : Complex_Number ) return Term is
! 91:
! 92: res : Term;
! 93:
! 94: begin
! 95: Copy(c,res);
! 96: for i in 1..res.dg(res.dg'last-1) loop -- evaluate s
! 97: res.cf := res.cf*s;
! 98: end loop;
! 99: res.dg(res.dg'last-1) := 0;
! 100: for i in 1..res.dg(res.dg'last) loop -- evaluate t
! 101: res.cf := res.cf*t;
! 102: end loop;
! 103: res.dg(res.dg'last) := 0;
! 104: return res;
! 105: end Eval;
! 106:
! 107: function Eval ( p : Poly; s,t : Complex_Number ) return Poly is
! 108:
! 109: res : Poly := Null_Poly;
! 110:
! 111: procedure Eval_Term ( ct : in Term; continue : out boolean ) is
! 112:
! 113: et : Term := Eval(ct,s,t);
! 114:
! 115: begin
! 116: Add(res,et);
! 117: Clear(et);
! 118: continue := true;
! 119: end Eval_Term;
! 120: procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
! 121:
! 122: begin
! 123: Eval_Terms(p);
! 124: return res;
! 125: end Eval;
! 126:
! 127: function Eval ( c : Standard_Complex_Poly_Matrices.Matrix;
! 128: s,t : Complex_Number )
! 129: return Standard_Complex_Poly_Matrices.Matrix is
! 130:
! 131: res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
! 132:
! 133: begin
! 134: for i in c'range(1) loop
! 135: for j in c'range(2) loop
! 136: if c(i,j) = Null_Poly
! 137: then res(i,j) := Null_Poly;
! 138: else res(i,j) := Eval(c(i,j),s,t);
! 139: end if;
! 140: end loop;
! 141: end loop;
! 142: return res;
! 143: end Eval;
! 144:
! 145: function Elim ( c : Term; s,t : Complex_Number ) return Term is
! 146:
! 147: res : Term;
! 148:
! 149: begin
! 150: res.dg := new Standard_Natural_Vectors.Vector'
! 151: (c.dg(c.dg'first..c.dg'last-2));
! 152: res.cf := c.cf;
! 153: for i in 1..c.dg(c.dg'last-1) loop -- evaluate s
! 154: res.cf := res.cf*s;
! 155: end loop;
! 156: for i in 1..c.dg(c.dg'last) loop -- evaluate t
! 157: res.cf := res.cf*t;
! 158: end loop;
! 159: return res;
! 160: end Elim;
! 161:
! 162: function Elim ( p : Poly; s,t : Complex_Number ) return Poly is
! 163:
! 164: res : Poly := Null_Poly;
! 165:
! 166: procedure Elim_Term ( ct : in Term; continue : out boolean ) is
! 167:
! 168: et : Term := Elim(ct,s,t);
! 169:
! 170: begin
! 171: Add(res,et);
! 172: Clear(et);
! 173: continue := true;
! 174: end Elim_Term;
! 175: procedure Elim_Terms is new Visiting_Iterator(Elim_Term);
! 176:
! 177: begin
! 178: Elim_Terms(p);
! 179: return res;
! 180: end Elim;
! 181:
! 182: function Elim ( c : Standard_Complex_Poly_matrices.Matrix;
! 183: s,t : Complex_Number )
! 184: return Standard_Complex_Poly_Matrices.Matrix is
! 185:
! 186: res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
! 187:
! 188: begin
! 189: for i in c'range(1) loop
! 190: for j in c'range(2) loop
! 191: if c(i,j) = Null_Poly
! 192: then res(i,j) := Null_Poly;
! 193: else res(i,j) := Elim(c(i,j),s,t);
! 194: end if;
! 195: end loop;
! 196: end loop;
! 197: return res;
! 198: end Elim;
! 199:
! 200: function Column_Localize ( top,bottom : Bracket;
! 201: locmap : Standard_Natural_Matrices.Matrix;
! 202: t : Term ) return Term is
! 203:
! 204: -- DESCRIPTION :
! 205: -- Applies the localization map to the term, eliminating those xij's
! 206: -- xij for which the corresponding entry in locmap is either 0 or 1.
! 207:
! 208: -- NOTE :
! 209: -- This localization assumes that t.dg(k) = 0 with k for which the
! 210: -- corresponding (i,j) with locmap(i,j) = 0.
! 211: -- The localization pattern is traversed columnwise.
! 212:
! 213: res : Term;
! 214: ndg : Standard_Natural_Vectors.Vector(t.dg'range);
! 215: cnt : natural := t.dg'first-1;
! 216: ind : natural := cnt;
! 217:
! 218: begin
! 219: for j in locmap'range(2) loop -- columnwise order of the variables
! 220: for i in top(j)..bottom(j) loop -- restricted range skips the zeros
! 221: ind := ind+1;
! 222: if locmap(i,j) = 2 -- skip the ones
! 223: then cnt := cnt + 1;
! 224: ndg(cnt) := t.dg(ind);
! 225: end if;
! 226: end loop;
! 227: end loop;
! 228: for i in ind+1..t.dg'last loop -- leave the lifting !
! 229: cnt := cnt+1;
! 230: ndg(cnt) := t.dg(i);
! 231: end loop;
! 232: res.cf := t.cf;
! 233: res.dg := new Standard_Natural_Vectors.Vector'(ndg(1..cnt));
! 234: return res;
! 235: end Column_Localize;
! 236:
! 237: function Column_Localize ( top,bottom : Bracket;
! 238: locmap : Standard_Natural_Matrices.Matrix;
! 239: p : Poly ) return Poly is
! 240:
! 241: -- DESCRIPTION :
! 242: -- Applies the localization map to the polynomial, eliminating
! 243: -- those xij's for which locmap(i,j) is either 0 or 1.
! 244:
! 245: res : Poly := Null_Poly;
! 246:
! 247: procedure Column_Localize_Term ( t : in Term; continue : out boolean ) is
! 248:
! 249: lt : Term := Column_Localize(top,bottom,locmap,t);
! 250:
! 251: begin
! 252: Add(res,lt);
! 253: Clear(lt.dg);
! 254: continue := true;
! 255: end Column_Localize_Term;
! 256: procedure Column_Localize_Terms is
! 257: new Visiting_Iterator(Column_Localize_Term);
! 258:
! 259: begin
! 260: Column_Localize_Terms(p);
! 261: return res;
! 262: end Column_Localize;
! 263:
! 264: function Column_Localize ( top,bottom : Bracket;
! 265: locmap : Standard_Natural_Matrices.Matrix;
! 266: p : Poly_Sys ) return Poly_Sys is
! 267:
! 268: res : Poly_Sys(p'range);
! 269:
! 270: begin
! 271: for i in p'range loop
! 272: res(i) := Column_Localize(top,bottom,locmap,p(i));
! 273: end loop;
! 274: return res;
! 275: end Column_Localize;
! 276:
! 277: function Column_Vector_Rep
! 278: ( top,bottom : Bracket;
! 279: cffmat : Standard_Complex_Matrices.Matrix )
! 280: return Standard_Complex_Vectors.Vector is
! 281:
! 282: dim : constant natural := Number_of_Variables(top,bottom);
! 283: res : Standard_Complex_Vectors.Vector(1..dim);
! 284: cnt : natural := 0;
! 285:
! 286: begin
! 287: for j in cffmat'range(2) loop
! 288: for i in top(j)..bottom(j) loop
! 289: cnt := cnt + 1;
! 290: res(cnt) := cffmat(i,j);
! 291: end loop;
! 292: end loop;
! 293: return res;
! 294: end Column_Vector_Rep;
! 295:
! 296: function Column_Vector_Rep ( locmap : Standard_Natural_Matrices.Matrix;
! 297: cffmat : Standard_Complex_Matrices.Matrix )
! 298: return Standard_Complex_Vectors.Vector is
! 299:
! 300: dim : constant natural := cffmat'length(1)*cffmat'length(2);
! 301: res : Standard_Complex_Vectors.Vector(1..dim);
! 302: cnt : natural := 0;
! 303:
! 304: begin
! 305: for j in cffmat'range(2) loop
! 306: for i in cffmat'range(1) loop
! 307: if locmap(i,j) = 2
! 308: then cnt := cnt + 1;
! 309: res(cnt) := cffmat(i,j);
! 310: end if;
! 311: end loop;
! 312: end loop;
! 313: return res(1..cnt);
! 314: end Column_Vector_Rep;
! 315:
! 316: function Column_Matrix_Rep
! 317: ( locmap : Standard_Natural_Matrices.Matrix;
! 318: cffvec : Standard_Complex_Vectors.Vector )
! 319: return Standard_Complex_Matrices.Matrix is
! 320:
! 321: res : Standard_Complex_Matrices.Matrix(locmap'range(1),locmap'range(2));
! 322: cnt : natural := 0;
! 323:
! 324: begin
! 325: for j in locmap'range(2) loop
! 326: for i in locmap'range(1) loop
! 327: if locmap(i,j) = 0
! 328: then res(i,j) := Create(0.0);
! 329: elsif locmap(i,j) = 1
! 330: then res(i,j) := Create(1.0);
! 331: else cnt := cnt + 1;
! 332: res(i,j) := cffvec(cnt);
! 333: end if;
! 334: end loop;
! 335: end loop;
! 336: return res;
! 337: end Column_Matrix_Rep;
! 338:
! 339: procedure Swap ( p : in out Poly; k,l : in natural ) is
! 340:
! 341: procedure Swap_Term ( t : in out Term; continue : out boolean ) is
! 342:
! 343: lval : natural := t.dg(l);
! 344:
! 345: begin
! 346: t.dg(l) := t.dg(k);
! 347: t.dg(k) := lval;
! 348: continue := true;
! 349: end Swap_Term;
! 350: procedure Swap_Terms is new Changing_Iterator(Swap_Term);
! 351:
! 352: begin
! 353: Swap_Terms(p);
! 354: end Swap;
! 355:
! 356: procedure Swap ( c : in out Standard_Complex_Poly_Matrices.Matrix;
! 357: k,l : in natural ) is
! 358: begin
! 359: for i in c'range(1) loop
! 360: for j in c'range(2) loop
! 361: if c(i,j) /= Null_Poly
! 362: then Swap(c(i,j),k,l);
! 363: end if;
! 364: end loop;
! 365: end loop;
! 366: end Swap;
! 367:
! 368: function Insert ( p : Poly; k : natural ) return Poly is
! 369:
! 370: res : Poly := Null_Poly;
! 371:
! 372: procedure Insert_Term ( t : in Term; continue : out boolean ) is
! 373:
! 374: rt : Term;
! 375:
! 376: begin
! 377: rt.cf := t.cf;
! 378: rt.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
! 379: rt.dg(t.dg'first..k-1) := t.dg(t.dg'first..k-1);
! 380: rt.dg(k) := 0;
! 381: rt.dg(k+1..rt.dg'last) := t.dg(k..t.dg'last);
! 382: Add(res,rt);
! 383: Clear(rt);
! 384: continue := true;
! 385: end Insert_Term;
! 386: procedure Insert_Terms is new Visiting_Iterator(Insert_Term);
! 387:
! 388: begin
! 389: Insert_Terms(p);
! 390: return res;
! 391: end Insert;
! 392:
! 393: function Insert ( c : Standard_Complex_Poly_Matrices.Matrix; k : natural )
! 394: return Standard_Complex_Poly_Matrices.Matrix is
! 395:
! 396: res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
! 397:
! 398: begin
! 399: for i in c'range(1) loop
! 400: for j in c'range(2) loop
! 401: if c(i,j) = Null_Poly
! 402: then res(i,j) := Null_Poly;
! 403: else res(i,j) := Insert(c(i,j),k);
! 404: end if;
! 405: end loop;
! 406: end loop;
! 407: return res;
! 408: end Insert;
! 409:
! 410: procedure Duplicate ( p : in out Poly; k,l : in natural ) is
! 411:
! 412: procedure Duplicate_Term ( t : in out Term; continue : out boolean ) is
! 413: begin
! 414: t.dg(l) := t.dg(k);
! 415: continue := true;
! 416: end Duplicate_Term;
! 417: procedure Duplicate_Terms is new Changing_Iterator(Duplicate_Term);
! 418:
! 419: begin
! 420: Duplicate_Terms(p);
! 421: end Duplicate;
! 422:
! 423: end Curves_into_Grassmannian;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to curves_into_grassmannian.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Schubert