Package $NotBuiltin -- some more CoCoA-4 function

//Export NewList; // BI  8 mar 2011
Export incr;
Export decr;
Export compts;
Export DF;
Export diff;
Export tail;
//Export Last; // BI  8 mar 2011
Export cyclotomic;
//Export CoeffListWRT; // BI 2012-12-20
Export indent;
Export PrintHCIndex;
Export PrintHGIndices;
Export IntersectList;
Export interreduce;
Export interreduced;
Export MinGensGeneral;
Export minimalize;
Export Bool01;
Export AsINT;
Export AsRAT;
Export ReducedGBasis;

----------------------------------------------------------------------

Define incr(ref N)
  If type(N) <> INT Then error("value must be integer"); EndIf;
  N := N+1;
EndDefine; -- incr
  
Define decr(ref N)
  If type(N) <> INT Then error("value must be integer"); EndIf;
  N := N-1;
EndDefine; -- decr
  
Define compts(v)
  If type(v) <> MODULEELEM Then error("Expected a MODULEELEM"); EndIf;
  Return [v[i] | i In 1..NumCompts(v)];
EndDefine; -- compts


Define Bool01(b)
  If type(b) <> BOOL Then error("Arg must be a boolean"); EndIf;
  If b Then Return 1; Else Return 0; EndIf;
EndDefine; -- Bool01


Define DF(G)
  If IsZero(G) Then Return G; EndIf;
  If GradingDim(RingOf(G))=0 Then Return G; EndIf;
  Return LF(G);
EndDefine; -- DF


Define AsINT(r)
  If type(r) = INT Then Return r; EndIf;
  If type(r) = RAT Then
    If den(r)=1 Then Return num(r);
    Else error("Cannot convert to INT");
    EndIf
  EndIf;
  If type(r) <> RINGELEM Then error("Arg must be INT, RAT or RINGELEM"); EndIf;

  If not(IsInteger(ref N, r)) Then error("Cannot convert to INT"); EndIf;
  Return N;
EndDefine; -- AsINT


Define AsRAT(r)
  TopLevel ZZ;
  If type(r) = INT Then discard := IsRational(ref q, r*one(ZZ)); Return q; EndIf;
  If type(r) = RAT Then Return r; EndIf;
  If type(r) <> RINGELEM Then error("Arg must be INT, RAT or RINGELEM"); EndIf;

  If not(IsRational(ref q, r)) Then error("Cannot convert to RAT"); EndIf;
  Return q;
EndDefine; -- AsRAT


Define indent(L, Opt level)
  If type(L)=STRING Then PrintLn L; Return; EndIf;
  If not(IsDefined(level)) Then level := 1; EndIf;
  IndentRec(L, level, "", "", "\n");
EndDefine; -- indent


Define IndentRec(L, level, spaces, before, after)
  If level>0 Then 
    If type(L) = IDEAL Then
      L := gens(L);
      If L=[] Then Print spaces, before, "ideal()", after;
      Else
	PrintLn spaces, before, "ideal(";
	For I := 1 To len(L)-1 Do
	  IndentRec(L[I], level-1, spaces+"  ", "", ",\n");
	EndFor;
	IndentRec(last(L), level-1, spaces+"  ", "", "\n");
	Print spaces, ")", after;
      EndIf;
      Return;
    Elif type(L) = MODULE Then
      L := gens(L);
      If L=[] Then Print spaces, before, "module()", after;
      Else
	PrintLn spaces, before, "SubmoduleRows(F, matrix([";
	For I := 1 To len(L)-1 Do
	  IndentRec(L[I], level-1, spaces+"  ", "", ",\n");
	EndFor;
	IndentRec(last(L), level-1, spaces+"  ", "", "\n");
	Print spaces, "]))", after;
      EndIf;
      Return;
    Elif type(L) = LIST Then
      If L=[] Then Print spaces, before, "[]", after;
      Else
	PrintLn spaces, before, "[";
	For I := 1 To len(L)-1 Do
	  IndentRec(L[I], level-1, spaces+"  ", "", ",\n");
	EndFor;
	IndentRec(last(L), level-1, spaces+"  ", "", "\n");
	Print spaces, "]", after;
      EndIf;
      Return;
    Elif type(L) = RECORD Then
      F := fields(L);
      If F=[] Then Print spaces, before, "Record[]", after;
      Else
	PrintLn spaces, before, "Record[";
	For I := 1 To len(F)-1 Do
	  IndentRec(L[F[I]], level-1, spaces+"  ", F[I]+" := ", ",\n");
	EndFor;
	IndentRec(L[last(F)], level-1, spaces+"  ", F[len(F)]+" := ", "\n");
	Print spaces, "]", after;
      EndIf;
      Return;
    EndIf;
  EndIf;
  Print spaces, before;
  If type(L)=STRING Then Print "\"",L,"\""; Else Print L; EndIf;
  Print after;
EndDefine; -- indent



Define diff(L1, L2)
  If type(L1) <> LIST Or type(L2) <> LIST Then
    error("diff: arguments must be lists");
  EndIf;
  Return [ X In L1 | not(X IsIn L2) ];
EndDefine; -- diff
  
Define tail(L)  Return last(L, len(L)-1); EndDefine;

