[BACK]Return to transformations.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift

Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/transformations.adb, Revision 1.1.1.1

1.1       maekawa     1: with unchecked_deallocation;
                      2: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
                      3: with Standard_Common_Divisors;           use Standard_Common_Divisors;
                      4: with Standard_Integer_Linear_Solvers;    use Standard_Integer_Linear_Solvers;
                      5:
                      6: package body Transformations is
                      7:
                      8:   type Transfo_tp is new matrix;
                      9:   procedure free is new unchecked_deallocation (Transfo_tp,Transfo);
                     10:
                     11: -- CONSTRUCTORS :
                     12:
                     13:   function Id ( n : natural ) return Transfo is
                     14:
                     15:     t : Transfo;
                     16:
                     17:   begin
                     18:     t := new Transfo_tp(1..n,1..n);
                     19:     for i in 1..n loop
                     20:       for j in 1..n loop
                     21:         t(i,j) := 0;
                     22:       end loop;
                     23:       t(i,i) := 1;
                     24:     end loop;
                     25:     return t;
                     26:   end Id;
                     27:
                     28:   function Create ( v : Vector; i : integer ) return Transfo is
                     29:
                     30:     t : Transfo;
                     31:     j : integer;
                     32:
                     33:   begin
                     34:     t := Id(v'last-v'first+1);
                     35:     if v(i) = 0
                     36:      then return t;
                     37:      else j := 0;
                     38:           for k in v'range loop
                     39:             if (v(k) /= 0) and (k /= i)
                     40:              then j := k;
                     41:             end if;
                     42:             exit when j /= 0;
                     43:           end loop;
                     44:          if j = 0
                     45:           then return t;
                     46:           else declare
                     47:                  a,b,k,l,d : integer;
                     48:                 begin
                     49:                  a := v(i); b := v(j);
                     50:                  gcd(a,b,k,l,d);
                     51:                  a := a/d;  b := b/d;
                     52:                  t(i,i) := k;  t(i,j) := l;
                     53:                  t(j,i) := -b; t(j,j) := a;
                     54:                 end;
                     55:                return t;
                     56:           end if;
                     57:     end if;
                     58:   end Create;
                     59:
                     60:   function Create ( v : VecVec ) return Transfo is
                     61:
                     62:     t : Transfo;
                     63:
                     64:   begin
                     65:     t := new Transfo_tp(v'range,v'range);
                     66:     for i in v'range loop
                     67:       for j in v(i)'range loop
                     68:        t(j,i) := v(i)(j);
                     69:       end loop;
                     70:     end loop;
                     71:     return t;
                     72:   end Create;
                     73:
                     74:   function Create ( m : Matrix ) return Transfo is
                     75:
                     76:     t : Transfo;
                     77:
                     78:   begin
                     79:     t := new Transfo_tp(m'range(1),m'range(2));
                     80:     t.all := transfo_tp(m);
                     81:     return t;
                     82:   end Create;
                     83:
                     84:   function Rotate ( v : Vector; i : integer ) return Transfo is
                     85:
                     86:     t : Transfo;
                     87:     j : integer;
                     88:
                     89:   begin
                     90:     t := Id(v'last-v'first+1);
                     91:     if v(i) = 0
                     92:      then return t;
                     93:      else declare
                     94:            w : Vector(v'range) := v;
                     95:            acc : Transfo := new Transfo_tp'(t.all);
                     96:          begin
                     97:            loop
                     98:              j := 0;
                     99:              for k in w'range loop
                    100:                if (w(k) /= 0) and (k /= i)
                    101:                  then j := k;
                    102:                 end if;
                    103:                exit when j /= 0;
                    104:               end loop;
                    105:              if j /= 0
                    106:               then declare
                    107:                      a,b,k,l,d : integer;
                    108:                     begin
                    109:                      a := w(i); b := w(j);
                    110:                      gcd(a,b,k,l,d);
                    111:                      a := a/d;  b := b/d;
                    112:                      t(i,i) := k;  t(i,j) := l;
                    113:                      t(j,i) := -b; t(j,j) := a;
                    114:                     end;
                    115:                    Apply(t,w);
                    116:                    Mult2(t,acc);
                    117:                    t(i,i) := 1; t(j,j) := 1;
                    118:                    t(i,j) := 0; t(j,i) := 0;
                    119:               end if;
                    120:              exit when j = 0;
                    121:             end loop;
                    122:            Clear(t);
                    123:            return acc;
                    124:           end;
                    125:     end if;
                    126:   end Rotate;
                    127:
                    128:   function Rotate ( v : Link_to_Vector; i : integer ) return Transfo is
                    129:   begin
                    130:     if v = null
                    131:      then return null;
                    132:      else return Rotate(v.all,i);
                    133:     end if;
                    134:   end Rotate;
                    135:
                    136:   procedure Transpose ( t : in out Transfo ) is
                    137:
                    138:   -- DESCRIPTION :
                    139:   --   Transposes the matrix of t.
                    140:
                    141:     temp : integer;
                    142:
                    143:   begin
                    144:     for i in t'range(1) loop
                    145:       for j in t'first(2)..(i-1) loop
                    146:        temp := t(i,j);
                    147:        t(i,j) := t(j,i);
                    148:        t(j,i) := temp;
                    149:       end loop;
                    150:     end loop;
                    151:   end Transpose;
                    152:
                    153:   function Build_Transfo ( v : Vector; i : integer ) return Transfo is
                    154:
                    155:     t : Transfo;
                    156:     ind : integer;
                    157:     zeroes : boolean;
                    158:
                    159:   begin
                    160:     t := Id(v'last-v'first+1);
                    161:     if v(i) /= 0
                    162:      then ind := i;
                    163:      else ind := v'last + 1;
                    164:           for k in v'range loop
                    165:             if v(k) /= 0
                    166:              then ind := k;
                    167:                  exit;
                    168:             end if;
                    169:           end loop;
                    170:     end if;
                    171:     if ind = v'last + 1
                    172:      then return t;
                    173:      else declare
                    174:            w : Vector(v'range) := v;
                    175:            acc : Transfo := new Transfo_tp'(t.all);
                    176:            j1,j2 : integer;
                    177:          begin
                    178:            if v(i) = 0
                    179:             then acc(i,i) := 0; acc(i,ind) := 1;
                    180:                  acc(ind,i) := 1; acc(ind,ind) := 0;
                    181:                  w(i) := w(ind); w(ind) := 0;
                    182:             end if;
                    183:             for k in v'first..i loop
                    184:              if w(k) /= 0
                    185:               then j1 := k;
                    186:                    exit;
                    187:               end if;
                    188:             end loop;
                    189:            if j1 /= i
                    190:              then zeroes := false;
                    191:                  loop
                    192:                    for k in (j1+1)..i loop
                    193:                      if w(k) /= 0
                    194:                       then j2 := k;
                    195:                            exit;
                    196:                       end if;
                    197:                     end loop;
                    198:                    declare
                    199:                      a,b,k,l,d : integer;
                    200:                     begin
                    201:                      a := w(j1); b := w(j2);
                    202:                      gcd(a,b,k,l,d);
                    203:                      a := a/d;  b := b/d;
                    204:                      t(j1,j1) := l;  t(j1,j2) := -k;
                    205:                      t(j2,j1) := a; t(j2,j2) := b;
                    206:                      w(j2) := d;
                    207:                     end;
                    208:                    Mult2(t,acc);
                    209:                    t(j1,j1) := 1; t(j2,j2) := 1;
                    210:                    t(j1,j2) := 0; t(j2,j1) := 0;
                    211:                    if j2 < i
                    212:                     then j1 := j2;
                    213:                     else exit;
                    214:                     end if;
                    215:                   end loop;
                    216:              else zeroes := true;
                    217:             end if;
                    218:             for k in reverse i..v'last loop
                    219:              if w(k) /= 0
                    220:               then j2 := k;
                    221:                    exit;
                    222:               end if;
                    223:             end loop;
                    224:            if j2 /= i
                    225:              then loop
                    226:                    for k in reverse i..(j2-1) loop
                    227:                      if w(k) /= 0
                    228:                       then j1 := k;
                    229:                            exit;
                    230:                       end if;
                    231:                     end loop;
                    232:                    declare
                    233:                      a,b,k,l,d : integer;
                    234:                     begin
                    235:                      a := w(j1); b := w(j2);
                    236:                      gcd(a,b,k,l,d);
                    237:                      a := a/d;  b := b/d;
                    238:                      t(j1,j1) := a;  t(j1,j2) := b;
                    239:                      t(j2,j1) := -l; t(j2,j2) := k;
                    240:                      w(j1) := d;
                    241:                     end;
                    242:                    Mult2(t,acc);
                    243:                    t(j1,j1) := 1; t(j2,j2) := 1;
                    244:                    t(j1,j2) := 0; t(j2,j1) := 0;
                    245:                    if j1 > i
                    246:                     then j2 := j1;
                    247:                     else exit;
                    248:                     end if;
                    249:                   end loop;
                    250:              elsif zeroes
                    251:                 then if w(i) < 0
                    252:                       then t(i,i) := -1;
                    253:                            Mult2(t,acc);
                    254:                       end if;
                    255:            end if;
                    256:            Clear(t); --otherwise segmentation fault...
                    257:            return acc;
                    258:           end;
                    259:     end if;
                    260:   end Build_Transfo;
                    261:
                    262:   function Build_Transfo ( v : Link_to_Vector; i : integer ) return Transfo is
                    263:   begin
                    264:     if v = null
                    265:      then return null;
                    266:      else return Build_Transfo(v.all,i);
                    267:     end if;
                    268:   end Build_Transfo;
                    269:
                    270: -- SELECTOR :
                    271:
                    272:   function Dimension ( t : Transfo ) return natural is
                    273:   begin
                    274:     return (t'last(1) - t'first(1) + 1);
                    275:   end Dimension;
                    276:
                    277:   function Sign ( t : Transfo ) return integer is
                    278:   begin
                    279:     if t /= null
                    280:      then declare
                    281:             a  : matrix(t'range(1),t'range(2)) := matrix(t.all);
                    282:           begin
                    283:             return Det(a);
                    284:           end;
                    285:      else return 0;
                    286:     end if;
                    287:   end Sign;
                    288:
                    289: -- OPERATIONS :
                    290:
                    291:   function Transpose ( t : Transfo ) return Transfo is
                    292:
                    293:     res : Transfo := new Transfo_tp(t'range(2),t'range(1));
                    294:
                    295:   begin
                    296:     for i in res'range(1) loop
                    297:       for j in res'range(2) loop
                    298:         res(i,j) := t(j,i);
                    299:       end loop;
                    300:     end loop;
                    301:     return res;
                    302:   end Transpose;
                    303:
                    304:   function Invert ( t : Transfo ) return Transfo is
                    305:
                    306:     n : constant natural := Dimension(t);
                    307:     triang : Matrix(t'range(1),t'range(2)) := matrix(t.all);
                    308:     l : Matrix(t'range(1),t'range(1));
                    309:     res : Transfo := new Transfo_tp(t'range(1),t'range(2));
                    310:     x,b : Vector(1..n) := (1..n => 0);
                    311:     fail : boolean;
                    312:
                    313:   begin
                    314:     Upper_Triangulate(l,triang);
                    315:     for j in t'range(2) loop
                    316:       for i in l'range(1) loop  -- image of jth basis vector
                    317:        b(i) := l(i,j);
                    318:       end loop;
                    319:       Solve1(triang,x,b,fail);
                    320:       for i in t'range(1) loop
                    321:        res(i,j) := x(i);
                    322:       end loop;
                    323:     end loop;
                    324:     return res;
                    325:   end Invert;
                    326:
                    327:   function "*" ( t1,t2 : Transfo ) return Transfo is
                    328:
                    329:     r : Transfo;
                    330:
                    331:   begin
                    332:     r := new Transfo_tp(t1'range(1),t2'range(2));
                    333:     r.all := t1.all*t2.all;
                    334:     return r;
                    335:   end "*";
                    336:
                    337:   procedure Mult1 ( t1 : in out Transfo; t2 : in Transfo ) is
                    338:   begin
                    339:     Mul1(t1.all,t2.all);
                    340:   end Mult1;
                    341:
                    342:   procedure Mult2 ( t1 : in Transfo; t2 : in out Transfo ) is
                    343:   begin
                    344:     Mul2(t1.all,t2.all);
                    345:   end Mult2;
                    346:
                    347:   function "*" ( t : Transfo; v : Vector ) return Vector is
                    348:
                    349:     r : Vector(t'range(1));
                    350:
                    351:   begin
                    352:     r := t.all*v;
                    353:     return r;
                    354:   end "*";
                    355:
                    356:   function "*" ( t : Transfo; v : Link_to_Vector ) return Link_to_Vector is
                    357:   begin
                    358:     if v = null
                    359:      then return v;
                    360:      else declare
                    361:             res : Link_to_Vector := new Vector(t'range(1));
                    362:           begin
                    363:             res.all := t.all*v.all;
                    364:             return res;
                    365:           end;
                    366:     end if;
                    367:   end "*";
                    368:
                    369:   procedure Apply ( t : in Transfo; v : in out Vector ) is
                    370:   begin
                    371:     v := t.all*v;
                    372:   end Apply;
                    373:
                    374:   procedure Apply ( t : in Transfo; v : in out Link_to_Vector ) is
                    375:   begin
                    376:     if v /= null
                    377:      then Apply(t,v.all);
                    378:     end if;
                    379:   end Apply;
                    380:
                    381:   function "*" ( t : Transfo; v : Standard_Complex_Vectors.Vector )
                    382:                return Standard_Complex_Vectors.Vector is
                    383:
                    384:     res : Standard_Complex_Vectors.Vector(v'range);
                    385:
                    386:   begin
                    387:     for j in t'range(2) loop
                    388:       res(j) := Create(1.0);
                    389:       for i in t'range(1) loop
                    390:         res(j) := res(j)*v(i)**t(i,j);
                    391:       end loop;
                    392:     end loop;
                    393:     return res;
                    394:   end "*";
                    395:
                    396:   procedure Apply ( t : in Transfo;
                    397:                     v : in out Standard_Complex_Vectors.Vector ) is
                    398:   begin
                    399:     v := t*v;
                    400:   end Apply;
                    401:
                    402: -- DESTRUCTOR :
                    403:
                    404:   procedure Clear ( t : in out Transfo ) is
                    405:   begin
                    406:     free(t);
                    407:   end Clear;
                    408:
                    409: end Transformations;

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