from gc import collect R = RealIntervalField() Pi = R(pi) L2 = R(log(2)) eR = R(e) S. = PolynomialRing(QQ,24) T. = PolynomialRing(R,16) #------------------------------------------------------------ # power series maxcutoff = 185 """maxcutoff needs to be bigger than 2*cutoff+1 for processcorner (because of quotients) here, just needs to be large enough to apply lemma 4.5 to get all coefficients up to a precision of about 15 digits """ LpxS = sum([-x^j/j for j in range(1,maxcutoff+1)]) LpyS = sum([-y^j/j for j in range(1,maxcutoff+1)]) def d(n): """compute the function d used in the elliptic functions """ return 2*l2-sum([2/((2*m-1)*(2*m)) for m in range(1,n+1)]) A1xS = sum([d(n)/2^(4*n)*binomial(2*n,n)^2*x^n for n in range(0,maxcutoff+1)]) A2xS = -sum([1/2^(4*n+1)*binomial(2*n,n)^2*x^n for n in range(0,maxcutoff+1)]) A3xS = 1+sum([1/2^(4*n+1)*binomial(2*n,n)*binomial(2*n+1,n)*(d(n)-1/((2*n+1)*(2*n+2)))*x^(n+1) for n in range(0,maxcutoff)]) A4xS = -sum([1/2^(4*n+2)*binomial(2*n,n)*binomial(2*n+1,n)*x^(n+1) for n in range(0,maxcutoff)]) A1yS = sum([d(n)/2^(4*n)*binomial(2*n,n)^2*y^n for n in range(0,maxcutoff+1)]) A2yS = -sum([1/2^(4*n+1)*binomial(2*n,n)^2*y^n for n in range(0,maxcutoff+1)]) A3yS = 1+sum([1/2^(4*n+1)*binomial(2*n,n)*binomial(2*n+1,n)*(d(n)-1/((2*n+1)*(2*n+2)))*y^(n+1) for n in range(0,maxcutoff)]) A4yS = -sum([1/2^(4*n+2)*binomial(2*n,n)*binomial(2*n+1,n)*y^(n+1) for n in range(0,maxcutoff)]) #-------------------------- def split(tb): """split off the Log(y)^2 bounds """ if len(tb) == 0: return tb, [[]] elif len(tb[0]) < 3: return tb, [ [] for i in range(len(tb)) ] elif len(tb[0]) == 3: return [ [ tb[i][j] for j in range(2) ] for i in range(len(tb)) ], [ [ tb[i][2] ] for i in range(len(tb)) ] else: raise Exception("log(y)^3 exists!") def num(pol): return T(pol.subs(x=xR,y=yR,Lx=LxR,Ly=LyR,piv=Pi,l2=L2)) def trunc(f,v,e): """truncate single variable f up to degree e in v """ g = 0 for j in range(e+1): g = g + f.coefficient({v:j})*v^j return g def clean(f,v1,v2,e1,e2): """truncate f up to terms of degree e1 and e2 in v1 and v2 resp """ g = 0 for j in range(e1+1): for k in range(e2+1): g = g + f.coefficient({v1:j,v2:k})*v1^j*v2^k return g def condense(tb, errb): return sum(tb[j][k]*errb*LxR^j*LyR^k for k in range(len(tb[0])) for j in range(len(tb))) def varquoterrorbound(x0,y0,mxexp,myexp,cutoff,lowbd): """ evaluate an error bound for power series remainder of ( f(x,y) - f(y,y) )/(x-y) starting from order cutoff+1 assuming the power series f is divisible by x^mxexp and y^myexp and the coeffs all bounded by 1 to be precise, return the bound for power series remainder divided by x^mxexp * y^myexp return a polynomial output (peel off some powers of x and y from cutoff) """ x1 = R(R(x0).upper()) y1 = R(R(y0).upper()) sm = yR^lowbd*(x1^mxexp*(1+(1-x1)*mxexp)-x1^(cutoff+1)*(1+(1-x1)*(1+cutoff)))/(1-x1)^2*y1^(cutoff+1-lowbd)/(1-y1) \ + xR^lowbd*(y1^myexp*(1+(1-y1)*myexp)-y1^(cutoff+1)*(1+(1-y1)*(1+cutoff)))/(1-y1)^2*x1^(cutoff+1-lowbd)/(1-x1) \ + xR^lowbd*y1^lowbd*(x1*y1)^(cutoff+1-lowbd)*(1+(1-x1*y1)*(cutoff+1))/((1-x1*y1)*(1-x1)*(1-y1)) sm = sm/(x1^mxexp*y1^myexp) return sm*R((-1,1)) def reducepoly(p,x0,y0): """reduce the poly to total degree at most 1 in x and y """ x1 = R((0,x0)) y1 = R((0,y0)) q = 0 lydeg = p.degree(LyR) lxdeg = p.degree(LxR) for i in range(lxdeg+1): for j in range(lydeg+1): p0 = p.coefficient({LxR:i, LyR:j}) q0 = 0 for m in p0.monomials(): c = p0.monomial_coefficient(m) dx = m.degree(xR) dy = m.degree(yR) if dx >= dy and dx+dy > 1: m0 = xR*x1^(dx-1)*y1^dy q0 = q0 + c*m0 elif dy > dx and dx+dy > 1: m0 = x1^dx*yR*y1^(dy-1) q0 = q0 + c*m0 else: q0 = q0 + c*m q = q + q0*LxR^i*LyR^j return q def checktinycorner(p,x0,y0,corner="SW",col="black"): """check with x0 and y0 small. convert main term to linear with reducepoly, and error term to linear with varquoterrorbound then check signs of terms multiplying log(x)^i and log(y)^j are correct. finally plug in the lower bound for log(x) and log(y) now the expression is linear, so minimum must be at one of the four vertices. plug in to each and check positive. """ q = reducepoly(p,x0,y0) lxdeg = p.degree(LxR) lydeg = p.degree(LyR) flg = true for i in range(lxdeg+1): for j in range(lydeg+1): if i+j > 0: c = q.coefficient({LxR:i, LyR:j})*(-1)^(i+j) if c.coefficient({xR:1}) < 0 or c.coefficient({yR:1}) < 0 or c.coefficient({xR:0,yR:0}) < 0: print("coefficient of ",i,j) print("c = ",c) flg = false q1 = q.subs(LxR=log(R(x0)), LyR=log(R(y0))) flg = flg and q1.subs(xR=0,yR=0) >= 0 and q1.subs(xR=0,yR=y0) >= 0 and q1.subs(xR=x0,yR=0) >= 0 and q1.subs(xR=x0,yR=y0) >= 0 if corner == "SW": CPlist.append([[0,x0],[0,y0],col]) elif corner == "SE": CPlist.append([[1-x0,1],[0,y0],col]) elif corner == "NW": CPlist.append([[0,x0],[1-y0,1],col]) elif corner == "NE": CPlist.append([[1-x0,1],[1-y0,1],col]) return flg def varerrorbound(x1,y1,mxexp,myexp,cutoff,lowbd): """evaluate an errorbound for power series remainder of f(x,y) starting from order cutoff+1 assuming the power series f is divisible by x^mxexp and y^myexp and the coeffs all bounded by 1 to be precise, return the bound for power series remainder divided by x^mxexp * y^myexp return a polynomial output (peel off some powers of x and y from cutoff) """ return ( yR^lowbd*y1^(cutoff+1-lowbd-myexp) + xR^lowbd*x1^(cutoff+1-lowbd-mxexp) - yR^lowbd*x1^(cutoff+1-mxexp)*y1^(cutoff+1-lowbd-myexp) )/((1-x1)*(1-y1))*R((-1,1)) def errorbound1varpol(v0,var,mxexp,cutoff,deg): v1 = R(R(v0).upper()) return v1^(cutoff+1-mxexp-deg)/(1-v1)*var^deg*R((-1,1)) def logquot(x,y): if x != y: return (log(x)-log(y))/(x-y) else: return 1/x def logquotinterval(x,y): """assume x and y as intervals don't overlap. this uses that logquot is monotonic decreasing in x and in y """ if x.lower() == 0 or y.lower() == 0: return R((logquot(R(x.upper()),R(y.upper())),+infinity)) else: return R((logquot(R(x.upper()),R(y.upper())),logquot(R(x.lower()),R(y.lower())))) def subsH(pol,var,val): """horner evaluation """ s = 0 for j in range(pol.degree(var),-1,-1): s = s*val + pol.coefficient({var:j}) return s def subsHtrunc(pol,var,val,v,mdeg): s = 0 for j in range(pol.degree(var),-1,-1): g = s*val + pol.coefficient({var:j}) s = 0 for k in range(mdeg+1): s = s + g.coefficient({v:k})*v^k return s def subsH2var(pol,v1,v2,x1,x2): """horner evaluation of a poly in v1 and v2 at values x1 and x2 """ s = 0 for j in range(pol.degree(v1),-1,-1): s1 = 0 pol1 = pol.coefficient({v1:j}) for k in range(pol1.degree(v2),-1,-1): s1 = s1*x2 + pol1.coefficient({v2:k}) s = s*x1 + s1 return s def subsHxy(pol,x0,y0): """horner evaluation of a poly in xR and yR """ return subsH2var(pol,xR,yR,x0,y0) def evalRdir(pol,x0,y0,logterm=0,dir=True): """evaluate poly in x,y,Lx,Ly in the real interval field R better to use this when x is smaller than y or can go to 0 since it is an expansion in powers of Lx (which can approach -infty) assuming the power series is divisible by x^mxexp and y^myexp logterm is the coeff of (Lx-Ly) and has to be divided by x-y. it's a function of x we should only use it if x and y intervals don't overlap this assumes degree in Lx and in Ly is at least 1 when logterm is nonzero (so that n=1,m=0 and n=0,m=1 terms are both counted) dir is True for x and False for y """ if logterm != 0: assert pol.degree(LxR) >= 1 assert pol.degree(LyR) >= 1 s = 0 d = 0 if dir else 1 xR0 = R(x0) yR0 = R(y0) lxR0 = log(xR0) lyR0 = log(yR0) if dir: Lv1 = LxR; Lv2 = LyR; lv1R0 = lxR0; lv2R0 = lyR0 else: Lv1 = LyR; Lv2 = LxR; lv1R0 = lyR0; lv2R0 = lxR0 for n in range(pol.degree(Lv1),-1,-1): sv2 = 0 cflv1 = pol.coefficient({Lv1:n}) for m in range(pol.degree(Lv2),-1,-1): cflv1lv2 = cflv1.coefficient({Lv2:m}) s0 = subsHxy(cflv1lv2,xR0,yR0) if n == 1-d and m == d: if logterm != 0: s0 = s0 + subsHxy(logterm,xR0,yR0).subs(LxR=lxR0)/(xR0-yR0) if n == d and m == 1-d: if logterm != 0: s0 = s0 - subsHxy(logterm,xR0,yR0).subs(LxR=lxR0)/(xR0-yR0) sv2 = sv2*lv2R0 + s0 s = s*lv1R0 + sv2 return R(s) def serbd(j): """this is the bound of lemma 4.1 for products of j power series """ if j < 1: return 1 else: return 2^j*j^j/exp(j-1) def termbound(pol): """provide a bound on the absolute value of every coefficient of a polynomial in A1x,A2x,A3x,A4x,Lpx,Lx,A1y,A2y,A3y,A4y,Lpy,Ly,x,y this is vector of vectors, whose [i][j] entry provides the bound for the coeff of Lx^i * Ly^j """ dlx = pol.degree(Lx) dly = pol.degree(Ly) tb = [ [ 0 for k in range(dly+1) ] for j in range(dlx+1) ] for j in range(dlx+1): for k in range(dly+1): termbd = 0 Z = pol.coefficient({Lx:j,Ly:k}) for t in Z.monomials(): totx = t.degree(A1x)+t.degree(A2x)+t.degree(A3x)+t.degree(A4x)+t.degree(Lpx) toty = t.degree(A1y)+t.degree(A2y)+t.degree(A3y)+t.degree(A4y)+t.degree(Lpy) coef = Z.monomial_coefficient(t) Zt = coef*t termbd = termbd + serbd(totx-1)*serbd(toty-1)*abs(Zt.subs(piv=Pi,l2=L2,\ A1x=(3/2),A2x=(1/2),A3x=1,A4x=(1/4),Lpx=1,\ A1y=(3/2),A2y=(1/2),A3y=1,A4y=(1/4),Lpy=1,x=1,y=1)) tb[j][k] = R(termbd) return tb def quot(pol,cutoff): """take quotient of a polynomial f(x,y) which vanishes identically along diagonal by x-y """ quotlist = [ sum([x^(k-1-m)*y^m for m in range(k)] ) for k in range(2*cutoff+2) ] ydeg = pol.degree(y) qpol = 0 for k in range(ydeg+1): qpol = qpol - pol.coefficient({y:k})*quotlist[k] # assert clean( qpol*(x-y) - pol, x, y, cutoff, cutoff ) == 0 return qpol def CP(f,a,b,c,d,corner="SW",col="black"): """this follows the subdivide-and-conquer strategy do we need the last case? or just fold it into penultimate """ var('g') global CPlist print("Checking positivity for ["+str(a)+","+str(b)+"] x ["+str(c)+","+str(d)+"].") if verb: print("Showing fraction checked and number of remaining cases:") vlist = [[a,b,c,d]] plist = [] s = 0; while len(vlist)!= 0: if verb: print(R(s/((b-a)*(d-c))).lower(),len(vlist)) vl0 = vlist[0] a0 = vl0[0] b0 = vl0[1] c0 = vl0[2] d0 = vl0[3] del vlist[0] g = f(R(a0,b0),R(c0,d0)) if (R(g.upper()) <= 0): print(a0,b0,c0,d0) return False elif (R(g.lower()) <= 0): if (d0-c0) > 2*(b0-a0): my = (c0+d0)/2 vlist.append([a0,b0,c0,my]) vlist.append([a0,b0,my,d0]) elif (b0-a0) > 2*(d0-c0): mx = (a0+b0)/2 vlist.append([a0,mx,c0,d0]) vlist.append([mx,b0,c0,d0]) else: mx = (a0+b0)/2 my = (c0+d0)/2 vlist.append([a0,mx,c0,my]) vlist.append([a0,mx,my,d0]) vlist.append([mx,b0,c0,my]) vlist.append([mx,b0,my,d0]) else: s = s + (b0-a0)*(d0-c0) plist.append([(a0+b0)/2,(c0+d0)/2]) if verb: print("Done.") timesofar() if corner == "SW": CPlist.append([[a,b],[c,d],col]) elif corner == "SE": CPlist.append([[1-b,1-a],[c,d],col]) elif corner == "NW": CPlist.append([[a,b],[1-d,1-c],col]) elif corner == "NE": CPlist.append([[1-b,1-a],[1-d,1-c],col]) return True def timesofar(): print("Time so far:",ceil(time.time()-start_time),"seconds.") def totaltime(): print("Total time:",ceil(time.time()-start_time),"seconds.") Lxderiv = 1/x Lpxderiv = 1/(x-1) Lyderiv = 1/y Lpyderiv = 1/(y-1) A1xderiv = A1x/(2*(1-x)) - A2x/x - A3x/(2*x*(1-x)) A2xderiv = A2x/(2*(1-x)) - A4x/(2*x*(1-x)) A3xderiv = A1x/(2*(1-x)) - A3x/(2*(1-x)) - A4x/x A4xderiv = A2x/(2*(1-x)) - A4x/(2*(1-x)) A1yderiv = A1y/(2*(1-y)) - A2y/y - A3y/(2*y*(1-y)) A2yderiv = A2y/(2*(1-y)) - A4y/(2*y*(1-y)) A3yderiv = A1y/(2*(1-y)) - A3y/(2*(1-y)) - A4y/y A4yderiv = A2y/(2*(1-y)) - A4y/(2*(1-y)) def ellmonomialderiv(monom,var): """differentiate a monomial in Lx,Ly,Lpx,Lpy,A1x,A2x,A3x,A4x,A1y,A2y,A3y,A4y (with a coefficient in QQ[piv,l2](x,y) ). var is one of x and y """ mon = S(monom) Lxd = mon.degree(Lx) Lyd = mon.degree(Ly) Lpxd = mon.degree(Lpx) Lpyd = mon.degree(Lpy) A1xd = mon.degree(A1x) A2xd = mon.degree(A2x) A3xd = mon.degree(A3x) A4xd = mon.degree(A4x) A1yd = mon.degree(A1y) A2yd = mon.degree(A2y) A3yd = mon.degree(A3y) A4yd = mon.degree(A4y) if mon == 0: return 0 der = 0 fact = x*y*(1-x)*(1-y) if var == x: if Lxd > 0: der = der + fact*Lxd*Lxderiv/Lx if Lpxd > 0: der = der + fact*Lpxd*Lpxderiv/Lpx if A1xd > 0: der = der + fact*A1xd*A1xderiv/A1x if A2xd > 0: der = der + fact*A2xd*A2xderiv/A2x if A3xd > 0: der = der + fact*A3xd*A3xderiv/A3x if A4xd > 0: der = der + fact*A4xd*A4xderiv/A4x if var == y: if Lyd > 0: der = der + fact*Lyd*Lyderiv/Ly if Lpyd > 0: der = der + fact*Lpyd*Lpyderiv/Lpy if A1yd > 0: der = der + fact*A1yd*A1yderiv/A1y if A2yd > 0: der = der + fact*A2yd*A2yderiv/A2y if A3yd > 0: der = der + fact*A3yd*A3yderiv/A3y if A4yd > 0: der = der + fact*A4yd*A4yderiv/A4y der = der*mon/fact der = der + derivative(mon,var) return der def ellderivpoly(pol,var): """compute derivative of a poly in Lx,Ly,Lpx,Lpy,A1x,A2x,A3x,A4x,A1y,A2y,A3y,A4y with coeffs in rational function field QQ[piv,l2](x,y) """ if pol.parent() == ZZ or pol.parent() == QQ: return 0 der = 0 fact = x*y*(1-x)*(1-y) for m in pol.monomials(): der = der + pol.monomial_coefficient(m)*ellmonomialderiv(m,var)*fact der = der/fact derd = denominator(der) dern = numerator(der) if parent(derd) == ZZ or parent(derd) == QQ: derc = derd else: derc = derd.content() return (dern/derc)/(derd/derc) def ellderiv(rat,var): """compute derivative of a rational function """ ratn = numerator(rat) ratd = denominator(rat) return (ratd*ellderivpoly(ratn,var)-ratn*ellderivpoly(ratd,var))/ratd^2 def absbound(pol,cutoff,eps): """given a poly in LxR, LyR, xR, yR and a bound eps for xR, yR return a matrix whose [i][j] entry is a bound on absolute value of p(i,j)-p_trunc(i,j) where p(i,j) is the poly multiplying LxR^i * LyR^j and p_trunc(i,j) is p(i,j) truncated up to degree cutoff in xR and yR. in particular if cutoff = -1 then we get a bound on the absolute value of p(i,j) """ dlx = pol.degree(LxR) dly = pol.degree(LyR) eb = [ [ 0 for k in range(dly+1) ] for j in range(dlx+1) ] for j in range(dlx+1): for k in range(dly+1): polxy = pol.coefficient({LxR:j,LyR:k}) sm = 0 for i1 in range(polxy.degree(xR)+1): for i2 in range(polxy.degree(yR)+1): if ((i1 > cutoff) or (i2 > cutoff)): sm = sm + abs(R(polxy.coefficient({xR:i1,yR:i2})))*eps^(i1+i2) eb[j][k] = sm return eb def varerrorbounduv(pol,u0,v0,cutx,cuty,lowbd): """evaluate an error bound for poly remainder starting from order cutoff+1 return a polynomial """ dlx = pol.degree(LxR) dly = pol.degree(LyR) eb = 0 for j in range(dlx+1): for k in range(dly+1): polxy = pol.coefficient({LxR:j,LyR:k}) sm = 0 for i1 in range(polxy.degree(uR)+1): for i2 in range(polxy.degree(vR)+1): if (i1 > cutx): sm = sm + abs(R(polxy.coefficient({uR:j,vR:k})))*u0^(i1-lowbd)*uR^lowbd*v0^i2*R((0,1)) elif (i2 > cuty): sm = sm + abs(R(polxy.coefficient({uR:j,vR:k})))*u0^i1*v0^(i2-lowbd)*vR^lowbd*R((0,1)) eb = eb + sm*R((-1,1))*LxR^j*LyR^k return eb def subsHuv(pol,u0,v0): """horner evaluation of a poly in uR and vR """ return subsH2var(pol,uR,vR,u0,v0) # let Kh = K(1/2), Eh = E(1/2). # Then A1(1/2), etc. can be determined # let gq = gamma(1/4) gq = R(gamma(1/4)) g3q = R(sqrt(2))*Pi/gq Kh = gq^2/(4*sqrt(Pi)) Eh = (4*g3q^2 + gq^2)/(8*sqrt(Pi)) A2h = -Kh/Pi A4h = (Eh-Kh)/Pi A1h = A2h*(L2-Pi) A3h = Eh + A4h*L2 # now we're using the change of variable x = (1+u-v)/2, y = (1-u-v)/2. def partial(pol,var): """ here var is one of x,y,u,v """ if var == x: return ellderiv(pol,x) if var == y: return ellderiv(pol,y) if var == u: return (ellderiv(pol,x) - ellderiv(pol,y))/2 if var == v: return (-ellderiv(pol,x) - ellderiv(pol,y))/2 print("bad variable!") return 0 def evalatcenter(rat): ratn = numerator(rat) ratd = denominator(rat) ratn0 = ratn.subs(Lx=-L2,Ly=-L2,Lpx=-L2,Lpy=-L2,A1x=A1h,A2x=A2h,A3x=A3h,A4x=A4h,A1y=A1h,A2y=A2h,A3y=A3h,A4y=A4h,x=R(1/2),y=R(1/2),piv=Pi,l2=L2) ratd0 = ratd.subs(Lx=-L2,Ly=-L2,Lpx=-L2,Lpy=-L2,A1x=A1h,A2x=A2h,A3x=A3h,A4x=A4h,A1y=A1h,A2y=A2h,A3y=A3h,A4y=A4h,x=R(1/2),y=R(1/2),piv=Pi,l2=L2) return ratn0/ratd0 def plugin(pol,cutoffx,cutoffy): # need to truncate all the power series from maxcutoff to cutoff first LpxSc = trunc(LpxS,x,cutoffx) A1xSc = trunc(A1xS,x,cutoffx) A2xSc = trunc(A2xS,x,cutoffx) A3xSc = trunc(A3xS,x,cutoffx) A4xSc = trunc(A4xS,x,cutoffx) LpySc = trunc(LpyS,y,cutoffy) A1ySc = trunc(A1yS,y,cutoffy) A2ySc = trunc(A2yS,y,cutoffy) A3ySc = trunc(A3yS,y,cutoffy) A4ySc = trunc(A4yS,y,cutoffy) polS = subsHtrunc(pol,Lpx,LpxSc,x,cutoffx) polS = subsHtrunc(polS,A1x,A1xSc,x,cutoffx) polS = subsHtrunc(polS,A2x,A2xSc,x,cutoffx) polS = subsHtrunc(polS,A3x,A3xSc,x,cutoffx) polS = subsHtrunc(polS,A4x,A4xSc,x,cutoffx) polS = subsHtrunc(polS,Lpy,LpySc,y,cutoffy) polS = subsHtrunc(polS,A1y,A1ySc,y,cutoffy) polS = subsHtrunc(polS,A2y,A2ySc,y,cutoffy) polS = subsHtrunc(polS,A3y,A3ySc,y,cutoffy) polS = subsHtrunc(polS,A4y,A4ySc,y,cutoffy) return polS def evalRdiruv(pol,x0,y0,ctrx,ctry,logterm=0,dir=True): """dir is True for x and False for y """ s = 0 d = 0 if dir else 1 xR0 = R(x0) yR0 = R(y0) uR0 = xR0-ctrx vR0 = yR0-ctry lxR0 = log(xR0) lyR0 = log(yR0) lydeg = pol.degree(LyR) lxdeg = pol.degree(LxR) if dir: Lv1 = LxR; Lv2 = LyR; lv1R0 = lxR0; lv2R0 = lyR0 else: Lv1 = LyR; Lv2 = LxR; lv1R0 = lyR0; lv2R0 = lxR0 for n in range(pol.degree(Lv1),-1,-1): sv2 = 0 cflv1 = pol.coefficient({Lv1:n}) for m in range(pol.degree(Lv2),-1,-1): cflv1lv2 = cflv1.coefficient({Lv2:m}) s0 = subsHuv(cflv1lv2,uR0,vR0) if n == 1-d and m == d: if logterm != 0: s0 = s0 + subsHxy(logterm,xR0,yR0).subs(LxR=lxR0)/(xR0-yR0) if n == d and m == 1-d: if logterm != 0: s0 = s0 - subsHxy(logterm,xR0,yR0).subs(LxR=lxR0)/(xR0-yR0) sv2 = sv2*lv2R0 + s0 s = s*lv1R0 + sv2 return R(s)