"Load the symbolic integration table for Carlson's elliptic integrals."

"The table was generated using equation 3.5 in the paper"

"Toward Symbolic Integration of Elliptic Integrals."

"This equation is DERIVE function R35(q,j,h,a,b,e) in c1998.mth."

LOAD("c1999.mth")

true

PrecisionDigits:=64

64

"This is an example of calculating the length of a Bezier curve."

"This is a Bezier curve."

v1:=[x=a3*t^3+a2*t^2+a1*t+x0,y=b3*t^3+b2*t^2+b1*t+y0]

[x=a1*t+a2*t^2+a3*t^3+x0,y=b1*t+b2*t^2+b3*t^3+y0]

"The curve starts at location (x0,y0) when t=0."

"The curve finishes at location (x3,y3) when t=1."

v2:=LIM(LIM(LIM(v1,x,x3),y,y3),t,1)

[x3=a1+a2+a3+x0,y3=b1+b2+b3+y0]

"The slope and the rate of change of slope at t=0 is determined by the contro~
l point (x1,y1)."

v3:=VECTOR(RHS(k),k,DIF(v1,t))

[a1+2*a2*t+3*a3*t^2,b1+2*b2*t+3*b3*t^2]

v4:=LIM([3*(x1-x0),3*(y1-y0)]=v3,t,0)

[3*(x1-x0)=a1,3*(y1-y0)=b1]

"The slope and the rate of change of slope at t=1 is determined by the contro~
l point (x2,y2)."

v5:=LIM([3*(x3-x2),3*(y3-y2)]=v3,t,1)

[3*(x3-x2)=a1+2*a2+3*a3,3*(y3-y2)=b1+2*b2+3*b3]

"Solve for a1, a2, a3, b1, b2, and b3 in terms of x0, x1, x2, x3, y0, y1, y2,~
 and y3."

v6:=APPEND(v2,v4,v5)

[x3=a1+a2+a3+x0,y3=b1+b2+b3+y0,3*(x1-x0)=a1,3*(y1-y0)=b1,3*(x3-x2)=a1+2*a2+3*~
a3,3*(y3-y2)=b1+2*b2+3*b3]

v7:=(SOLUTIONS(v6,[a1,a2,a3,b1,b2,b3])) SUB 1

[3*(x1-x0),3*(x0-2*x1+x2),-x0+3*x1-3*x2+x3,3*(y1-y0),3*(y0-2*y1+y2),-y0+3*y1-~
3*y2+y3]

"INT(int1,t,0,1) is the length of a Bezier curve. It is an elliptic integral.~
"

'DIF(s,t,1)='SQRT(DIF(x,t,1)^2+DIF(y,t,1)^2)

DIF(s,t,1)=SQRT(DIF(x,t,1)^2+DIF(y,t,1)^2)

length='INT(SQRT(DIF(x,t,1)^2+DIF(y,t,1)^2),t,0,1)

length=INT(SQRT(DIF(x,t,1)^2+DIF(y,t,1)^2),t,0,1)

int1:=SQRT(v3 . v3)

SQRT(a1^2+2*a1*t*(2*a2+3*a3*t)+4*a2^2*t^2+12*a2*a3*t^3+9*a3^2*t^4+(b1+2*b2*t+~
3*b3*t^2)^2)

"Pick a actual Bezier curve to check the calculations."

"In this example Bezier curve [[x0,y0],[x1,y1],[x2,y2],[x3,y3]]=v8."

v8:=[[6,8],[1,10],[7,3],[4,4]]

[[6,8],[1,10],[7,3],[4,4]]

"Solve for [a1,a2,a3,b1,b2,b3] in this example."

v10:=LIM(LIM(LIM(LIM(LIM(LIM(LIM(LIM(v7,x0,v8 SUB 1 SUB 1),y0,v8 SUB 1 SUB 2)~
,x1,v8 SUB 2 SUB 1),y1,v8 SUB 2 SUB 2),x2,v8 SUB 3 SUB 1),y2,v8 SUB 3 SUB 2),~
x3,v8 SUB 4 SUB 1),y3,v8 SUB 4 SUB 2)

[-15,33,-20,6,-27,17]

int2:=INT(LIM(LIM(LIM(LIM(LIM(LIM(int1,a1,v10 SUB 1),a2,v10 SUB 2),a3,v10 SUB~
 3),b1,v10 SUB 4),b2,v10 SUB 5),b3,v10 SUB 6),t,0,1)

3*INT(SQRT(689*t^4-1492*t^3+1076*t^2-292*t+29),t,0,1)

"refint is a check value calulated to 64 digits."

refint:=1809305842082148221456654802739005603491766996504376040840721629/2500~
00000000000000000000000000000000000000000000000000000000000

1809305842082148221456654802739005603491766996504376040840721629/250000000000~
000000000000000000000000000000000000000000000000000

"Do numeric integration to show this check value is correct."

APPROX(int2-refint,12)

-1.42458839316*10^(-11)

