Package $points -- computations with ideals of points

Export GenericPoints;
Export IdealOfPoints_OLD;

Define About()
  PrintLn "    Author: J Abbott";
  PrintLn "    Date: 1 July 1999  (corrected use of \"End\"  4 July 2002)  ";
EndDefine;

-----------------------------------------------------------------------------
-- This is only tedious interface code; the real stuff is written in C and
-- is directly accessible via the "built-in" function "BuchbergerMoeller".
-- Be careful calling that function as it does not check its argument;
-- it returns either an empty list (meaning that duplicate points were
-- supplied in the input) or a two element list whose first element is the
-- GBasis and the second is the list of separators.

Define BuchbergerMoeller(PolyR, Pts)
  error("BuchbergerMoeller: Please be patient, we're working on it");
EndDefine; -- BuchbergerMoeller


Define BuchbergerMoellerArgCheck(Ref L, PolyR, FnName)
  If L = [] Then Return; EndIf;
  ErrMsg := FnName+": argument must be a list of points with rational coords";
  If Type(L) <> LIST Then Error(ErrMsg); EndIf;
  Tmp := MakeSet([Type(X) | X In L]);
  If Tmp <> [LIST] Then Error(ErrMsg); EndIf;
  LenPt := MakeSet([Len(X) | X In L]);
  If Len(LenPt) <> 1 Then Error(FnName+": points must lie in the same space"); EndIf;
  LenPt := LenPt[1];
  If NumIndets(PolyR) < LenPt Then Error(FnName+": current ring has too few indeterminates"); EndIf;
  For I:=1 To Len(L) Do
    For J:=1 To LenPt Do
      If Type(L[I][J]) = RINGELEM And L[I][J] = LC(L[I][J]) Then
        NoPrint := IsRational(L[I][J], LC(L[I][J]));
      EndIf;
      If not(Type(L[I][J]) IsIn [INT,RAT]) Then Error(ErrMsg); EndIf;
    EndFor;
  EndFor;
EndDefine;


Define CommonDenoms(PolyR, ref L)
  If L = [] Then Return 1; EndIf;
  P := Characteristic(PolyR);
  Nvars := Len(L[1]);
  Ans := NewList(Nvars);
  Npoints := Len(L);
  For I:=1 To Nvars Do
    LCM := 1;
    For J:=1 To Npoints Do
      If Type(L[J][I]) = RAT Then
        LCM := LCM(LCM, L[J][I].Den);
      EndIf;
    EndFor;
    Ans[I] := LCM;
    If P <> 0 And Mod(LCM, P) = 0 Then Error("Point has zero denominator"); EndIf;
    For J:=1 To Npoints Do
      L[J][I] := LCM*L[J][I]; -- assumes automatic retraction to INT
      If P <> 0 Then L[J][I] := Mod(L[J][I], P); EndIf;
    EndFor;
  EndFor;
  Return Ans;
EndDefine;


Define IdealOfPoints_OLD(PolyR, P)
  If P = [] Then Return Ideal(1); EndIf;
  $points.BuchbergerMoellerArgCheck(ref P, PolyR, "IdealOfPoints");
  Denoms := $points.CommonDenoms(PolyR, ref P);
  Tmp := BuchbergerMoeller(1, P);
  If Tmp = [] Then Error("IdealOfPoints: points must be distinct"); EndIf;
  Tmp[2] := []; -- to avoid wasting memory
  Rescale := [[Indet(PolyR,J), Denoms[J]*Indet(PolyR,J)] | J In 1..Len(Denoms)];
  GB := Subst(Tmp[1], Rescale);
  GB := [Monic(F) | F In GB];
  Tmp := []; -- to avoid wasting memory
  I := Ideal(GB);
  I.GBasis := GB;
  Return I;
EndDefine;


Define IdealOfProjectivePoints(PolyR, P)
  If P = [] Then Return Ideal(1); EndIf;
  $points.BuchbergerMoellerArgCheck(P, PolyR, "IdealOfProjectivePoints");
  Denoms := $points.CommonDenoms(PolyR, ref P);
  Tmp := BuchbergerMoeller(0, P);
  If Tmp = [] Then Error("IdealOfProjectivePoints: points must be distinct (and with at least one non-zero coord)"); EndIf;
  Tmp[2] := []; -- to avoid wasting memory
  Rescale := [[Indet(PolyR,J), Denoms[J]*Indet(PolyR,J)] | J In 1..Len(Denoms)];
  GB := Subst(Tmp[1], Rescale);
  GB := [Monic(F) | F In GB];
  Tmp := []; -- to avoid wasting memory
  I := Ideal(GB);
  I.GBasis := GB;
  Return I;
EndDefine;


Define SeparatorsOfPoints(PolyR, P)
  If P = [] Then Return []; EndIf;
  $points.BuchbergerMoellerArgCheck(P, PolyR, "SeparatorsOfPoints");
  Pcopy := P;
  Denoms := $points.CommonDenoms(PolyR, ref Pcopy);
  Tmp := BuchbergerMoeller(1, Pcopy);
  If Tmp = [] Then Error("SeparatorsOfPoints: points must be distinct"); EndIf;
  Tmp[1] := []; -- to avoid wasting memory
  Rescale := [[Indet(I), Denoms[I]*Indet(I)] | I In 1..Len(Denoms)];
  S := Subst(Tmp[2], Rescale);
  Tmp := []; -- to avoid wasting memory
  For I:=1 To Len(S) Do S[I] := S[I]/Eval(S[I], P[I]); EndFor;
  Return S;
