login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit
 Topics FrontPage SandBox SandBox Speed <-- You are here. 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.

 Subject:   Be Bold !! ( 14 subscribers )