"Factor the polynomial to put the integral in a standard form."

v11:=LIM(LIM(LIM(LIM(LIM(LIM(v3 . v3,a1,v10 SUB 1),a2,v10 SUB 2),a3,v10 SUB 3~
),b1,v10 SUB 4),b2,v10 SUB 5),b3,v10 SUB 6)

9*(689*t^4-1492*t^3+1076*t^2-292*t+29)

v12:=FACTORS(FACTOR(v11,Complex))

[[t-SQRT(SQRT(6005)/1378+23377/474721)-373/689+#i*(7/689-SQRT(SQRT(6005)/1378~
-23377/474721)),1],[t-SQRT(SQRT(6005)/1378+23377/474721)-373/689+#i*(SQRT(SQR~
T(6005)/1378-23377/474721)-7/689),1],[t+SQRT(SQRT(6005)/1378+23377/474721)-37~
3/689-#i*(SQRT(SQRT(6005)/1378-23377/474721)+7/689),1],[t+SQRT(SQRT(6005)/137~
8+23377/474721)-373/689+#i*(SQRT(SQRT(6005)/1378-23377/474721)+7/689),1],[620~
1,1]]

"km is the scale factor."

km:=SQRT(v12 SUB 5 SUB 1)

3*SQRT(689)

IA(m):=

IA(m)

"Now do partial fraction expansion using the equation 2.19 in the paper"

"Toward Symbolic Integration of Elliptic Integrals."

"The DERIVE function used is D219A(i,m,n,h) in c1999.mth."

"mm=m in the elliptic integral I(m) defined in equation 2.17."

mm:=[1,1,1,1]

[1,1,1,1]

"am is a in equation 2.17."

am:=VECTOR(-LIM(v12 SUB i SUB 1,t,0),i,1,4)

[SQRT(SQRT(6005)/1378+23377/474721)+373/689+#i*(SQRT(SQRT(6005)/1378-23377/47~
4721)-7/689),SQRT(SQRT(6005)/1378+23377/474721)+373/689+#i*(7/689-SQRT(SQRT(6~
005)/1378-23377/474721)),-SQRT(SQRT(6005)/1378+23377/474721)+373/689+#i*(SQRT~
(SQRT(6005)/1378-23377/474721)+7/689),-SQRT(SQRT(6005)/1378+23377/474721)+373~
/689-#i*(SQRT(SQRT(6005)/1378-23377/474721)+7/689)]

"bm is b in equation 2.17."

bm:=[-1,-1,-1,-1]

[-1,-1,-1,-1]

"In equation 2.19 set i=1, n=4, and h=4."

"The meaning of n and h are in defined in equation 2.17."

v13:=D219A(1,mm,4,4)

(16*b1^2*b2^2*b3^2*b4^2*AA(2*e1+e2+e3+e4)-4*b1^2*b2*b3*b4*(b2*(b3*d14+b4*d13)~
+b3*b4*d12)*AA(e1+e2+e3+e4)-2*b1^2*(b2^2*(3*b3^2*d14^2-2*b3*b4*d13*d14+3*b4^2~
*d13^2)-2*b2*b3*b4*d12*(b3*d14+b4*d13)+3*b3^2*b4^2*d12^2)*AA(e2+e3+e4)+d12*d1~
3*d14*(b2^2*(3*b3^2*d14^2-2*b3*b4*d13*d14+3*b4^2*d13^2)-2*b2*b3*b4*d12*(b3*d1~
4+b4*d13)+3*b3^2*b4^2*d12^2)*IA(-e1)-3*(b2^3*(b3*d14+b4*d13)*(b3^2*d14^2-2*b3~
*b4*d13*d14+b4^2*d13^2)-b2^2*b3*b4*d12*(b3^2*d14^2-2*b3*b4*d13*d14+b4^2*d13^2~
)-b2*b3^2*b4^2*d12^2*(b3*d14+b4*d13)+b3^3*b4^3*d12^3)*IA(e1)-2*b2*b3*b4*d12*d~
13*d14*(b2*(b3*d14+b4*d13)+b3*b4*d12)*IA(0))/(48*b1^3*b2^2*b3^2*b4^2)

"This is the symbolic result for a general case where [a1, a2, a3, a4]"

"and [b1, b2, b3, b4] are not known."

"The definition of d12 is in the equations 2.5, so d12=a1*b2-a2*b1."

"The DERIVE function D219C(v13,am,bm) substitutes the known values of am and ~
bm"

"into the expression v13. D219C(v,a,b) is in c1999.mth."

v14:=APPROX(D219C(v13,am,bm))

-AA(2*e1+e2+e3+e4)/3+410316503325994767334443019837177*AA(e1+e2+e3+e4)/379016~
2049822001357891628018744114+93410*AA(e2+e3+e4)/1424163-871740743946931484830~
55525093209*IA(-e1)/189613005024257193739742847756397599+360962*IA(e1)/327082~
769+45644084902941782834148180234469*IA(0)/2995452391143679918893525304481139~
8+170067266642230038208757659614900*#i*AA(e1+e2+e3+e4)/6945770629516546831603~
733513290637-150040412236723122079234056376066*#i*IA(0)/467341208540795325680~
88840565462875+29542178408859098661921356142802*#i*IA(-e1)/144171263579593605~
59333975090108167

