f := (1-t)*(1-s)*(1-u)/(1-2*(s+t+u)+3*(s*t+t*u+u*s)-4*s*t*u): F := normal(eval(f, {t=t/s, u=x/t})/s/t): q := denom(F): q1 := q/s/t: disc := factor(discrim(q1, t)): # Same as input to RectangularSystem.mpl. P1 := 2*(s-1)*(3*s^2-6*s+2) + (6*x*s^3-2*s^3-10*x*s^2+s^2-4*x^2*s+10*x*s+3*x^2-4*x) * dx + 2*s*(3*s-2)*(s-1)^2 * ds: P2 := -2*(-19*s^2-9*x+13*s^3+7*s-16*x*s^2+24*x*s) * dx + disc * dx^2: p0 := -eval(P1, {ds=0, dx=0})/coeff(P1, ds, 1): p1 := -coeff(P1, dx, 1)/coeff(P1, ds, 1): q := -coeff(P2, dx, 1)/coeff(P2, dx, 2): Alg := Ore_algebra:-diff_algebra([dx, x], [ds, s]): with(Groebner): MOrd := MonomialOrder(Alg, plex(ds, dx)): GB := Basis([P1, P2], MOrd): Q := phi[0](s) + phi[1](s)*dx: ds_Q := collect(diff(Q, s) + Q*ds, {dx, ds}, distributed); d := 3: add(eta[i] * NormalForm(dx^i, GB, MOrd), i=0..d) - ( phi[1](s) * NormalForm(dx*ds, GB, MOrd) + phi[0](s) * NormalForm(ds, GB, MOrd) + eval(ds_Q, ds = 0) ): zz := collect(%, dx, normal); M := 2: N := 7: eval(zz, { phi[0] = unapply( add(w0[i]*s^i, i=0..M)/(s-1), s ), phi[1] = unapply( add(w1[i]*s^i, i=0..N)/2/(s-1)^2/disc, s ) }): collect(%, dx, numer@normal): collect(%, {dx, s}, distributed): sys := {coeffs(%, {dx, s})}: sol := SolveTools:-Linear(sys, {seq(eta[i], i=0..d), seq(w0[i], i=0..M), seq(w1[i], i=0..N)}): # This complicated-looking final normalisation is session independent. eval([[seq(eta[i], i=0..d)], [seq(w0[i], i=0..M)], [seq(w1[i], i=0..N)]], sol); %[1][d+1]: factor( eval(%%, op(indets(%) minus {x}) = denom(%)) ): sign(%[1][d+1]): SOL := factor( expand(% * %%) ); # phi[0]: phi[0]:= factor( add(SOL[2][i+1]*s^i, i=0..M) ) / (s-1); # phi[1]: gamma_ := add(SOL[3][i+1]*s^i, i=0..N); degree(gamma_, x), degree(gamma_, s); phi[1] := gamma_ / (2*(s-1)^2*disc): # Check that P - ds Q is 0 modulo the ideal I = . P := add(SOL[1][i]*dx^(i-1), i=1..nops(SOL[1])); Q := phi[0] + phi[1]*dx: Groebner:-Reduce(P - Ore_algebra:-skew_product(ds, Q, Alg), [P1, P2], MOrd);