Define IntersectList(...)
  If len(ARGV) = 0 Then error("IntersectList: need at least 1 arg"); EndIf;
  If len(ARGV) > 1 Then L := ARGV; Else L := ARGV[1]; EndIf;
  ans := L[1];
  For i := 2 To len(L) Do
    ans := intersect(ans, L[i]);
  EndFor;
  Return ans;
EndDefine; -- IntersectList

Define interreduced(L)
  TmpL := diff(L, [0]);
  If TmpL=[] Then Return TmpL; EndIf;
  While True Do
    SortBy(ref TmpL, Func(f,g) Return LT(f)<LT(g); EndFunc);
    NewL := [];
    For I := 1 To len(TmpL) Do
      R := NR(TmpL[I], NewL);
      If R <> 0 Then append(ref NewL, R); EndIf;
    EndFor;
    If EqSet(TmpL, NewL) Then Return NewL; EndIf;
    TmpL := NewL;
  EndWhile;
EndDefine; -- interreduced

Define interreduce(ref L)  L := interreduced(L); EndDefine;

Define MinGensGeneral(X)
  L := gens(X);
  If type(X)=IDEAL Then 
    For I := 1 To len(L) Do
      J := ideal(WithoutNth(L,I));
      If L[I] IsIn J Then L[I] := 0*L[I]; EndIf; -- removes redundant generator
    EndFor;
  EndIf;
  If type(X)=MODULE Then 
    For I := 1 To len(L) Do
      J := submodule(ModuleOf(X), WithoutNth(L,I));
      If L[I] IsIn J Then L[I] := 0*L[I]; EndIf; -- removes redundant generator
    EndFor;
  EndIf;
  Return [V In L | not(IsZero(V))];
EndDefine;


Define minimalize(ref I)
  I := minimalized(I);
EndDefine; -- minimalize

  
Define cyclotomic(K, x)
  If shape([K, x]) <> [INT, RINGELEM] Then
    error("Args must be an INT and an indeterminate");
  EndIf;
  If not(IsIndet(x)) Then
    error("Second arg must be an indeterminate");
  EndIf;
  If K = 1 Then Return x-1; EndIf;                        -- K is 1
  L := divisors(K);
  If L = [1] Then Return sum([x^I| I In 0..(K-1)]); EndIf;-- K is prime
  M := MaxDiv(L);
  If len(M) = 1 Then Return (x^K-1)/(x^M[1]-1); EndIf;    -- K is prime power
  Return GCDList([ (x^K-1)/(x^D-1) | D In M ]);               -- general case
EndDefine; -- cyclotomic


Define divisors(K)
  Return [ I In 1..(K-1) | mod(K, I) = 0 ];
EndDefine; -- divisors

Define MaxDiv(L)
  LenL := len(L);
  Return [ L[I] | I In 1..LenL   And
	  [ X In last(L, LenL-I) | mod(X, L[I])=0]=[] ];
EndDefine; -- MaxDiv


Define PrintHGIndices(L)
  SortBy(Ref L, Func(A,B) Return A.C>B.C; EndFunc);
  indent(first([X In L | X.T <> ""], 3));
  PrintLn "..... -> total citations: ", sum([X.C | X In L]);
  PrintLn "H-index: ", max([ I In 1..len(L) | L[I].C >= I]);
  PrintLn "G-index: ",
  max([ I In 1..len(L) | sum([X.C | X In [L[A]|A In 1..I]]) >= I^2 ]);
EndDefine; -- PrintIndices


Define PrintHCIndex(L,YYYY)
  LL := L;
  For i:=1 To len(LL) Do
    LL[i].S := (4/(YYYY-LL[i].Y+1))*LL[i].C;
  EndFor;
  SortBy(Ref LL, Func(A,B) Return A.S>B.S; EndFunc);
  indent(first([X In LL | X.T <> ""], 3));
  PrintLn "..... -> total citations/age: ",
  DecimalStr(sum([X.C | X In LL]/(YYYY-1900-min([X.Y|X In LL]))));
  PrintLn "HC-index: ", max([ I In 1..len(LL) | LL[I].S >= I]);
EndDefine; -- PrintHCIndex


Define ReducedGBasis(I)
  Return [monic(F) | F In interreduced(GBasis(I))];
EndDefine; -- ReducedGBasis


Define LinKer(M)
  K := RingOf(M);
  If not(IsField(K)) Then error("NYI: Not yet implemented"); EndIf;
  IndetNames := [ Record[head := "x", indices := [N]] | N In 1..NumCols(M) ];
  Kx := NewPolyRing(K, IndetNames);
  Sys := apply(CanonicalHom(K,Kx), M) * ColMat(indets(Kx));
  I := ideal(GetCol(Sys,1));
  Params := diff(indets(Kx), LT(I));
  Basis := [];
  Foreach X In Params Do
    J := I + ideal(Kx,diff(Params,[X])) + ideal(X-1);
    append(ref Basis, [NF(indet(Kx,N), J) | N In 1..NumIndets(Kx)]);
  EndForeach;
  LC0 := Func(F) If F=0 Then Return zero(K); Else Return LC(F); EndIf; EndFunc;
  Return [[ LC0(F) | F In Vec] | Vec In Basis];
EndDefine; -- LinKer

EndPackage;