"AA(e1+e2+e3+e4), AA(e2+e3+e4), and AA(2*e1+e2+e3+e4) are defined in equation~
 3.1."

"The DERIVE function AF(m,h,n,a,b,x,y) is in c1998.mth."

v15:=APPROX(AF([1,1,1,1],4,4,am,bm,1,0))

-3356419125071115776193163025630323/39634043647108890510864781516973327-#i/88~
45809875073415866358024419432147575757815309383336946182132291926

v16:=APPROX(AF([0,1,1,1],4,4,am,bm,1,0))

-1054818717697944464841667031558144/1138033847611793597106052032459343-481231~
794788388669044871219842429*#i/1338157964006727444110339662373480

v17:=APPROX(AF([2,1,1,1],4,4,am,bm,1,0))

-342732992658261421087414103787037/1768283676398045307074502575383106-3457522~
2640213610632975337162470*#i/5558218721047529300643628647692921

"Substitute the values of AA(e1+e2+e3+e4), AA(e2+e3+e4), and AA(2*e1+e2+e3+e4~
)"

"into the expression v14."

aaa:=[[e1+e2+e3+e4,v15],[e2+e3+e4,v16],[2*e1+e2+e3+e4,v17]]

[[e1+e2+e3+e4,-3356419125071115776193163025630323/396340436471088905108647815~
16973327-#i/88458098750734158663580244194321475757578153093833369461821322919~
26],[e2+e3+e4,-1054818717697944464841667031558144/113803384761179359710605203~
2459343-481231794788388669044871219842429*#i/13381579640067274441103396623734~
80],[2*e1+e2+e3+e4,-342732992658261421087414103787037/17682836763980453070745~
02575383106-34575222640213610632975337162470*#i/55582187210475293006436286476~
92921]]

AA(v):=(SELECT(IF(k SUB 1=v,1,0)=1,k,aaa)) SUB 1 SUB 2

APPROX(v14)

-0.0004597473384462261342413422130124312594647705268614391818799336994*IA(-e1~
)+0.001103579993234067307287593618237957377693595348032534236005565918*IA(e1)~
+0.001523779347583442126669540565907030340197307684489641415782174681*IA(0)-0~
.005353798479946789192685998100549122019216521277799806300854944187-0.0235874~
2887952136191306518829221358973364120729208992186182820527*#i-0.0032105110676~
03098767155636814318360138587410413111765150706050460*#i*IA(0)+0.002049103106~
636056374227004570560577340560268333471399221867933705*#i*IA(-e1)

"This is the results in terms of basic integrals I(-e1), I(e0), and I(e1)."

"e1 and e0 are defined in equation 2.18."

"Find the basic integrals in the table."

IS([IA(-e1),IA(0),IA(e1)],4)

[I428(a,b,x,y),I414(a,b,x,y),I429(a,b,x,y)]

"The DERIVE functions I428(a,b,x,y), I414(a,b,x,y), and I429(a,b,x,y) are in ~
c1998.mth."

"They are defined by equation 4.28, 4.14, and 4.29 in the paper."

"There is a revision to these equations in a latter paper"

"Reduction Theorems for Elliptic Integrands with Square Root of two Quadratic~
 Factors."

"Evaluate the basic integrals."

bi:=APPROX([I428(am,bm,1,0),I414(am,bm,1,0),I429(am,bm,1,0)])

[11220931622987071630752252037347423/542145479908018290171078323685157-522887~
00383515498226616152991954231*#i/1274657689271415753913537136464525,233244277~
0581540920367352720732871/193701306992308274126247860964963,89167717648651703~
75328086037331619/2252616204237948884922678304191219+548138824240512803451026~
079028593*#i/619712947005436412478777754565267]

"Now subsitute the values of the basic integrals into the expression v5 and"

"compare that to the check value."

iaa:=[[-e1,bi SUB 1],[0,bi SUB 2],[e1,bi SUB 3]]

[[-e1,11220931622987071630752252037347423/542145479908018290171078323685157-5~
2288700383515498226616152991954231*#i/1274657689271415753913537136464525],[0,~
2332442770581540920367352720732871/193701306992308274126247860964963],[e1,891~
6771764865170375328086037331619/2252616204237948884922678304191219+5481388242~
40512803451026079028593*#i/619712947005436412478777754565267]]

IA(v):=(SELECT(IF(k SUB 1=v,1,0)=1,k,iaa)) SUB 1 SUB 2

APPROX(km*v14-refint)

-1.485879888166669708140215357405205816568774728300889934160608663*10^(-62)-4~
.051277923685633254437564424508368724508509310805302476802034486*10^(-65)*#i

"This result matches the check value to 62 places."