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 }
}
aldor
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.