read "tools.txt";


OpeningOut := proc(Slot:: list)
          local m,s,w,z,a,b;
	options trace;
          a := evalf(Slot[1]); b := evalf(Slot[2]);
          if not(type(a,numeric)) or not(type(b,numeric)) then
              ERROR("Give real numeric input values.")
          fi;
          if a = b then ERROR("Slots must have positive length") fi;
          if a > b then a,b := b,a; fi;
          m,s := (b+a)/2,(b-a)/4;
          unapply(piecewise(z<=m-2*s,(m+z)/2 - Re(sqrt(((m-z)/2)^2-s^2)),z>=m+2*s, (m+z)/2 + Re(sqrt(((m-z)/2)^2-s^2)),0),z);
    end;

# Here Slot is a list as follows [first end-point,second-end point, opened or not (0/1), generation]
#

#OpeningIn := proc(A,B:: numeric)
##          local m,s,w,z,a,b;
#          a := evalf(A); b := evalf(B);
#          if not(type(a,numeric)) or not(type(b,numeric)) then
#              ERROR("Give real numeric input values.")
#          fi;
#          if a = b then ERROR("Input numbers must be different") fi;
#          if a > b then a,b := b,a; fi;
#          m,s := (b+a)/2,(b-a)/4;
#          unapply(piecewise(z<=m-2*s,(m+z)/2 + sqrt(((m-z)/2)^2-s^2),z>=m+2*s, (m+z)/2 - sqrt(((m-z)/2)^2-s^2)),z);
#    end;



Iteration := proc(L::list(numeric),g::posint,d::posint)
     local Points, i,k,l,m,n, NofSlots;
     global Slots;
     m := nops(L);
     if not(type(m/2,integer)) then
       ERROR("The input must be a list with an even number of numeric elements.");
     fi;
     Points := sort(L);
     Slots := array(1..g,1..d);
     for i from 1 to m/2 do
        Slots[1,i] :=[evalf(Points[2*i-1]),evalf(Points[2*i]),1];
        od;
        n:=0;
     for i from 1 to m/2 do
     for k from 1 to m/2 do
        Slots[2,n+i] := [OpeningOut(Slots[1,k])(Slots[1,i][1]),OpeningOut(Slots[1,k])(Slots[1,i][2]),1];
        for l from 1 to n+i-1 do
        Slots[2,l]:= [OpeningOut(Slots[1,k])(Slots[1,l][1]),OpeningOut(Slots[1,k])(Slots[1,l][2]),1];
        od;
        od;
        n := k*m/2
        od;

     Slots;
     end;


Iteration2X := proc(L::list(numeric))
     local Slots, TempSlots, Points, i,k,m;
     m := nops(L);
     if not(type(m/2,integer)) then
       ERROR("The input must be a list with an even number of numeric elements.");
     fi;
     Points := sort(L);
     Slots[1] := Points;
     for i from 1 to m/2 do
        for k from 1 to m do
         Points[k] := evalf(OpeningOut(evalf(Slots[1][2*i-1]),evalf(Slots[1][2*i]))(evalf(Slots[1][k])));
        od;
        od;
        TempSlots[2] := Array() ;
     for i from 1 to m/2 do
        for k from 2*i + 1 to m/2
        do TempSlots[2]:= [op(TempSlots[2]),evalf(OpeningOut(evalf(Slots[1][2*i-1]),evalf(Slots[1][2*i])(Slots[1][k])))];
        od;
        od;
     Slots[2] := sort(TempSlots[2]);
     for i from 1 to nops(Slots[1])/2
        do
          for k from 1 to m
          do Points[k] := evalf(OpeningOut(evalf(Slots[1][2*i-1])),evalf(Slots[1][2*i])(evalf(Points[k])));
          od;
         od;
     for i from 1 to nops(Slots[2])/2
        do
          for k from 1 to nops(Points)
          do Points[k] := evalf(OpeningOut(evalf(Slots[2][2*i-1])),evalf(Slots[2][2*i])(evalf(Points[k])));
          od;
         od;
     eval(Points);
     end;


IterationManual := proc(L::list(numeric))
     local   Slots4, i, j,k,m,l,oo,ee,Points;
     global Slots1,Points1, Points2, Points3, Slots2, Slots3;
options trace;
     Points := L;
     Slots1 := L;
     Slots2 := list();
     Slots3 := list();
     Points1 := list();
     m := nops(L);
     if m<>6 then ERROR("Give six different branch point of an algebraic curve of genus 2"); fi;
     for i from 1 to m do
       for k from 1 to m do
       if j<>m and L[j]=L[k] then ERROR("Branch points must be different"); fi;
       od;
       od;
       for i  from 1 to m/2 do
       print("Round");print(i);
       print(Slots1[2*i-1]); print(Slots1[2*i]);
       oo := OpeningOut([Re(evalf(Slots1[2*i-1])),Re(evalf(Slots1[2*i]))]);
print(oo);
trace(oo);
       print(Slots1);
       print(Slots1[2*i-1]); print("oo"); print(oo(Slots1[2*i-1]));
       print(Slots1[2*i]);print("oo"); print(oo(Slots1[2*i]));
       print("Onko oikein");
       print(oo);
       Slots1 := map(oo,Slots1);
       print(Slots1);
       Slots2 := map(oo,Slots2);
       Slots3 := map(oo,Slots3);
       Points1 := map(oo,Points);
       ee := CreateEllMoe(Slots1[2*i-1],Slots1[2*i],Pi);
       for j from 2*i+1 to m do
             Slots2 := sort([op(Slots2),ee(Slots1[j])]);
           od;
       #if i = 2 then
       #   Slots3 := sort(map(ee,Slots2));
       #fi;
       #if i=3 then
       #  Slots3 := [op(Slots3),ee(evalf(Slots2[5])),ee(evalf(Slots2[6]))];
       #fi;
       #print("Round ends");
       #print(i);
       #print(Slots3);
       od;
       Slots3 := sort(Slots3);
       Slots2 := sort(Slots2);
       Points;
    end;