Edit detail for SandBox Speed revision 1 of 1

changed:
-
\begin{aldor}
#include "axiom.as"

INT ==> Integer; 
NNI ==> NonNegativeInteger;
P ==> UnivariatePolynomial(x ,F);
FP ==> Fraction P ;
VFP ==> Vector FP;
LP ==> List P;
LFP ==> List FP;
MFP ==> Matrix FP; 

PREC ==> SparseTable( INT , P, 0) ;
ListHelp(R:Ring): with {
        
        clean: (R, List R) -> List R ;
        ++ clean(z, lst) [*,z,*,*,z,z,z,z] -> [*,z,*,*] 
        cleanZero:List R -> List R ;
        ++ equivalent to clean(0,lst)

}  == add {
 
        -- Get rid of the trailing z
        clean(z:R, lst:List R): List R == {
                null(lst) => lst ;
                headl: R := first lst ;
                -- Not "z" 
                headl ~= z => cons(headl, clean(z, rest lst )) ;
                -- "z"
                restList:List R := clean(z ,rest lst );
                null restList => restList ;
                cons( headl, restList) ;
        }
        
        cleanZero(lst: List R): List R == {
                clean(0$R, lst) -- $
        }
}   
UDFS(x:Symbol,F:Field): with {   

    +:(%,%) -> % ;
    *:(%,%) -> % ;         

   -- Should I create a new intead ? 
   coerce: LFP -> % ;
   coerce: % -> LFP ;
   coerce: % -> OutputForm ;        

    linearise: % -> % ;
    linearise: % -> List P ;


    order: % -> INT ;
   degree: % -> INT ;

    homogenous?: % -> Boolean ;    
      

    diffComb: (% ,INT) -> List VFP ;
      
    precursive: % ->  INT ;

        
   -- Should not be plublic just for testing    
    pochhammer: (P ,INT)-> P ; -- could be generalize
           
           findSolutionDESpace: (Matrix FP, INT) -> List VFP ;


}  == add {
    Rep ==> LFP; -- List Fractional UnivaratePolynomial

   import from INT, NNI ;
   import from F, P, LFP, VFP, List VFP, List P;
   import from PrintPackage ;
   import from Vector F;
   import from Vector VFP ;
   import from MFP;
   import from List MFP ;
   import from PREC ;
   import from Boolean ; -- Needed ???

   -- This function need to be somehow robust
   coerce (f:LFP):% == {
      cleanF: LFP := cleanZero(f)$ListHelp(FP) ;
      null cleanF => error "List need not to be empty" ;           
      per cleanF; 
   }
           
   coerce (f:%):LFP == rep f;

   coerce (f:%):OutputForm == {
      a1 := rep f :: OutputForm ;
      a2 := message("D ")$OutputForm;
      hconcat(a2,a1) }

   order (f:%):INT  == {
      cleanF: LFP := cleanZero( rep f )$ListHelp( FP) ; --$
      (#(rep f))::INT - 2   }   
   degree (f:%):INT  == {
      a1 := [degree( i ) for i in linearise( f) ]@List NonNegativeInteger ;
      print (a1::OutputForm) ;
      reduce( max, a1,0) :: Integer }
   homogenous?( f:% ): Boolean == {
      frep := rep f;
      #frep = 0 => error "No definied";
      frep.1 = 0 => true ;
      false ;
   }
   -- As make a diff eq homegenous increase the order by one 
   -- I wildly speculte the following code
   orderNonHomogenous(f:%):INT == {
      res: INT := order f ;
      homogenous? f => res + 1 ;
      res ;
   }
   orderNH ==> orderNonHomogenous ;
       

   linearise (f:%):%  == {
      frep := rep f ; 
      denomLcm := lcm [ denom i for i in frep ] ;
      res:LFP := [(denomLcm * i)  for i in frep ] ;
      per res ;
   }  

   linearise (f:%):List P  == {
      res: LFP := rep( linearise f ) ;
      [numer i for i in res ] ;
   }
   nextDiffComb( fK:VFP, fD:VFP):VFP  == { 
      sizeVec: INT := #fK::NNI::INT ;
      -- Leibniz rule     
      fKp1:VFP := map( (i:FP):FP +-> D( i)$FP, fK) ;
      --we start at 3 because degre c and f are never affected
      for i in 3..sizeVec repeat { 
            fKp1.i := fKp1.i + fK.(i-1) } 

      -- we substitute the last degree       
      fKp1 := fKp1 + fD*fK.sizeVec ;
   }
   -- We can really be smarter about it and having a constructive function 
   -- give a liste of f(n)= q_n +sum d-1,0 f(i) up to deg
   diffComb( f:% , deg:INT ):List VFP == {
      frep := rep f;
      -- we use vector representation for speed
      ordf:INT := order( f) ;
      sizef := ordf +2 ; 
      sizeVec := sizef - 1 ; -- The last degree is implicit
      deg < ordf => 
          error "Rec has higher degree than deg asked" ; 
      lastf:FP  := frep.last ;
      fD:VFP := vector [ (- frep.i / lastf)::FP  for i in 1..sizeVec ] ;
      -- TD : Use vector
      res:List VFP := reverse [ [(if i = n then 1 else 0) for i in 1..sizeVec] $VFP for n in 2..(sizef-1) ] ;
      res := cons( fD, res) ;
      sizeRes:INT := #res :: INT ;
      -- Implement Leibniz rule, + substitution of max power 
      -- I guess there should be smarter think todo as use ODPOL ...
      new:VFP := fD ;
      if deg >= sizeRes  then {
            for w in sizeRes..deg repeat {
                  prev:VFP := new ;
                  new:VFP := nextDiffComb( prev , fD ) ; 
                  res := cons( new,res)  }}
      print( res::OutputForm) ;
      reverse res 
   }
    
   findSolutionDESpace( m:Matrix FP, min:INT):List VFP == {
      -- We do the heavy lifting only once
      recFindSolution(rowEchelon m, min ) ;
   }

   recFindSolution( m:Matrix FP, min:INT):List VFP == {
      print( min :: OutputForm );
      min < 1 => error "Something wrong except largeur min" ;
      -- m is the largest possible we get its size
      rowSize:INT  := nrows(m)::INT ;  
      colMaxSize:INT := ncols(m)::INT;
      -- Maybe I should just give an empty list ?
      min > colMaxSize => error "We found nothing" ;

      -- We create the submatrix
      subM:Matrix FP :=subMatrix(m,1,rowSize, 1, min) ;
      if( nullity(subM) = 0) then {
              recFindSolution(m, min +1 ); }
      else {
              nullSpace subM ; }
   }


   -- we contruct the desired matrix
   (f:%) + (g:%):% == {
       ordh:INT := orderNH(f) + orderNH(g) ;
       fdiff:List LFP := [entries i for i in diffComb( f, ordh )];
       gdiff:List LFP := [entries i for i in diffComb( g, ordh )];
       mfdiff: Matrix FP := transpose matrix fdiff;
       mgdiff: Matrix FP := transpose matrix gdiff;

       m: Matrix FP := new((ordh+1)::NNI ,(ordh+2)::NNI ,0 );
       m.(1,1) := 1;

       -- the constant row is the only one with mix from f and g
       -- We need a list to recreate a matrix
       -- coercing to matrix will creave a colum matrix
       cstRow: LFP := entries(row(mfdiff ,1) + row(mgdiff ,1));
       setsubMatrix!(m,1,2, matrix( [cstRow] ));
       -- now we extrat the constant row
       submf: Matrix FP := subMatrix( mfdiff, 2 ,nrows(mfdiff)::INT, 1,ncols(mfdiff)::INT);
       submg: Matrix FP := subMatrix( mgdiff, 2 ,nrows(mgdiff)::INT, 1,ncols(mgdiff)::INT);
       -- two different setSubmatrix should be faster than the concat
       setsubMatrix!(m,2,2, vertConcat(submg,submf));
       solm: List VFP := nullSpace m;
       print( m::OutputForm ) ;
       print( solm::OutputForm ) ;

       sol2: List VFP  := findSolutionDESpace(m, max( orderNH f, orderNH g ));
       print( sol2::OutputForm ) ; 

   -- DUMMY SOLUTION
       f}

   -- I have a huge pb with the order !
   -- Look like I need a larger degree if it not homegenous
   (f:%) * (g:%):% == {
      ordh: INT := (orderNH(f) )*(orderNH(g)+1) ;

       df: Vector VFP := vector diffComb(f, ordh) ;
      dg: Vector VFP := vector diffComb(g, ordh) ;

      -- Thats a huge one liner
      -- I should split it
      --we use outer product for computing all the possibility( fg f'g f''g ect easely
      lDfg: List List FP := [parts reduce(+,[ (binomial(k,i)*outerProduct(df.(i+1), dg.(k-i+1))) for i:INT in 0..k]) for k:INT in 0..ordh ] ;

      mpDfg: Matrix FP := transpose matrix lDfg;
      print(mpDfg ::OutputForm ) ;
      -- the constant (non homegenous) term
      cTerm: Matrix FP := new(nrows(mpDfg),1,0) ;
      cTerm.(1,1) := 1 ;
      mDfg: Matrix FP := horizConcat(cTerm,mpDfg) ;
      res: List VFP := nullSpace mDfg ; -- Field are Integral Domain 
      print( res::OutputForm ) ;
      -- per (nullSpace mDfg).1 
      f}  
        
   pochhammer(n:P, k:INT):P == {
      k < 0 => error "not defintied for negatif power";
      k = 0 => 1 ;
      res :P := n ;
      for i in 1..(k-1) repeat {
              res := res*(n+i::P) ;
      }
      res;
   }

   precursive( f:%): INT == {
      -- should be homogenized before hand
      not (homogenous? f) => error "Expext homogenous only" ; 
      -- We are going to put our future p-rec there
      -- I hope the cost is not too bad
      precTable : SparseTable( INT , P, 0) := table() ;
      -- we should also homegenize
      linf: List P := linearise (f) ;
      -- we are going to use nested for loops so we have the index
      -- Maybe pattern matching will have be of interest with more advanced structure
      -- first iteration is for the degree 0 
      -- instead of n it will be x that could be fix
     
      for j in 0..( order f) repeat {
              -- We put each polynome in a vector only way to iterate I found
              polyj: P := linf.(j+2) ; -- term c and compensation
   --                      degPoly:NNI := degree polyj ;
   --                      vpolyj: Vector F := vectorise( polyj, degPoly +1 );
   --                      print( vpolyj::OutputForm) ;
              -- we iterate on each degree there is a difference betwen the degre and the position in the vector
              for monomialPolyj in monomials polyj  repeat {
                      -- The degree of the monomial
                      k:INT := degree monomialPolyj :: INT ;
                      coeffK:F := content monomialPolyj ;
                      -- we use x instead of p
                      -- x := monomial(1,1)
                      precTable.(j - k) := coeffK * pochhammer( monomial(1,1)$P + (1-k)::P,j) ;
              } 
      }
      print( precTable :: OutputForm ); 
      12 }

}
   

\end{aldor}

#include "axiom.as"

INT ==> Integer; 
NNI ==> NonNegativeInteger;
P ==> UnivariatePolynomial(x ,F);
FP ==> Fraction P ;
VFP ==> Vector FP;
LP ==> List P;
LFP ==> List FP;
MFP ==> Matrix FP; 

PREC ==> SparseTable( INT , P, 0) ;
ListHelp(R:Ring): with {

        clean: (R, List R) -> List R ;
        ++ clean(z, lst) [*,z,*,*,z,z,z,z] -> [*,z,*,*] 
        cleanZero:List R -> List R ;
        ++ equivalent to clean(0,lst)

}  == add {

        -- Get rid of the trailing z
        clean(z:R, lst:List R): List R == {
                null(lst) => lst ;
                headl: R := first lst ;
                -- Not "z" 
                headl ~= z => cons(headl, clean(z, rest lst )) ;
                -- "z"
                restList:List R := clean(z ,rest lst );
                null restList => restList ;
                cons( headl, restList) ;
        }

        cleanZero(lst: List R): List R == {
                clean(0$R, lst) -- $
        }
}   
UDFS(x:Symbol,F:Field): with {   

    +:(%,%) -> % ;
    *:(%,%) -> % ;         

   -- Should I create a new intead ? 
   coerce: LFP -> % ;
   coerce: % -> LFP ;
   coerce: % -> OutputForm ;        

    linearise: % -> % ;
    linearise: % -> List P ;

    order: % -> INT ;
   degree: % -> INT ;

    homogenous?: % -> Boolean ;    

    diffComb: (% ,INT) -> List VFP ;

    precursive: % ->  INT ;

   -- Should not be plublic just for testing    
    pochhammer: (P ,INT)-> P ; -- could be generalize

           findSolutionDESpace: (Matrix FP, INT) -> List VFP ;

}  == add {
    Rep ==> LFP; -- List Fractional UnivaratePolynomial

   import from INT, NNI ;
   import from F, P, LFP, VFP, List VFP, List P;
   import from PrintPackage ;
   import from Vector F;
   import from Vector VFP ;
   import from MFP;
   import from List MFP ;
   import from PREC ;
   import from Boolean ; -- Needed ???

   -- This function need to be somehow robust
   coerce (f:LFP):% == {
      cleanF: LFP := cleanZero(f)$ListHelp(FP) ;
      null cleanF => error "List need not to be empty" ;           
      per cleanF; 
   }

   coerce (f:%):LFP == rep f;

   coerce (f:%):OutputForm == {
      a1 := rep f :: OutputForm ;
      a2 := message("D ")$OutputForm;
      hconcat(a2,a1) }

   order (f:%):INT  == {
      cleanF: LFP := cleanZero( rep f )$ListHelp( FP) ; --$
      (#(rep f))::INT - 2   }   
   degree (f:%):INT  == {
      a1 := [degree( i ) for i in linearise( f) ]@List NonNegativeInteger ;
      print (a1::OutputForm) ;
      reduce( max, a1,0) :: Integer }
   homogenous?( f:% ): Boolean == {
      frep := rep f;
      #frep = 0 => error "No definied";
      frep.1 = 0 => true ;
      false ;
   }
   -- As make a diff eq homegenous increase the order by one 
   -- I wildly speculte the following code
   orderNonHomogenous(f:%):INT == {
      res: INT := order f ;
      homogenous? f => res + 1 ;
      res ;
   }
   orderNH ==> orderNonHomogenous ;

   linearise (f:%):%  == {
      frep := rep f ; 
      denomLcm := lcm [ denom i for i in frep ] ;
      res:LFP := [(denomLcm * i)  for i in frep ] ;
      per res ;
   }  

   linearise (f:%):List P  == {
      res: LFP := rep( linearise f ) ;
      [numer i for i in res ] ;
   }
   nextDiffComb( fK:VFP, fD:VFP):VFP  == { 
      sizeVec: INT := #fK::NNI::INT ;
      -- Leibniz rule     
      fKp1:VFP := map( (i:FP):FP +-> D( i)$FP, fK) ;
      --we start at 3 because degre c and f are never affected
      for i in 3..sizeVec repeat { 
            fKp1.i := fKp1.i + fK.(i-1) } 

      -- we substitute the last degree       
      fKp1 := fKp1 + fD*fK.sizeVec ;
   }
   -- We can really be smarter about it and having a constructive function 
   -- give a liste of f(n)= q_n +sum d-1,0 f(i) up to deg
   diffComb( f:% , deg:INT ):List VFP == {
      frep := rep f;
      -- we use vector representation for speed
      ordf:INT := order( f) ;
      sizef := ordf +2 ; 
      sizeVec := sizef - 1 ; -- The last degree is implicit
      deg < ordf => 
          error "Rec has higher degree than deg asked" ; 
      lastf:FP  := frep.last ;
      fD:VFP := vector [ (- frep.i / lastf)::FP  for i in 1..sizeVec ] ;
      -- TD : Use vector
      res:List VFP := reverse [ [(if i = n then 1 else 0) for i in 1..sizeVec] $VFP for n in 2..(sizef-1) ] ;
      res := cons( fD, res) ;
      sizeRes:INT := #res :: INT ;
      -- Implement Leibniz rule, + substitution of max power 
      -- I guess there should be smarter think todo as use ODPOL ...
      new:VFP := fD ;
      if deg >= sizeRes  then {
            for w in sizeRes..deg repeat {
                  prev:VFP := new ;
                  new:VFP := nextDiffComb( prev , fD ) ; 
                  res := cons( new,res)  }}
      print( res::OutputForm) ;
      reverse res 
   }

   findSolutionDESpace( m:Matrix FP, min:INT):List VFP == {
      -- We do the heavy lifting only once
      recFindSolution(rowEchelon m, min ) ;
   }

   recFindSolution( m:Matrix FP, min:INT):List VFP == {
      print( min :: OutputForm );
      min < 1 => error "Something wrong except largeur min" ;
      -- m is the largest possible we get its size
      rowSize:INT  := nrows(m)::INT ;  
      colMaxSize:INT := ncols(m)::INT;
      -- Maybe I should just give an empty list ?
      min > colMaxSize => error "We found nothing" ;

      -- We create the submatrix
      subM:Matrix FP :=subMatrix(m,1,rowSize, 1, min) ;
      if( nullity(subM) = 0) then {
              recFindSolution(m, min +1 ); }
      else {
              nullSpace subM ; }
   }

   -- we contruct the desired matrix
   (f:%) + (g:%):% == {
       ordh:INT := orderNH(f) + orderNH(g) ;
       fdiff:List LFP := [entries i for i in diffComb( f, ordh )];
       gdiff:List LFP := [entries i for i in diffComb( g, ordh )];
       mfdiff: Matrix FP := transpose matrix fdiff;
       mgdiff: Matrix FP := transpose matrix gdiff;

       m: Matrix FP := new((ordh+1)::NNI ,(ordh+2)::NNI ,0 );
       m.(1,1) := 1;

       -- the constant row is the only one with mix from f and g
       -- We need a list to recreate a matrix
       -- coercing to matrix will creave a colum matrix
       cstRow: LFP := entries(row(mfdiff ,1) + row(mgdiff ,1));
       setsubMatrix!(m,1,2, matrix( [cstRow] ));
       -- now we extrat the constant row
       submf: Matrix FP := subMatrix( mfdiff, 2 ,nrows(mfdiff)::INT, 1,ncols(mfdiff)::INT);
       submg: Matrix FP := subMatrix( mgdiff, 2 ,nrows(mgdiff)::INT, 1,ncols(mgdiff)::INT);
       -- two different setSubmatrix should be faster than the concat
       setsubMatrix!(m,2,2, vertConcat(submg,submf));
       solm: List VFP := nullSpace m;
       print( m::OutputForm ) ;
       print( solm::OutputForm ) ;

       sol2: List VFP  := findSolutionDESpace(m, max( orderNH f, orderNH g ));
       print( sol2::OutputForm ) ; 

   -- DUMMY SOLUTION
       f}

   -- I have a huge pb with the order !
   -- Look like I need a larger degree if it not homegenous
   (f:%) * (g:%):% == {
      ordh: INT := (orderNH(f) )*(orderNH(g)+1) ;

       df: Vector VFP := vector diffComb(f, ordh) ;
      dg: Vector VFP := vector diffComb(g, ordh) ;

      -- Thats a huge one liner
      -- I should split it
      --we use outer product for computing all the possibility( fg f'g f''g ect easely
      lDfg: List List FP := [parts reduce(+,[ (binomial(k,i)*outerProduct(df.(i+1), dg.(k-i+1))) for i:INT in 0..k]) for k:INT in 0..ordh ] ;

      mpDfg: Matrix FP := transpose matrix lDfg;
      print(mpDfg ::OutputForm ) ;
      -- the constant (non homegenous) term
      cTerm: Matrix FP := new(nrows(mpDfg),1,0) ;
      cTerm.(1,1) := 1 ;
      mDfg: Matrix FP := horizConcat(cTerm,mpDfg) ;
      res: List VFP := nullSpace mDfg ; -- Field are Integral Domain 
      print( res::OutputForm ) ;
      -- per (nullSpace mDfg).1 
      f}  

   pochhammer(n:P, k:INT):P == {
      k < 0 => error "not defintied for negatif power";
      k = 0 => 1 ;
      res :P := n ;
      for i in 1..(k-1) repeat {
              res := res*(n+i::P) ;
      }
      res;
   }

   precursive( f:%): INT == {
      -- should be homogenized before hand
      not (homogenous? f) => error "Expext homogenous only" ; 
      -- We are going to put our future p-rec there
      -- I hope the cost is not too bad
      precTable : SparseTable( INT , P, 0) := table() ;
      -- we should also homegenize
      linf: List P := linearise (f) ;
      -- we are going to use nested for loops so we have the index
      -- Maybe pattern matching will have be of interest with more advanced structure
      -- first iteration is for the degree 0 
      -- instead of n it will be x that could be fix

      for j in 0..( order f) repeat {
              -- We put each polynome in a vector only way to iterate I found
              polyj: P := linf.(j+2) ; -- term c and compensation
   --                      degPoly:NNI := degree polyj ;
   --                      vpolyj: Vector F := vectorise( polyj, degPoly +1 );
   --                      print( vpolyj::OutputForm) ;
              -- we iterate on each degree there is a difference betwen the degre and the position in the vector
              for monomialPolyj in monomials polyj  repeat {
                      -- The degree of the monomial
                      k:INT := degree monomialPolyj :: INT ;
                      coeffK:F := content monomialPolyj ;
                      -- we use x instead of p
                      -- x := monomial(1,1)
                      precTable.(j - k) := coeffK * pochhammer( monomial(1,1)$P + (1-k)::P,j) ;
              } 
      }
      print( precTable :: OutputForm ); 
      12 }

}

   Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/7884349318222726611-25px001.as
      using Aldor compiler and options 
-O -Fasy -Fao -Flsp -lfricas -Mno-ALDOR_W_WillObsolete -DFriCAS -Y $FRICAS/algebra -I $FRICAS/algebra
      Use the system command )set compiler args to change these 
      options.
   The )library system command was not called after compilation.