EndDefine;


Define SeparatorsOfProjectivePoints(PolyR, P)
  If P = [] Then Return []; EndIf;
  $points.BuchbergerMoellerArgCheck(P, PolyR, "SeparatorsOfProjectivePoints");
  Pcopy := P;
  Denoms := $points.CommonDenoms(PolyR, ref Pcopy);
  Tmp := BuchbergerMoeller(0, Pcopy);
  If Tmp = [] Then Error("SeparatorsOfProjectivePoints: points must be distinct (and with at least one non-zero coord)"); EndIf;
  Tmp[1] := []; -- to avoid wasting memory
  Rescale := [[Indet(I), Denoms[I]*Indet(I)] | I In 1..Len(Denoms)];
  S := Subst(Tmp[2], Rescale);
  Return S;
EndDefine;


Define IdealAndSeparatorsOfPoints(PolyR, P)
  If P = [] Then Return Record[points := P, ideal := Ideal(1), separators := []]; EndIf;
  $points.BuchbergerMoellerArgCheck(P, PolyR, "IdealAndSeparatorsOfPoints");
  Pcopy := P;
  Denoms := $points.CommonDenoms(PolyR, ref Pcopy);
  Tmp := BuchbergerMoeller(1, Pcopy);
  If Tmp = [] Then Error("IdealAndSeparatorsOfPoints: points must be distinct"); EndIf;
  Rescale := [[Indet(J), Denoms[J]*Indet(J)] | J In 1..Len(Denoms)];
  S := Subst(Tmp[2], Rescale);
  For J:=1 To Len(S) Do S[J] := S[J]/Eval(S[J], P[J]); EndFor;
  Tmp[2] := []; -- to avoid wasting memory
  GB := Subst(Tmp[1], Rescale);
  GB := [Monic(F) | F In GB];
  Tmp := []; -- to avoid wasting memory
  I := Ideal(GB);
  I.GBasis := GB;
  Return Record[points := P, ideal := I, separators := S];
EndDefine; -- IdealAndSeparatorsOfPoints



Define IdealAndSeparatorsOfProjectivePoints(PolyR, P)
  If P = [] Then Return Record[points := P, ideal := Ideal(1), separators := []]; EndIf;
  $points.BuchbergerMoellerArgCheck(P, PolyR, "IdealAndSeparatorsOfProjectivePoints");
  Pcopy := P;
  Denoms := $points.CommonDenoms(PolyR, ref Pcopy);
  Tmp := BuchbergerMoeller(0, Pcopy);
  If Tmp = [] Then Error("IdealAndSeparatorsOfProjectivePoints: points must be distinct"); EndIf;
  Rescale := [[Indet(J), Denoms[J]*Indet(J)] | J In 1..Len(Denoms)];
  S := Subst(Tmp[2], Rescale);
  Tmp[2] := []; -- to avoid wasting memory
  GB := Subst(Tmp[1], Rescale);
  GB := [Monic(F) | F In GB];
  Tmp := []; -- to avoid wasting memory
  I := Ideal(GB);
  I.GBasis := GB;
  Return Record[points := P, ideal := I, separators := S];
EndDefine; -- IdealAndSeparatorsOfProjectivePoints



-----------------------------------------------------------------------------
-- Sundry other handy functions more or less related to those above.

Define Interpolate(Points, Values)
  If Len(Points) <> Len(Values) Then Error("Interpolate: different number of points and values"); EndIf;
  If Len(Points) = 0 Then Return 0; EndIf;
  S := $points.SeparatorsOfPoints(Points);
  Return Sum([Values[I]*S[I] | I In 1..Len(Points)]);
EndDefine;


Define GenericPoints(PolyR, N, Opt Range)
  TopLevel ERR;
  GenericPoints_ERR := "First argument must be a polynomial ring";
  If Type(PolyR) <> RING Then error(GenericPoints_ERR); EndIf;
  If not(IsPolyRing(PolyR)) Then error(GenericPoints_ERR); EndIf;
//  N := Cast(ARGV[1], INT);
  P := Characteristic(PolyR);
  If IsDefined(Range) Then
//    Range := Cast(ARGV[2], INT);
    If Range < 1 Then Range := 100; EndIf;
    If P <> 0 And Range >= P Then Range := P-1; EndIf;
  Else
    If P <> 0 Then Range := P-1;  Else Range := 100; EndIf;
  EndIf;
  If P <> 0 Then Lo := 0; Hi := Range; Else Lo := -Range; Hi := Range; EndIf;
  V := NumIndets(PolyR);
  Ans := NewList(N);
//  M := Identity(V);
  M := [[ 0 | J In 1..V] | K In 1..V];
  For J:=1 To NumIndets(PolyR) Do M[J,J] := 1; EndFor;
  For I:=1 To Min(V, N) Do
    Ans[I] := M[I];
  EndFor;
  If N <= V Then Return Ans; EndIf;
  Ans[V+1] := NewList(V, 1);
  For I:=V+2 To N Do
    Ans[I] := [random(Lo, Hi) | J In 1..V];
  EndFor;
  Return Ans;
EndDefine;



EndPackage;


