# Sequences from OEIS (http://www.research.att.com/~njas/sequences/) last downloaded 7/18/10 with(LinearAlgebra, Determinant, LinearSolve): # Determinant and LinearSolve are used in findgaugeandterm with(LRETools, hypergeomsols): read "/Users/heba/Desktop/Implementations/Giles_Levy_implementation/code/findgauge_v2.1.txt": # G:=proc(Ra,Rb,u,v,n) for GT: Ra in terms of u, Rb in terms of v read "/Users/heba/Desktop/Implementations/Giles_Levy_implementation/code/curvp_v1.3.txt": # This calls multrecs_v1.3 read "/Users/heba/Desktop/Implementations/Giles_Levy_implementation/code/Jsolver4.txt": read "/Users/heba/Desktop/Implementations/Giles_Levy_implementation/code/groups_with_pcurv_and_offset_v2.6.txt": # file with latest Sloane sequences _All_2nd_irr := `/M/m/_All_2nd_irr`: # only gets read if needed # multrecs:=proc(L1,L2,u,n),input 2 lists of coeffs, each minus u(n+ord) (output in terms of the input u and n) for j to 7 do # Just in case _u or _n is assigned # and we want _v global (for Liouv case) and _twoFone # _t and _T are from J's code if not type(eval([_u, _v, _n, _tau, _twoFone, _t, _T][j]), symbol) then # print(cat("Warning: ", ['_u', '_v', '_n', '_tau', '_twoFone', '_t', '_T'][j], " has been unassigned")); assign(['_u', '_v', '_n', '_tau', '_twoFone', '_t', '_T'][j], ['_u', '_v', '_n', '_tau', '_twoFone', '_t', '_T'][j]) fi od; # we'll repeat this every time findrel is called (see below) # ***************************************************************************************************************************************************************** findrel := proc(input, rname, ics) # input rec eqn, name (e.g. u(n)), and (optional) a set of sequential IC's (e.g. {u(12)=3, u(13)=55, ...}) # IMPORTANT: This procedure assumes that eqn is already of form u(n+2)...u(n) not something like u(n)...u(n-2) local numprimes, i, inits, opts, inp, invar, ind, k, L, two, LL, arbitrary, Det, inptau, result, j, sub, unsub, temp1, temp2, temp3, temp4, start, smaller, answer; global _u, _v, _n, _tau, _twoFone, _info_for_findrel, _t, _T, results_with_pcurv; for j to 7 do # Just in case _u or _n is assigned # and we want _v global (for Liouv case) and _twoFone # _t and _T are from J's code if not type(eval([_u, _v, _n, _tau, _info_for_findrel, _twoFone, _t, _T][j]), symbol) then # print(cat("Warning: ", ['_u', '_v', '_n', '_tau', '_twoFone', '_t', '_T'][j], " has been unassigned")); assign(['_u', '_v', '_n', '_tau', '_info_for_findrel', '_twoFone', '_t', '_T'][j], ['_u', '_v', '_n', '_tau', '_info_for_findrel', '_twoFone', '_t', '_T'][j]) fi od; numprimes := nops(_Sloane_reps[1][2]); # Tells us for how many primes we check p-curvature to get invar below opts := [1$4,0$6]; # defaults (1-4 are Liouv, 2F1, special fns, and Sloane) if nargs = 3 then # options are in ics so if there are options but no inits := ics; if membertype(list, inits, 'opts') then # i.e. options were input (options are a list within ics) opts, inits := inits[opts], subsop(opts=NULL, inits); opts := [op(opts), seq(1, j=nops(opts)+1..10)] fi; inits := subs(op(0, rname) = _u, inits) else inits := {} fi; inp := subs({op(0, rname) = _u, op(1, rname) = _n}, input); if hastype(inp, equation) then inp := lhs(inp)-rhs(inp) fi; if not (has(inp, _u(_n+2)) and has(inp, _u(_n))) then return NULL fi; if coeff(inp, _u(_n+1)) = 0 then return NULL fi; j := hypergeomsols(inp, _u(_n), {}); if j <> 0 then j := convert( eval(j, {_u = op(0, rname), _n = op(1, rname)}), string); return `input has hypergeometric solution: ` || j fi; L := []; LL := NULL; # LL will hold IC's in the form [offset, [sequence]], e.g. for u(3)=4, u(1)=5, u(2)=2 we get LL = [1, [5,2,4]] if nops(inits) > 0 then # one IC is okay if the previous one is arbitrary, i.e. this one determines the sequence for j in inits do temp1 := isolate(j, _u); L := [op(L), [op(1, lhs(temp1)), rhs(temp1)]] od; smaller := (L1,L2) -> evalb(L1[1] <= L2[1]); L := sort(L, smaller); LL := [L[1][1], [seq(L[j][2], j=1..nops(L))]]; if L[-1][1] - L[1][1] + 1 <> nops(L) then userinfo(5, findrel, `Must input all initial conditions in range (e.g. can't input _u(2), _u(0) without _u(1))`); LL := NULL fi; fi; # this next line must be done before we divide through by coeff u(n+2) lest we lose some "content," e.g. an n multiplying all terms if LL <> NULL then arbitrary := find_arbitrary(inp, LL); else arbitrary := NULL fi; inp := collect(inp/coeff(inp, _u(_n+2)), _u); Det := NormalizeRat(coeff(inp, _u(_n))/coeff(inp, _u(_n+2)), _n); inptau := subs(seq(_u(_n+j)=_tau^j, j=0..2), inp); # use this in 2F1 check, Liouvillian check, and in call to J's code # ******** Liouvillian check ******** temp1 := `if`(opts[1] = 1, FindLiouvillian(inptau, _tau, _n), NULL); if temp1 <> NULL then indets(rhs(temp1[2])) minus {n}; temp1 := [add(coeff(temp1[1], _tau, j)*_v(_n+j), j=0..2)=0, subs(op(0, %[1]) = _v, temp1[2])]; if LL <> NULL then _info_for_findrel[1] := findtwo([rhs(temp1[2]), temp1[1]], [LL[1], inp, LL[2]], arbitrary); # I'll be using this info in external test file temp1 := `if`(info_for_findrel[1] = NULL, NULL, [op(temp1), _info_for_findrel[1][2]]); # this case:findtwo returns offset, 2 from new seq, old seq _info_for_findrel[1] := `if`(temp1 <> NULL, subs(_n = op(1, rname), [_info_for_findrel[1]]), NULL) fi; if temp1 <> NULL then temp1 := subs(_u = op(0, rname), _n = op(1, rname), temp1); return op(temp1) fi fi; # ******** 2F1 check ******** answer := `if`(opts[2] = 1, Find2F1(inptau, _tau, _n), NULL); if answer <> NULL and answer <> "FAIL" then temp3 := subs(seq(_twoFone(_n+i) = subs(_n=_n+i, rhs(answer[2])), i=0..1), answer[1]); temp3 := optimize(temp3, answer[3]); if answer[3][3] = 1 then # if c <> 1 then we have two lin ind solutions (at this point: if c in Z then c = 1) unprotect(_twoFone); _twoFone := '_twoFone'; if LL <> NULL then temp2 := result_ICs([subs(_c=1, temp3)], LL, arbitrary); if temp2 <> NULL then return rname = subs(_n=op(1, rname), temp2[1]), op(1, rname) >= temp2[2] fi else WARNING("not all solutions found"); # lower case 'n' b/c WARNING("message") prints: Warning, message return rname = subs(_n=op(1, rname), temp3) fi else unprotect(_twoFone); _twoFone := '_twoFone'; temp4 := subs(seq(_twoFone(_n+i) = subs(_n=_n+i, rhs(answer[4])), i=0..1), answer[1]); temp4 := optimize(temp4, answer[5]); if LL <> NULL then temp2 := result_ICs([subs(_c=1, temp3), subs(_c=1, temp4)], LL, arbitrary); if temp2 <> NULL then return rname = subs(_n=op(1, rname), temp2[1]), op(1, rname) >= temp2[2] fi else return rname = subs(_n=op(1, rname), _c=_c[1], temp3) + subs(_n=op(1, rname), _c=_c[2], temp4) fi fi fi; # *********** Call to J's code begins *********** inptau := collect(inptau/coeff(inptau, _tau, 2), _tau); # try temp1 := `if`(opts[3] = 1, solver(inptau), NULL); # catch: # temp1 := NULL # finally if temp1 = 'FAIL' then temp1 := NULL fi; # end: result := NULL; temp2 := NULL; if temp1 <> NULL then if LL <> NULL then temp2 := result_ICs(temp1, LL, arbitrary); if temp2 <> NULL then result := subs(_n=op(1, rname), temp2[1]) else return NULL fi fi; return rname = `if`(result = NULL, add(_c[j]*temp1[j], j=1..nops(temp1)), result), `if`(result <> NULL, op(1, rname) >= temp2[2], NULL) fi; # *********** Call to Database *********** if opts[4] = 1 then ind := op(2,ArrayDims(_Sloane_reps)) else # else we skip database work (and may as well skip p-curvature computations ind, numprimes := 0,0 fi; invar := []; for j from 1 to numprimes do try invar := [op(invar), curvp(inp, ithprime(j), _u, _n)] catch: invar := [ op(invar), [] ] end od; for k to ind do i := _Sloane_reps[k]; if i[2] <> invar then next fi; two := 1; if nops(i[3]) > 1 and nops(inits) > 1 then two := 2 fi; # does this representative (i=_Sloane_reps[k]) have > 1 solution/is input a sequence or an equation answer := [findgaugeandterm([inp, _u(_n), LL ], [ seq( ' i[5][j], _u(_n), [ i[6][j], i[7][j] ], i[4][j] ', j=1..two) ])]; if answer <> [] and two = 1 and nops(i[3]) > 1 then # No IC's but there are 2 reps temp1 := findgaugeandterm([inp, _u(_n), LL ], [ seq( ' i[5][j], _u(_n), [ i[6][j], i[7][j] ], i[4][j] ', j=2..2) ]); answer := [subs(_c=_c[1], op(answer)), subs(_c=_c[2], temp1)] fi; if answer <> [] then for j in {seq(op(0, j), j in indets(answer) minus {_n})} do unprotect(j); assign(j, subs(_n=op(1, rname), eval(j))); protect(j) od; if nops(i[3]) = 1 and LL = NULL then # There is only 1 rep and no ICs were input (if ICs were input then there is a match -> 1 solution) WARNING("not all solutions found") # lower case 'n' b/c WARNING("message") prints: Warning, message fi; temp1 := answer[-1]; if type(opts[5], symbol) then # i.e. want info for other entries in the Sloane grouping if type(_All_2nd_irr, symbol) then read _All_2nd_irr # keep from reading more than once (read findrel_v2* resets it to symbol) fi; assign(opts[5], [[`Related (by recurrence, not necessarily by sequence) Sloane entries are given as:`,\ _All_2nd_irr[0]], seq(_All_2nd_irr[j], j in i[1][1])]); fi; if type(temp1, list) then return rname = add(subs(_n=op(1, rname), answer[j]), j=1..nops(answer)-1), op(answer[-1]) else return rname = add(subs(_n=op(1, rname), answer[j]), j=1..nops(answer)) fi fi od; NULL end: # ***************************************************************************************************************************************************************** result_ICs := proc(temp1, LL, arbitrary) # temp1 is a list of two solutions (w/o arb constants and 2nd can be = 0), LL is ICs in form [offset, [seq]] local start, L, temp2, j, temp3, result, temp4, k, fns, sub, x, unsub, co; start, L, temp2, j, temp3, result := LL[1], LL[2], 0, 1, [], false; # start, L are the offset and the seq, temp2 is num of successful, j is index, temp3 is list while temp2 < 2 and j <= nops(L) do # need two terms to be sure even though nops(temp1) (= # undetermined constants) may be < 2 if not j in arbitrary then try # This (temp4 below) doesn't just find division by 0 but inconsistencies as well, e.g. when only 1st IC doesn't match found result temp4 := simplify(L[j] = subs(_n=j-1+start, add(_c[k]*temp1[k], k=1..nops(temp1)))); if has(temp4, [seq(_c[k], k=1..nops(temp1))]) then temp3 := [op(temp3), temp4]; temp2 := temp2 + 1 fi catch: end fi; j := j + 1; if j < nops(L) and L[j-1] = 0 then while L[j] = 0 and j < nops(L) do j := j + 1 od fi od; if temp2 > 1 then temp4, temp2 := {op(indets(temp1) minus {_n})}, []; fns := [hypergeom, LegendreP, LegendreQ, GegenbauerC, WhittakerM, WhittakerW, BesselJ, BesselY, BesselI, BesselK, JacobiP, LaguerreL]; for k to nops(temp4) do if hasfun(temp4[k], fns) then temp2 := [op(temp2), temp4[k]] fi od; sub := [seq(temp2[k] = x[k], k = 1..nops(temp2))]; unsub := [seq(x[k] = temp2[k], k = 1..nops(temp2))]; solve(temp3, [seq(_c[k], k=1..nops(temp1))]); if % <> [] then temp3 := add(simplify(rhs(%[1][k]))*temp1[k], k=1..nops(temp1)); temp3 := subs(sub, temp3); temp3 := collect(primpart(temp3, [seq(x[k], k=1..nops(temp2))], 'co'), [seq(x[k], k=1..nops(temp2))]); temp3 := collect(subs(unsub, temp3), temp2, factor); result := true fi fi; `if`(result, [co*temp3, start+j-3], NULL) # start=offset, j starts at 1 and increases at least 2 end: # ***************************************************************************************************************************************************************** find_arbitrary := proc(rec, LL) # arbitrary union off-limits: all j in Z s.t. 0 <= j < start local start, L, j, arbitrary, eqn, coeffroots, co,k; # Important to send in the (original input) rec eqn unaltered, i.e. no primpart or normalizing after dividing by coeff u(n+2). # If a term of seq is never present in eqn then it's arbitrary, e.g. _u(15) in (_n-15)*_u(_n) + (_n-14)*_u(_n+1) + (_n-13)*_u(_n+2). # Keep in mind, also, that for the first and second term of a sequence we will need fewer conditions (one and two, respectively). # It is important to use this in subroutines, e.g. if u(k) is arbitrary (but not u(k+1)) then u(k+1) defines the seq and we will need # u(k+1) and u(k+2) (if not arbitrary) to determine linear dependence. # *** find_arbitrary finds not just arbitrary, but terms that make a subsequent term one of the defining terms of the sequence, # e.g. if u(0) = u(1) = 0 and this is not the zero sequence then u(0), u(1) would be included in the result, # so any zeros in a pair in {0, arbitrary} are included at the beginning and after >= 2 arbitrary terms (i.e. new beginnings of the seq). start, L := LL[1], LL[2]; arbitrary := {seq(j, j=0..start-1)}; eqn := collect( primpart(rec,[_u(_n),_u(_n+1),_u(_n+2)], 'co'), _u); rec,eqn,co,LL; coeffroots := [seq( select(type, {solve( co*coeff(eqn, _u(_n+j) ), _n)}, 'integer'), j=0..2)]; if start in coeffroots[1] then arbitrary := arbitrary union {start} fi; if start+1 in coeffroots[1] and start in coeffroots[2] then arbitrary := arbitrary union {start + 1} fi; for j in {seq(op(coeffroots[k]), k=1..3)} do if j > start + 1 and j in coeffroots[1] and j-1 in coeffroots[2] and j-2 in coeffroots[3] then arbitrary := arbitrary union {j} fi od; for j to nops(L) - 1 do if (j-1+start in arbitrary or L[j] = 0) and (j+start in arbitrary or L[j+1] = 0) then arbitrary := arbitrary union {j - 1 + start, j + start} fi od; # `if`(arbitrary={},NULL,arbitrary) arbitrary end: # ***************************************************************************************************************************************************************** FindLiouvillian := proc(L, tau, n) local a0,a1,a2,v,detL,sa0,sa1,sa2,c3,c2,c1,c0,v1,v2,v3,LT,R,eqn13, isz,d_values,g,sub2,sub5,A,sg,a,b,den,Lhat,i; a0, a1, a2 := seq(coeff(L, tau, i), i=0..2); if a0=0 or degree(L,tau)<>2 then error "not of order 2" elif a1=0 then return L, _u(_n) = _v(_n) # Step 1 in the paper fi; # Step 2. sub2 := _u(_n+2) = -1/a2 * (a0 * _u(_n) + a1 * _u(_n+1)); # -------------- Step 3 ---------------------------- # The symmetric square of L is Ls2 := c3*tau^3 + c2*tau^2 + c1*tau + c0; # where the c0..c3 below are given by formula (4) in the paper: # sa0, sa1, sa2 := seq(eval(i, _n=_n+1), i=[a0,a1,a2]); c3 := a1*sa2^2*a2; c2 := sa1*a2*(-sa1*a1+sa0*a2); c1 := -sa0*a1*(-sa1*a1+sa0*a2); c0 := -sa1*sa0*a0^2; detL := a0/a2; v1 := -detL; # Note: if you use + detL then you v2 := subs(_n=_n+1, -detL) * v1; # can use R below to solve equations v3 := subs(_n=_n+2, -detL) * v2; # of the form LCLM(tau-r1, tau-r2). # Term product of Ls2 with tau + 1/detL LT := primpart(v3*c3*tau^3 + v2*c2*tau^2 + v1*c1*tau + c0, tau); R := LREtools[ratpolysols](add(coeff(LT, tau, i)*_u(_n+i), i=0..3), _u(_n), {}, output=basis); if R=NULL or R={} or R=[] or R=0 then # Under the assumption that L has no hypergeometric solutions # we can now conclude that there are no Liouvillian solutions. return NULL fi; R := R[1]; # ---------------- Step 4 ------------------- eqn13 := (g*_u(_n) + _u(_n+1))^2 = d0*_u(_n)^2 + d1*_u(_n+1)^2 + d2*_u(_n+2)^2; isz := expand(lhs(eqn13) - subs(sub2, rhs(eqn13))); # Used sub2 from Step 2. isz := coeff(isz, _u(_n), 2), coeff(isz, _u(_n+1),2), coeff(coeff(isz,_u(_n),1),_u(_n+1),1); d_values := solve({isz}, {d0, d1, d2}); # ---------------- Step 5 ------------------- # Compute the substitution given in Step 5 sub5 := _u(_n)^2 = R*S(_n), _u(_n+1)^2 = -eval(R, _n=_n+1) * detL * S(_n), _u(_n+2)^2 = eval(R, _n=_n+2) * eval(detL, _n=_n+1) * detL * S(_n); A := normal(subs(d_values, subs(sub5, rhs(eqn13))) / S(_n)); A := collect(A, g, normal); # ---------------- Step 6 ------------------- g := SolveQuadEquation(A, g, _n); # ---------------- Step 7 ------------------- sg := eval(g, _n=_n+1); a, b := normal(a1/a2), normal(a0/a2); den := normal(g*sg - g*a + b); if den=0 then error "Input is reducible, a factor is", tau+g fi; Lhat := tau^2 + normal(b*eval(den, _n=_n+1)/den); Lhat, _u(_n) = normal((sg-a)/den)*_v(_n) - _v(_n+1)/den end: # Solve a quadratic equation and make sure that the square root gets simplified to a rational function. SolveQuadEquation := proc(A, g, _n) local a,b,c,d,sqd,i; c,b,a := seq(coeff(A,g,i),i=0..2); d := collect(b^2-4*a*c, _n, factor); d := sqrfree(d,_n); if member(1, [seq(i[2] mod 2, i=d[2])]) then error "math error, discriminant is not a square" fi; sqd := sqrt(d[1])*mul(i[1]^(i[2]/2), i=d[2]); normal( (-b+sqd)/(2*a) ) end: # ***************************************************************************************************************************************************************** # redid below procedure to pass in results from find_arbitrary and fixed so it doesn't require sequential terms (not necessary to test linear dependence) # and so we no longer worry about zeros as long as they are not arbitrary and we don't have them as both terms # this, of course, means the part for Liouvillian is rewritten as well (originally was only taking a pair, e.g. odd then the next even) findtwo := proc(h, L, arbitrary) # L = [offset, eqn, sequence], arbitrary are those terms of the sequence for input that are ... arbitrary # h = [rhs of transformation, eqn for dependent of transformation (optional)] # Recall that Database returns trans in terms of procedures (e.g. A000085) local counter, positions, new_terms, beginning, rec, sequence, i, cond, condL, a, ans, eqns, fns; counter, positions, new_terms := 0, [], []; beginning, rec, sequence := op(L); rec := solve(rec, _u(_n+2)); for i from beginning to beginning + 40 do if i > nops(sequence) + beginning - 1 then # we were unable to generate additional terms in the try statement below return NULL fi; if nargs = 3 and i in arbitrary then # if arbitrary is passed into findtwo and i is in there then skip this number next fi; cond := not has(h[1], _v) and i - beginning >= nops(sequence) and not ({i-2, i-1} subset arbitrary); # not Liouvillian condL := has(h[1], _v) and i - beginning >= nops(sequence)-1 and not (i-2 in arbitrary); # Liouvillian part needs an extra term hence the -1 if (cond or condL) and not (i in arbitrary) then # if we've run out of terms of sequence that were input try sequence := [op(sequence), eval(rec, [_u(_n+1)=sequence[-1], _u(_n)=sequence[-2], _n=i-2])] catch "numeric exception: division by zero": return NULL end try fi; try a[i] := eval(h[1], _n=i); if not has(a[i], _v) then # i.e. if in Database (not in the form that FindLiouvillian returns) if type(a[i], realcons) then positions := [op(positions), i]; new_terms := [op(new_terms), a[i]]; counter := counter + 1; # if counter = 2 and a[i-1] = 0 and a[i] = 0 then # counter := 1; # next # else # counter := counter + 1 # fi # else # counter := 1; # next fi else # i.e. if Liouvillian eqns, fns:=[seq(subs(_n=i+j, sequence[i-beginning+1+j] = h[1]), j=0..1)], {seq(_v(i+j), j=0..1)}; eqns := subs(_v(i+2) = solve(subs(_n=i, h[2]), _v(i+2)), eqns); ans := solve(eqns, fns); if ans <> [] and ans <> NULL then #assign(op(ans)); #return i, [_v(i), _v(i+1)], sequence return i, ans, sequence fi fi; catch "numeric exception: division by zero": # counter := 1; next end try; if counter = 2 then # return i-1, a[i-1], a[i], sequence # i.e. returns offset and two terms from new seq and old (i.e. input) seq return positions, new_terms, sequence # i.e. returns two offsets and terms from new seq and also the old (i.e. input) seq # positions means position in the sequence, i.e. not position in the list # e.g. given [3,[1,2,3]] position = 3 references 1 fi; od; return NULL end: # ***************************************************************************************************************************************************************** findMatrix := proc(recs, offsets, list, _n, arbitrary) local temp1, temp2, substitutions, first, g, h, counter, start, i, a, b, c, trya, tryb, tryc, positions; temp1 := eval( solve(recs[1], _u(_n+2)), _n=_n-2 ); temp2 := eval( solve(recs[1], _u(_n+2)), _n=k-2 ); substitutions := {'eqn' = subs('_u'=replace, temp2), 'eqn2' = subs('_u'=replace, temp1), 'lst' = list, 'offsetH' = offsets[1]}; assign(first, subs(substitutions, eval(AssignName))); # AssignName creates a function, first(n), for the input equation (recs[1]) g, h := recs[2], recs[3]; counter := 1; start := max(op(offsets)); a,b,c,positions := [],[],[],[]; for i from start to start + 50 do if i in arbitrary then next fi; try tryc := eval(first(_n), _n=i); trya := eval(g, _n=i); tryb := eval(h, _n=i) catch "numeric exception: division by zero": counter := 1; next end try; c := [op(c), tryc]; a := [op(a), trya]; b := [op(b), tryb]; positions := [op(positions), i]; counter := counter + 1; if counter = 3 then return Matrix([[a[1], b[1]], [a[2], b[2]]]), , positions[1] fi; od end: # ***************************************************************************************************************************************************************** AssignName := proc(k) global _n; # _n will be replaced with the input variable (which needs to be global, i.e. the same as the one in here) local seq, temp1; option remember; # proc used to be called H that's why I have such names as offsetH seq := lst; _n := '_n'; # w/o this had trouble with, e.g. the external call n:=3; A...(n); b/c that n:=3 made all the _n's which become n's equal to 3 if has(k, _n) then if not type(k-_n, integer) then return 'procname'(k) else if k-_n > 1 then return collect(eval(collect(subs(_n=k, subs(replace = 'procname', 'eqn2')), 'procname', normal)), 'procname', normal) else return 'procname'(k) fi fi fi; if not type(k, integer) or k < offsetH then return 'procname'(k) fi; if k < offsetH + nops(seq) then return seq[k + 1 - offsetH] fi; if (not has(k,_n)) and nops(seq) < 2 then return "FAIL" fi; temp1 := eval(collect(subs(_n=k, subs(replace = 'procname', 'eqn2')), 'procname', normal)); simplify(temp1, GAMMA) end: # ***************************************************************************************************************************************************************** NormalizeIrrMonicPoly := proc(f,x) # returns shifted poly such that average of real part of roots is in (-1, 0] local k,d; if lcoeff(f,x)<>1 then error "input not monic" fi; k := degree(f,x); d := floor(coeff(f,x,k-1)/k); expand(subs(x=x-d,f)) end: # ***************************************************************************************************************************************************************** NormalizeRat := proc(f,x) local facts,ff,v,i,k; membertype(`!`, [op(f)], k); facts := 1; ff := f; while type(k,integer) do facts := facts * _n! * _n^( op(op(ff)[k]) - _n ); # I should have trouble with something like (_n!)^2, but (_n-2)! is fine ff := ff/op(ff)[k]; k := 'k'; membertype(`!`, [op(ff)], k); od; v := factors(ff); v[1] * mul(NormalizeIrrMonicPoly(i[1],x)^i[2],i=v[2]) * facts end: # ***************************************************************************************************************************************************************** Find2F1 := proc(Linp, tau, _n) local L, Det, sfDet, d, dhat, k, i, alpha1, alpha2, c, b0, b1, b1hat, Lhat, dt, possac, tmp, rts, a, ysols, seriesorder, yratio1, yratio2, ser1, ser2, c01, c11, c02, c12, j, c0j, c1j, z, bhat, bhatset, b, Lbase, Lu, gttrans,m; # ---------------- Step 1 ------------------- L := collect(primpart(Linp, tau), tau, normal); if 2*degree(coeff(L, tau, 1), _n) > degree(coeff(L, tau, 2), _n) + degree(coeff(L, tau, 0), _n) or not type((degree(coeff(L, tau, 2), _n) - degree(coeff(L, tau, 0), _n))/2, integer) then return "FAIL" fi; # ---------------- Step 2 ------------------- Det := NormalizeRat(coeff(L, tau, 0)/coeff(L, tau, 2), _n); # shift Normalized Det := denom(Det)*numer(Det); sfDet := sqrfree(Det); d, dhat := 1, 1; k := sfDet[1]; for i in sfDet[2] do if i[2] mod 2 = 1 then d,dhat := d*i[1], dhat*i[1]^((i[2]-1)/2) else d,dhat := d*i[1], dhat*i[1]^(i[2]/2) fi od; if degree(d, _n) = 2 then alpha1, alpha2 := solve(d, _n); c := alpha2 - alpha1 elif d = 1 then c := 1 else return "FAIL" fi; # ---------------- Step 3 ------------------- b0, b1 := coeff(L, tau, 0)*subs(_n=_n+1, dhat)*dhat/coeff(L, tau, 2), -coeff(L, tau, 1)*subs(_n=_n+1, dhat)/coeff(L, tau, 2); b1hat := NormalizeRat(b1, _n); Lhat := tau^2 + b1*tau + b0; dt := denom(normal(coeff(L, tau, 1)/coeff(L, tau, 2))); # variable name dt = "denominator of tau" # ---------------- Step 4 ------------------- if c <> 1 then possac := {[-alpha1, c]} # Theorem 4.4 shows that either case 1a or 1b is sufficient else # variable name possac = "possible a and c values" # B/c of the way we're finding gt-trans we work in Q so need only consider linear factors of dt that are in Q possac, rts := {},{}; for i in factors(dt)[2] do # Avoid repeats of a mod Z while preserving a tmp := NormalizeIrrMonicPoly(i[1], _n); possac := possac union `if`(degree(i[1], _n) = 1 and not tmp in rts, {[i[1]-_n-1, 1]}, NULL); rts := rts union {tmp}; od; fi; # ---------------- Step 5 ------------------- for i in possac do a,c := i[1],i[2]; ysols := [ (-b1 + sqrt(b1^2-4*b0))/2, (-b1 - sqrt(b1^2-4*b0))/2 ]; seriesorder := 6; # Maple has Order=6 by default try yratio1 := series(ysols[2]/ysols[1], _n=infinity, seriesorder); yratio2 := series(ysols[1]/ysols[2], _n=infinity, seriesorder); catch: try seriesorder := seriesorder + 5; yratio1 := series(ysols[2]/ysols[1], _n=infinity, seriesorder); yratio2 := series(ysols[1]/ysols[2], _n=infinity, seriesorder); catch: try seriesorder := seriesorder + 5; yratio1 := series(ysols[2]/ysols[1], _n=infinity, seriesorder); yratio2 := series(ysols[1]/ysols[2], _n=infinity, seriesorder); catch: return NULL end try end try end try; ser1 := collect(simplify(convert(yratio1, polynom),symbolic),_n,factor); ser2 := collect(simplify(convert(yratio2, polynom),symbolic),_n,factor); if not degree(ser1, _n) = 0 or not degree(ser1, _n) = 0 then return "FAIL" fi; c01, c11 := coeff(ser1, _n, 0), coeff(ser1, _n, -1); c02, c12 := coeff(ser2, _n, 0), coeff(ser2, _n, -1); for j to 2 do z, c0j, c1j := `if`(j=1, op([1 - c02, c01, c11]), op([1 - c01, c02, c12])); if z=0 then return NULL fi; bhat := (c + c1j/c0j)/2; bhatset := `if`(c11=0, {c/2}, {bhat, bhat + 1/2}); # check for the case c=2b for b in bhatset do Lbase := (z-1)*(a+_n+1)*_u(_n+2) + (2*(a+_n+1)-c-(a+_n+1)*z+b*z)*_u(_n+1) + (c-a-_n-1)*_u(_n); Lu := add(coeff(L, tau, m)*_u(_n+m), m=0..2); gttrans := findgaugeandterm([Lu, _u(_n)], [Lbase, _u(_n), _twoFone]); if gttrans <> NULL then if c = 1 then return gttrans, '_twoFone(_n)' = hypergeom( [a+_n,b], [c], z), [a,b,c,z] else return gttrans, '_twoFone(_n)' = hypergeom( [a+_n,b], [c], z), [a,b,c,z], '_twoFone(_n)' = GAMMA(a+_n+1-c)/GAMMA(a+_n)*hypergeom( [a+_n+1-c,b+1-c], [2-c], z), [a+1-c, b+1-c, 2-c, z] fi fi od od od end: # ***************************************************************************************************************************************************************** ridconstant := proc(result, params) # gets rid of constants (other than _c) for optimize local a, b, c, z, sub, unsub, newresult, h, pp, const, i; a,b,c,z := op(params); sub := [hypergeom([a+_n+1, b], [c], z) = x^2, hypergeom([a+_n, b], [c], z) = x]; unsub := [x^2 = hypergeom([a+_n+1, b], [c], z), x = hypergeom([a+_n, b], [c], z)]; newresult := subs(sub, result); newresult := simplify(newresult, GAMMA); h := content(newresult, x, 'pp'); const := `if`(has({op(newresult)}, -1), -1, 1); for i in op(h) do if (not has(i,_n)) and (not has(i,_c)) then const := const*i fi od; h/const*(subs(op(unsub), collect(pp, x, factor))) end: # ***************************************************************************************************************************************************************** optimize := proc(initresult, initparams) # This tries to find nicer parameters for 2F1 local result, params, Lu2, a, b, c, z, Lshifted, better, improve, i, aneg, u2, id, hg, tauid, sub, newresult, k,j; global _u; result := ridconstant(initresult, initparams); params := initparams; Lu2 := -((a+_n+1)*(2-z)+b*z-c)/((a+_n+1)*(z-1))*_u(_n+1) - (c-a-_n-1)/((a+_n+1)*(z-1))*_u(_n); Lshifted := [ #[abc, -1, hypergeom([a+_n, b], [c], z) = (c-1)/((b-1)*z)*(hypergeom([a+_n, b-1], [c-1], z) - hypergeom([a+_n-1, b-1], [c-1], z))], [a, -1, hypergeom([a+_n, b], [c], z) = (-2*a+a*z-2*_n+_n*z+2-z-b*z+c)/((z-1)*(a+_n-1))*hypergeom([a+_n-1, b], [c], z) + (-c+a+_n-1)/((z-1)*(a+_n-1))*hypergeom([a+_n-2, b], [c], z)], [a, +1, hypergeom([a+_n, b], [c], z) = -(-2*a+a*z-2*_n+_n*z-2+z-b*z+c)/(1-c+a+_n)*hypergeom([a+_n+1, b], [c], z) + (z-1)*(a+_n+1)/(1-c+a+_n)*hypergeom([a+_n+2, b], [c], z)], [b, -1, hypergeom([a+_n, b], [c], z) = (b-1-a-_n)*hypergeom([a+_n, b-1], [c], z)/(b-1) + (a+_n)*hypergeom([b-1, a+_n+1], [c], z)/(b-1)], [b, +1, hypergeom([a+_n, b], [c], z) = (a+_n)*(z-1)/(b+1-c)*hypergeom([a+_n+1, b+1], [c], z) - (c-a-_n-b-1)/(b+1-c)*hypergeom([b+1, a+_n], [c], z)], [c, -1, hypergeom([a+_n, b], [c], z) = -(c-1)*(a+_n+(b-c+1)*z)*hypergeom([a+_n, b], [c-1], z)/((a+_n-c+1)*(b-c+1)*z) + (c-1)*(a+_n)*(1-z)*hypergeom([b, a+_n+1], [c-1], z)/((a+_n-c+1)*(b-c+1)*z)], [c, +1, hypergeom([a+_n, b], [c], z) = (a+_n)/c*hypergeom([a+_n+1, b], [c+1], z) - (a+_n-c)/c*hypergeom([b, a+_n], [c+1], z)], [z, -z + z/(z-1), hypergeom([a+_n, b], [c], z) = (1-z)^(-a-_n)*hypergeom([a+_n, c-b], [c], z/(z-1))], [z, -z + z/(z-1), hypergeom([a+_n, b], [c], z) = (1-z)^(-b)*hypergeom([c-a-_n, b], [c], z/(z-1))], [a, -2, hypergeom([b, a+_n], [c], z) = (6+3*a^2-3*a^2*z+6*a*_n-3*a*c+a^2*z^2-3*a*z^2+3*_n^2-3*_n^2*z-3*_n*c-3*_n*z^2+_n^2*z^2+2*z^2+3*b*z^2-z*c+b^2*z^2+2*a*z^2*_n-2*a*z^2*b-2*_n*z^2*b-6*a*_n*z+4*a*b*z+a*z*c+4*_n*b*z+_n*z*c-2*b*z*c+c^2-9*_n+4*c-6*z-9*a+9*a*z+9*_n*z-6*b*z)*hypergeom([b, a+_n-2], [c], z)/((z-1)^2*(a+_n-1)*(a+_n-2))+(-2*a+a*z-2*_n+_n*z+2-z-b*z+c)*(-c+a-2+_n)*hypergeom([b, a-3+_n], [c], z)/((z-1)^2*(a+_n-1)*(a+_n-2))], [a, +2, hypergeom([b, a+_n], [c], z) = (6+3*a^2-3*a^2*z+6*a*_n-3*a*c+a^2*z^2+3*a*z^2+3*_n^2-3*_n^2*z-3*_n*c+3*_n*z^2+_n^2*z^2+2*z^2-3*b*z^2+2*z*c+b^2*z^2+2*a*z^2*_n-2*a*z^2*b-2*_n*z^2*b-6*a*_n*z+4*a*b*z+a*z*c+4*_n*b*z+_n*z*c-2*b*z*c+c^2+9*_n-5*c-6*z+9*a-9*a*z-9*_n*z+6*b*z)*hypergeom([b, a+_n+2], [c], z)/((2-c+a+_n)*(1-c+a+_n))-(-2*a+a*z-2*_n+_n*z-2+z-b*z+c)*(z-1)*(a+_n+2)*hypergeom([b, a+3+_n], [c], z)/((2-c+a+_n)*(1-c+a+_n))], [b, -2, hypergeom([b, a+_n], [c], z) = -(a+_n)*(-b*z+2*b+z-2+a*z+_n*z-c)*hypergeom([b-2, a+_n+1], [c], z)/((z-1)*(b-2)*(b-1))+(-2-3*b*z+3*a*z+3*_n*z+2*b*a-b^2+2*b*_n+b^2*z+a^2*z+_n^2*z-a*c-_n*c-2*b*a*z-2*b*_n*z+2*a*_n*z-2*_n+2*z-2*a+3*b)*hypergeom([b-2, a+_n], [c], z)/((z-1)*(b-2)*(b-1))], [b, +2, hypergeom([b, a+_n], [c], z) = (a+_n)*(z-1)*(a*z+_n*z+2-z-b*z+2*b-c)*hypergeom([b+2, a+_n+1], [c], z)/((-b-1+c)*(-b-2+c))+(2+a*z+_n*z+2*b*a+b^2+2*b*_n+a^2*z+_n^2*z-a*c-_n*c+2*a*_n*z-a*z*c-_n*z*c+c^2-2*c*b+2*_n-3*c+2*a+3*b)*hypergeom([b+2, a+_n], [c], z)/((-b-1+c)*(-b-2+c))], [c, -2, hypergeom([b, a+_n], [c], z) = (c-1)*(c-2)*(5*a*z+5*_n*z+a^2*z+_n^2*z+a*c+_n*c+2*b*a*z+2*b*_n*z+2*a*_n*z-3*a*z*c-3*_n*z*c-2*b*z^2*c+b^2*z^2+3*b*z^2+2*z^2-3*z^2*c+z^2*c^2-2*_n-2*a)*hypergeom([b, a+_n], [c-2], z)/((a+_n-c+2)*(-b-2+c)*z^2*(a+_n-c+1)*(-b-1+c))+(c-1)*(c-2)*(a+_n)*(z-1)*(_n*z-2*z*c+b*z+a*z+3*z-2+c)*hypergeom([b, a+_n+1], [c-2], z)/((a+_n-c+2)*(-b-2+c)*z^2*(a+_n-c+1)*(-b-1+c))], [c, +2, hypergeom([b, a+_n], [c], z) = -(a+_n)*(-z+b*z-2*z*c+a*z+_n*z+c)*hypergeom([b, a+_n+1], [c+2], z)/(c*(c+1)*(z-1))+(a+_n-c-1)*(a*z+_n*z-z*c+c)*hypergeom([b, a+_n], [c+2], z)/(c*(c+1)*(z-1))] ]; # This section is just verifying that I typed in the identities (Lshifted) correctly # This won't work for the z relations without evalf (and even then it'll be off by something like 1*10^DIGITS) # shouldbezero := []; # for i in Lshifted2 do # shouldbezero := [op(shouldbezero), simplify(eval(lhs(i[3])) - eval(rhs(i[3])), hypergeom)] # od; # return op(shouldbezero); better := 1; while better <> 0 do better := 0; improve := [0,0,0,0]; for i in Lshifted do member(i[1], [a,b,c,z,abc], 'k'); if i[1] = c or i[1] = abc then # Making sure c won't become a non-pos integer (function won't exist) if type(params[3] + i[2], nonposint) then next fi fi; # *************** below part is in progress (and probably unnecessary) ************************************* # if i[1] = abc then # u2 := subs({a=a+i[2], b=b+i[2], c=c+i[2]}, Lu2); # id := i[-1]; # hg := eval(hypergeom([a+n, b], [c], z), {a=a+i[2], b=b+i[2], c=c+i[2]}); # a,b,c,z := op(params); # hg, Lu2, u2, id := eval(hg), eval(Lu2), eval(u2), eval(id); # _u := j -> subs(n=j, hg); # tauid := eval(id, n=n+1); # newresult := eval(result, [id, tauid]); # else if i[1] = z then # Making sure we're in range where identity is valid and don't want to leave radius of convergence # *************************************** Maybe leave off the `or ...' below ************************************************** if (type(params[4], realcons) and params[4] >= 1) or abs(params[4]) >= 1/2 then next fi; hg := eval(hypergeom([a+_n, b], [c], z), i[1]=i[1]+i[2]); a,b,c,z := op(params); _u := j -> subs(_n=j, eval(hg)); id := eval(i[-1]); tauid := eval(id, _n=_n+1); newresult := simplify(eval(result, [id, tauid])); # Need simplify for stuff like (4/3)^_n*(3/4)^(1/2-_n) else aneg := `if`(i[1] = a and i[2] < 0, true, false); u2 := subs(i[1]=i[1]+i[2], Lu2); id := i[-1]; hg := eval(hypergeom([a+_n, b], [c], z), i[1]=i[1]+i[2]); a,b,c,z := op(params); hg, Lu2, u2, id := eval(hg), eval(Lu2), eval(u2), eval(id); _u := j -> subs(_n=j, hg); tauid := eval(id, _n=_n+1); sub := _u(_n+2) = u2; tauid := `if`(aneg, collect(eval(eval(tauid, id)), hypergeom, factor), collect(eval(eval(tauid, sub)), hypergeom, factor)); newresult := eval(result, [id, tauid]); fi; # fi; sub := `if`(k=5, {seq(j=params[j] + 1, j=1..3)}, {k = params[k] + eval(i[2])}); newresult := ridconstant(newresult, subsop(op(sub), params)); # if numboccur(newresult, hypergeom) <= 2 and length(result) - length(newresult) > improve[-1] then # The numboccur is an aesthetic consideration if length(result) - length(newresult) > improve[-1] then improve := [k, eval(i[2]), newresult, length(result) - length(newresult)]; # need eval for cases with z better := 1 fi; a,b,c,z,_u := 'a','b','c','z','_u'; Lu2 := -((a+_n+1)*(2-z)+b*z-c)/((a+_n+1)*(z-1))*_u(_n+1) - (c-a-_n-1)/((a+_n+1)*(z-1))*_u(_n); od; if better <> 0 then k, result := improve[1], improve[3]; params := subsop(k=params[k] + improve[2], params); fi; od; return result end: # ***************************************************************************************************************************************************************** # findgaugeandterm only gets called after p-curvature is checked # Warning: findgaugeandterm doesn't work well now w/o a name for 2nd sequence # findgaugeandterm returns NULL not FAIL findgaugeandterm := proc(input1, input2) # need to have name(s) for input2 - it assigns to -1 entry in list # input [rec eqn, name ( e.g. u(n) ), [offset, [sequence]] (optional)], [rec eqn, name, [offset, [sequence]] (optional), rec name (e.g. 'A000172'), ... (other rec is optional)] if IC's in input1 then must have in input2 # findgaugeandterm will also assign input2 to its name local recurrence, findsolution, i, LL, temp1, rec, offset, arbitrarylist, det, right, something, arbitrary, which, ratio, mpos, mneg, v, fg, temp2, pos, substitutions, fctrs, hg, left, j, indx, terms, mat, sequence, c, a, M, two, MAT, x, y, k; if nops(input2) = 8 then recurrence := [ input1, input2[1..4], input2[5..8] ] else recurrence := [ input1, input2 ]; fi; if nops(recurrence[1]) = 3 then findsolution := true else findsolution := false; fi; offset := [seq(0, i=1..nops(recurrence))]; for i to nops(recurrence) do LL[i] := []; temp1 := subs({op(0, recurrence[i][2]) = _u, op(1, recurrence[i][2]) = _n}, recurrence[i][1]); if hastype(temp1, equation) then temp1 := lhs(temp1)-rhs(temp1) fi; if (i = 1 and findsolution) or (i > 1 and nops(recurrence[i]) = 4) then offset[i], LL[i] := op(recurrence[i][3]); arbitrarylist[i] := find_arbitrary(temp1, [offset[i], LL[i]]) fi; rec[i] := temp1/coeff(temp1, _u(_n+2)); det[i] := simplify(NormalizeRat(coeff(rec[i], _u(_n)) , _n)); od; right, something := [rec[1], 0, 0], false; c := NULL; for which from 2 to 3 do # which references the (up to two) recs in input2 (2 and 3 since 1 references input1) if which = 3 and nops(recurrence) <> 3 then next fi; if LL[1] <> [] then # i.e. there are ICs arbitrary := `union`(arbitrarylist[1], arbitrarylist[which]) fi; ratio := root( det[1]/det[which], 2, symbolic ); # root b/c for f,g of order m,n: detfg=detf^n*detg^m for r the 1st order guy: deta=detb*detr^2 if type(ratio, ratpoly(rational)) then mpos := multrecs([ coeff(rec[which], _u(_n+1)), coeff(rec[which], _u(_n)) ], [-ratio], _u, _n); mneg := multrecs([ coeff(rec[which], _u(_n+1)), coeff(rec[which], _u(_n)) ], [ratio], _u, _n); fg := findgauge( rec[1], subs(_u=v, mpos), _u, v, _n); pos := 2; if fg <> "FAIL" then pos := 1; something := true else fg := findgauge( rec[1], subs(_u=v, mneg), _u, v, _n); if fg <> "FAIL" then pos := 0; something := true fi fi; if pos <> 2 then fg := eval(fg, _c[1]=1); temp1 := eval( solve(rec[which], _u(_n+2)), _n=_n-2 ); temp2 := eval( solve(rec[which], _u(_n+2)), _n=k-2 ); unprotect(recurrence[which][-1]); if nops(recurrence[which]) = 4 then substitutions:={'eqn' = subs('_u'=replace, temp2), 'eqn2' = subs('_u'=replace, temp1), 'lst' = LL[which], 'offsetH' = offset[which]}; assign(recurrence[which][-1], subs(substitutions, eval(AssignName))) else substitutions := {'eqn' = subs('_u'=replace, temp2), 'eqn2' = subs('_u'=replace, temp1), 'lst' = [], 'offsetH' = -infinity}; assign(recurrence[which][-1], subs(substitutions, eval(AssignName))) fi; protect(recurrence[which][-1]); # hg := LREtools[hypergeomsols](a(_n+1) + (-1)^pos*ratio*a(_n), a(_n), {a(1)=1}); # Above is no good: if ratio includes (_n-1) then will return NULL since a(1) leads to something undefined # the below lines are pretty much the same as, but quicker than, unTP := LREtools[hypergeomsols](_u(n+1) + det*_u(n), _u(n), {}) fctrs := factors((-1)^pos*ratio); hg := (-fctrs[1])^_n; fctrs := fctrs[2]; for i in fctrs do if degree(i[1]) = 1 then hg := subs(n=_n, hg*GAMMA(i[1])^i[2]) # subs is silly but some_n's somehow changed to n's else hg := hg*mul(GAMMA(j[1])^j[2], j in convert(evala(AFactors(i[1])), radical)[2])^i[2] fi od; temp1 := subs(v(_n+1) = (-1)^(pos+1)*ratio*recurrence[which][-1](_n+1), v(_n) = recurrence[which][-1](_n), fg); left, right[which] := lhs(temp1), hg*rhs(temp1); # At this point we have our answer (w/o ICs): left = a linear comb of right[1] and right[2] # of course it still gets simplified and transformed into nicer result elsewhere if LL[1] <> [] then # i.e. are there ICs two := findtwo([right[which]], [offset[1], rec[1], LL[1]], arbitrary); if two <> NULL then indx[which], terms[which], sequence := two; # sequence is original (i.e. input) sequence if nops(sequence) > nops(LL[1]) then LL[1] := sequence fi; # findtwo already gives it as position in seq (i.e. not postion in list) so below line changes it to position in list indx[which] := indx[which] - [offset[1],offset[1]] + [1,1]; mat := Matrix([terms[which], [LL[1][indx[which][1]], LL[1][indx[which][2]]]]); if Determinant(mat) = 0 then if terms[which][1] <> 0 then c := LL[1][indx[which][1]]/terms[which][1] elif terms[which][2] <> 0 then c := LL[1][indx[which][2]]/terms[which][2] fi fi fi; if c <> NULL then return c*right[which], [n >= indx[which][1]] fi fi; if something and not findsolution then temp1 := collect(right[which], recurrence[which][-1], factor); return _c*temp1/icontent(temp1) # no point in having extra constants if I'm including an arbitrary constant fi fi fi od; if something and (right[2] <> 0 and right[3] <> 0) then arbitrary := `union`(seq(arbitrarylist[i], i=1..3)); MAT := findMatrix(right, offset, LL[1], _n, arbitrary); c := NULL; if MAT = NULL then userinfo(5, findrel, `Try more ICs or try no ICs`); return NULL elif Determinant(MAT[1]) <> 0 then c := solve({MAT[1][1,1]*x + MAT[1][1,2]*y = MAT[2][1], MAT[1][2,1]*x + MAT[1][2,2]*y = MAT[2][2]}, [x,y]) fi; if c <> NULL then return rhs(c[1][1])*right[2] + rhs(c[1][2])*right[3], [n >= MAT[3]] fi; # c := LinearSolve(MAT[1], MAT[2]) # fi; # # if type(c, Vector) then # return c[1]*right[2] + c[2]*right[3] # fi; if nops(input2) <> 8 then userinfo(5, findrel, `No solution found, but for the same recurrence equation with _u(k), _u(k+1) =`, terma[2], termb[2], `where k =`, indx[2], `Your sequence is`, left = right[2]); return NULL fi fi end: