{VERSION 3 0 "IBM INTEL NT" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "2 D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" 0 21 "" 0 1 0 0 0 1 0 0 0 0 2 0 0 0 0 }{CSTYLE "" -1 256 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 259 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "Cou rier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 " " 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 } {PSTYLE "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 6 6 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 4 4 0 0 0 0 0 0 -1 0 }{PSTYLE " Heading 3" 4 5 1 {CSTYLE "" -1 -1 "" 1 12 0 0 0 0 1 0 0 0 0 0 0 0 0 } 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 }{PSTYLE "" 2 6 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 2 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 12 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 } {PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }} {SECT 0 {PARA 18 "" 0 "" {TEXT -1 36 "Eigenstructure of the Kobbelt Sc heme" }}{PARA 19 "" 0 "" {TEXT -1 47 "Denis Zorin, Stanford University , February 1998" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 610 "In this worksheet, we explore the eigenstructure of the subdivision matrix for Kobbelt's subdivision scheme. We compute guar anteed intervals for the magnitude of the largest eigenvalue, and for mulas for computing corresponding eigenvector. The result is a file w ith three interval arithemtics C functions computing the magnitude of the largest eigenvalue for a large range of valences, and two functio ns to compute eigenvectors for the eigenvalue with the largest magnitu de. In addition, we estimate the range of eigenvalues for large vale nces, which allows us to analyze C1-continuity for all valences." }} {SECT 1 {PARA 3 "" 0 "" {TEXT -1 9 "Utilities" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 85 "with(share):\nwith(plots):\nreadshare(intpak, numeric s):\nreadlib(optimize):\nreadlib(C):" }}{PARA 6 "" 1 "" {TEXT -1 70 "S ee ?share and ?share,contents for information about the share library " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "This file contains some funct ions for code generation, and a fix for interval arithmetics " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "read `subdivmatrix-util.mpl`;" }}} }{SECT 0 {PARA 3 "" 0 "" {TEXT -1 41 "Subdivision matrix of the Kobbe lt scheme" }}{PARA 0 "" 0 "" {TEXT -1 159 "Define the blocks of the D FT-transformed subdivision matrix; perform some tests to check if the \+ matrix was defined correctly. We use the following parameters:" } {XPPEDIT 18 0 "alpha" "6#%&alphaG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "beta" "6#%%betaG" }{TEXT -1 125 " are the coeffficients of the 4-poi nt scheme (we consider the case when the coefficients are 9/16 and -1/ 16 respectively), " }{XPPEDIT 18 0 "c = cos(2*Pi/K)" "6#/%\"cG-%$cosG 6#*(\"\"#\"\"\"%#PiGF*%\"KG!\"\"" }{TEXT -1 7 ", and " }{XPPEDIT 18 0 "omega = exp(2*Pi/K)" "6#/%&omegaG-%$expG6#*(\"\"#\"\"\"%#PiGF*%\"KG !\"\"" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 373 "We test the cor rectness of the matrix in two ways: first, we compute a submatrix exp licitly for the case K = 4, and check if the matrices agree; second, compute the eigenvectors and eigenvalues for K =4; in this case, th e matrix has to have eigenvalue 1/2, and te corresponding complex eige nvector should be a part of a regular quadrilateral grid in the comple x plane." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1305 "Kobbelt := mat rix([[alpha+4*beta*d-beta*(1+2*c),\n4*d*beta^2/alpha-beta^2*(conjugate (omega)^2+2*c+1)/alpha, beta, 0, 0,\n0, 0, 0, 0, 0, 0, 0],\n[4*beta*al pha*d+alpha^2*(1+omega)-(1+omega)*alpha*beta,\n4*beta^2*d-beta^2*(1+2* c)+2*alpha*beta*c+alpha^2,\n(1+omega)*alpha*beta, beta^2*omega+alpha*b eta, beta^2,\nbeta^2*conjugate(omega)+alpha*beta, 0, 0, 0, 0, 0, 0], [ 1, 0, 0, 0, 0,\n0, 0, 0, 0, 0, 0, 0], [alpha, alpha+beta*conjugate(ome ga), 0, 0, 0,\nbeta, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n[alpha*omega, alpha+beta*omega, 0, beta, 0, 0, 0, 0, 0, 0, 0, 0],\n[alpha, 0, alpha, 0, 0, 0, beta, 0, 0, 0, 0, 0],\n[beta^2*conjug ate(omega)+alpha^2+alpha*beta*omega,\nalpha*beta*conjugate(omega)+alph a^2, beta^2*omega+alpha^2, alpha^2,\nalpha*beta, alpha*beta*(1+conjuga te(omega)), alpha*beta, alpha*beta,\nbeta^2, 0, 0, beta^2*conjugate(om ega)], [beta*omega, alpha, 0, alpha,\n0, 0, 0, beta, 0, 0, 0, 0], [(1+ omega)*alpha*beta, alpha^2,\n(1+omega)*alpha*beta, alpha^2, alpha^2, a lpha^2, beta^2*(1+omega),\nalpha*beta, alpha*beta, beta^2, alpha*beta, alpha*beta], [beta, alpha,\n0, 0, 0, alpha, 0, 0, 0, 0, 0, beta],\n[a lpha*beta+alpha^2*omega+beta^2*omega^2, alpha^2+alpha*beta*omega,\nbet a^2+alpha^2*omega, (1+omega)*alpha*beta, alpha*beta, alpha^2,\nalpha*b eta*omega, beta^2*omega, 0,0,beta^2, alpha*beta]]);\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(KobbeltG-%'MATRIXG6#7.7.,(%&alphaG\"\"\"*&%%bet aGF,%\"dGF,\"\"%*&F.F,,&F,F,%#c|irG\"\"#F,!\"\",&*(F/F,F.F4F+F5F0*(F.F 4,(*$-%*conjugateG6#%&omegaGF4F,F3F4F,F,F,F+F5F5F.\"\"!F?F?F?F?F?F?F?F ?7.,(*(F.F,F+F,F/F,F0*&F+F4,&F,F,F>F,F,F,*(FDF,F+F,F.F,F5,**&F.F4F/F,F 0*&F.F4F2F,F5*(F+F,F.F,F3F,F4*$F+F4F,FE,&*&F.F4F>F,F,*&F+F,F.F,F,*$F.F 4,&*&F.F4F;F,F,FMF,F?F?F?F?F?F?7.F,F?F?F?F?F?F?F?F?F?F?F?7.F+,&F+F,*&F .F,F;F,F,F?F?F?F.F?F?F?F?F?F?7.F?F,F?F?F?F?F?F?F?F?F?F?7.*&F+F,F>F,,&F +F,*&F.F,F>F,F,F?F.F?F?F?F?F?F?F?F?7.F+F?F+F?F?F?F.F?F?F?F?F?7.,(FPF,F JF,*(F+F,F.F,F>F,F,,&*(F+F,F.F,F;F,F,FJF,,&FLF,FJF,FJFM*(F+F,F.F,,&F;F ,F,F,F,FMFMFNF?F?FP7.FYF+F?F+F?F?F?F.F?F?F?F?7.FEFJFEFJFJFJ*&F.F4FDF,F MFMFNFMFM7.F.F+F?F?F?F+F?F?F?F?F?F.7.,(FMF,*&F+F4F>F,F,*&F.F4F>F4F,,&F JF,FgnF,,&FNF,FcoF,FEFMFJFgnFLF?F?FNFM" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "Kobconst := \{ alpha = 9/16, beta = -1/16 \}; " }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%)KobconstG<$/%&alphaG#\"\"*\"#;/%%be taG#!\"\"F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "KobExpanded \+ := subs(Kobconst , evalm(Kobbelt));" }}{PARA 11 "" 0 "" {TEXT -1 0 "" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,KobExpandedG-%'MATRIXG6#7.7.,(#\" \"&\"\")\"\"\"%\"dG#!\"\"\"\"%%#c|irG#F.F-,*F/#F.\"#O*$-%*conjugateG6# %&omegaG\"\"##F1\"$W\"F3#F1\"#sF>F.#F1\"#;\"\"!FDFDFDFDFDFDFDFD7.,(F/# !\"*\"#k#\"#X\"$G\"F.F " 0 "" {MPLTEXT 1 0 137 " KobRegular := map( unapply( s ubs( c = 0,x), x), map( simplify, map( evalc, subs( \{ d = 0, omega = I, c = 0 \}, evalm(KobExpanded) )))):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "Manually computed matr ix for the regular case" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 284 "KobRegu larManual := matrix( [[ 9/16 - I^2/16, 0,-1/16,0,0,0], [ (81/256)*(I+1 ) - (9/256)*(I^2 - I), 81/256 + (1/256)*I^2, (-9/256)*(I+1), -9/256 + I/256, 1/256, -9/256 - (1/256)*I], \n[1,0,0,0,0,0], [9/16,9/16 +I/16, \+ 0,0,0, -1/16],[0,1,0,0,0,0],[ 9*I/16, 9/16 - I/16, 0,-1/16,0, 0] ]);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%1KobRegularManualG-%'MATRIXG6#7(7( #\"\"&\"\")\"\"!#!\"\"\"#;F-F-F-7(,&#\"#X\"$G\"\"\"\"%\"IGF3#F+F0,&#! \"*\"$c#F6F7F:,&F:F6F7#F6F,&F:F6F7#F/F<7(F6F-F-F-F-F-7(#\"\"*F0,&FC F6F7#F6F0F-F-F-F.7(F-F6F-F-F-F-7(,$F7FC,&FCF6F7F.F-F.F-F-" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 39 "Check agreement with the regular case ." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "norm( evalm(submatrix( KobRegula r, 1..6,1..6) - KobRegularManual));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #\"\"!" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 66 "Check if the eigenvector for the eigenvalue 1/2 is a regu lar grid:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "eigenvects( KobRegular );" }}{PARA 12 "" 0 "" {TEXT -1 0 "" }}{PARA 12 "" 1 "" {XPPMATH 20 " 6*7%#\"\"\"\"$c#F%<#-%'VECTORG6#7.\"\"!F,F,F,F,F,F,F,F,F%F,F,7%#!\"\" \"$G\"F%<#-F)6#7.F,F,F,F,F,F,F,#F%\"\")F%,&#\"#FF6F%%\"IG#!#FF6,$F:F/, $F:#F/F67%#F%\"\"#F%<#-F)6#7.F%,&F%F%F:F%FB,&FBF%F:F%,&FBF%F:FB,&F%F%F :FB\"\"$,&FKF%F:F%,&FKF%F:FB,&FKF%F:FK,&FBF%F:FK,&F%F%F:FK7%#F/\"#KF%< #-F)6#7.F,F,F,F,F,F,F,,$F:FAF:,&#FKFBF%F:FZF%FA7%#F/\"#;F%<#-F)6#7.F,F ,F,F,F,F,F%F%F%FGF:F:7%#F/\"#kFB<#-F)6#7.F,F,F,F,F,F,F,F=,$F:!\"%,&\" \"*F%F:!\"*\"\"%F%7%#F%FSFB<#-F)6#7.F,F%F,,&\"\"'F%F:!\"#FS,&F`pF%F:FB F,,&\"#=F%F:Fgo,&\"#!*F%F:!#=\"$V#,&FfpF%F:Fdp,&FdpF%F:Ffo7%F5FK<$-F)6 #7.F5F,F%,&#FKFhoF%F:#!\"$F6F,,&FbqF%F:FaqF8,&FKF%F:F/,&#\"#:F6F%F:#! \"&FhoF,,&FiqF%F:Fgq,&F/F%F:FK-F)6#7.F,F5F,,&#FKF6F%F:F5F%,&FarF%F:F?F ,,&FaqF%F:Far,&FgqF%F:FarF8,&FgqF%F:FbqF`q" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "Eigenvalues of the 0-th block for all valences:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 142 "KobZeroBlock := map( unapply( subs( c = \+ 1,'x'), 'x'), map( simplify, map( evalc, subs( \{ d = 1, omega = 1, c = 1 \}, evalm(KobExpanded) )))):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 182 "The largest eigenvalue is 1/4; we will see that the largest e igenvalue of one of the other blocks is always greater than 1/4, and t he dominant eigenvalue is never in the 0-th block." }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 24 "eigenvals(KobZeroBlock);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6.#\"\"\"\"$c##!\"\"\"#K#F'\"$G\"#F'\"#;#F'\"#kF-#F$\"\"% F/#F$F,F1F1F1" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 51 "Characteristic polynomial of the subdivision matrix" }}{PARA 0 "" 0 "" {TEXT -1 147 "Compute and factor the characteristic polynomial of the blocks of DFT -transformed subdivision matrix. The resulting polynomial is parameter ized by " }{XPPEDIT 18 0 "c = cos(2*Pi/K)" "6#/%\"cG-%$cosG6#*(\"\"#\" \"\"%#PiGF*%\"KG!\"\"" }{TEXT -1 384 ". The polynomial has a number o f small roots, which do not depend on c, and a factor of degree 6. F or illustrative purposes, and to guide us in the subsequent derivation s, we compute all the roots of the polynomial numerically, and plot th e magnitudes of the roots. Of course, neither the plot nor the compute d values cannot be used in the analysis without additional verificatio n." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 84 "We have already computed eig envalues for the 0th block, and we can assume that d = 0" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 235 "KobCharpolynom := subs( \{ s^3 = s*(1-c^2),s ^2 = 1 - c^2,s^4 = (1-c^2)^2\}, factor( collect(charpoly( map( simplif y, map( evalc, subs( s^2 = 1-c^2, map( simplify, subs( \{ d = 0, omega = c + I*s\}, eval(KobExpanded)))))),lambda),lambda))):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "KobCharpolynom := factor( map( simp lify, KobCharpolynom));" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/KobCharp olynomG,$*.,&%'lambdaG\"#K\"\"\"F*F*,&F(\"$G\"F*F*F*,&F*F*F(\"#k\"\"#, " 0 "" {MPLTEXT 1 0 106 "Ko bFactor6 := KobCharpolynom/( (32*lambda+1)*(128*lambda+1)*(1+64*lambda )^2*(256*lambda-1)*(16*lambda+1)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "KobFactor6 := collect(KobFactor6/lcoeff(KobFactor6,la mbda),lambda);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%+KobFactor6G,0*$%' lambdaG\"\"'\"\"\"*&,&%#c|irG#!\"$\"#k#!#:\"#;F)F)F'\"\"&F)*&,&F,#!\"* \"%C5#\"$(HF8F)F)F'\"\"%F)*&,(#!$N$\"%#>)F)*$F,\"\"##!\"(\"&%Q;F,#\"#@ \"%'4%F)F'\"\"$F)*&,&F,#F7FE#\"$$=\"&Ob'F)F)F'FBF)*&,&#!#X\"')GC&F)F,# \"\"*FTF)F'F)F)#F)\"(w&[5F)" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "Co mpute the eigenvalues; execution of this statement may take a while. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "EigenvalsList := seq([solve( su bs( c = evalf( (n+1e-10)/100),KobFactor6))], n = -100..100): " }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 66 "Convert the lists of eigenvalues t o the form suitable for plotting" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 110 " EigenvalsPlotLists := seq( [ seq( [-1 + (i-1)/100, abs(op(j, op( i,[EigenvalsList])))], i = 1..201)] ,j=1..6):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 124 "Plot of the magnitudes of eigenvalues as functions of c; to see approximately the magnitudes of the eigenvalues for a block " }{XPPEDIT 18 0 "m" "6#%\"mG" }{TEXT -1 39 " of the subdivision matr ix for valence " }{XPPEDIT 18 0 "K" "6#%\"KG" }{TEXT -1 27 ", draw a \+ vertical line at " }{XPPEDIT 18 0 "c = cos(2*Pi*m/K)" "6#/%\"cG-%$cosG 6#**\"\"#\"\"\"%#PiGF*%\"mGF*%\"KG!\"\"" }{TEXT -1 55 "; and find whe re it intersects the curves in the plot." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "display(seq(plot(op(i,[EigenvalsPlotLists]),color=bl ack), i = 1..6),color=black, axesfont=[TIMES,ITALIC,10], labels=[``,`` ]);" }}{PARA 13 "" 1 "" {GLPLOT2D 553 305 305 {PLOTDATA 2 "6+-%'CURVES G6$7ew7$$!\"\"\"\"!$\"1+++)***\\i:!#<7$$!1+++++++**!#;$\"1+++O?0o:F-7$ $!1+++++++)*F1$\"1+++nTlt:F-7$$!1+++++++(*F1$\"1+++9sIz:F-7$$!1+++++++ '*F1$\"1+++U?,&e\"F-7$$!1+++++++&*F1$\"1+++?&p2f\"F-7$$!1************* R*F1$\"1+++Y0e'f\"F-7$$!1+++++++$*F1$\"1+++UgW-;F-7$$!1+++++++#*F1$\"1 +++[pO3;F-7$$!1+++++++\"*F1$\"1+++\\UM9;F-7$$!1+++++++!*F1$\"1+++O*y.i \"F-7$$!1+++++++*)F1$\"1+++P?ZE;F-7$$!1+++++++))F1$\"1+++,YiK;F-7$$!1+ ++++++()F1$\"1+++Dx$)Q;F-7$$!1+++++++')F1$\"1+++=D6X;F-7$$!1+++++++&)F 1$\"1+++K,X^;F-7$$!1+++++++%)F1$\"1+++S<&yl\"F-7$$!1+++++++$)F1$\"1+++ u&=Vm\"F-7$$!1+++++++#)F1$\"1+++!)=&3n\"F-7$$!1,++++++\")F1$\"1+++[HXx ;F-7$$!1+++++++!)F1$\"1+++7J7%o\"F-7$$!1+++++++zF1$\"1+++XP'3p\"F-7$$! 1+++++++yF1$\"1+++iin(p\"F-7$$!1+++++++xF1$\"1+++C@c/F-7$$!1+++++++^F1$\"1 +++T%yb\">F-7$$!1+++++++]F1$\"1+++F8CD>F-7$$!1+++++++\\F1$\"1+++c,0N>F -7$$!1+++++++[F1$\"1+++![4]%>F-7$$!1+++++++ZF1$\"1+++&3C^&>F-7$$!1++++ +++YF1$\"1+++0!*Rl>F-7$$!1+++++++XF1$\"1+++Y&Re(>F-7$$!1+++++++WF1$\"1 +++48X')>F-7$$!1+++++++VF1$\"1+++.-C(*>F-7$$!1+++++++UF1$\"1+++zC@3?F- 7$$!1+++++++TF1$\"1+++aZP>?F-7$$!1+++++++SF1$\"1+++[StI?F-7$$!1+++++++ RF1$\"1+++5yHU?F-7$$!1+++++++QF1$\"1+++cR2a?F-7$$!1+++++++PF1$\"1+++M4 2m?F-7$$!1+++++++OF1$\"1+++RxHy?F-7$$!1+++++++NF1$\"1+++*)Rw!4#F-7$$!1 +++++++MF1$\"1+++%)*zM5#F-7$$!1+++++++LF1$\"1+++snX;@F-7$$!1+++++++KF1 $\"1+++FiqH@F-7$$!1+++++++JF1$\"1+++V6CV@F-7$$!1+++++++IF1$\"1+++S`2d@ F-7$$!1+++++++HF1$\"1+++vPAr@F-7$$!1+++++++GF1$\"1+++%o-d=#F-7$$!1++++ +++FF1$\"1+++Y(H0?#F-7$$!1+++++++EF1$\"1+++fUs:AF-7$$!1+++++++DF1$\"1+ ++otIJAF-7$$!1+++++++CF1$\"1+++ABIZAF-7$$!1+++++++BF1$\"1+++uZtjAF-7$$ !1+++++++AF1$\"1+++pJj!G#F-7$$!1+++++++@F1$\"1+++d\"H!)H#F-7$$!1++++++ +?F1$\"1+++j\"efJ#F-7$$!1+++++++>F1$\"1+++H+YMBF-7$$!1+++++++=F1$\"1++ +Q)zNN#F-7$$!1+++++++alDF-7$$FjqF-$\"1 +++[.&ef#F-7$$!1,++++++qF-$\"1+++M=QGEF-7$$F^xF-$\"1+++aGF-7$$FhglF-$\"1+++$G5C$HF-7$F*$\"1+++PU([7$F-7$$\"1+++++ ++5F-$\"1+++q<2HJF-7$$\"1+++++++?F-$\"1+++^^=LJF-7$$\"1+++++++IF-$\"1+ ++B+MPJF-7$$\"1+++++++SF-$\"1+++qi`TJF-7$$\"1+++++++]F-$\"1+++nPxXJF-7 $$\"1+++++++gF-$\"1+++yB0]JF-7$$\"1,++++++qF-$\"1+++k>PaJF-7$$\"1+++++ ++!)F-$\"1+++oBteJF-7$$\"1+++++++!*F-$\"1+++PM8jJF-7$$FejlF1$\"1+++$) \\dnJF-7$$\"1+++++++6F1$\"1+++Ko0sJF-7$$\"1+++++++7F1$\"1+++v(yl<$F-7$ $\"1+++++++8F1$\"1+++/19\"=$F-7$$\"1+++++++9F1$\"1+++)3Ud=$F-7$$\"1+++ ++++:F1$\"1+++!*HQ!>$F-7$$\"1+++++++;F1$\"1+++JI1&>$F-7$$\"1+++++++y*>$F-7$$\"1+++++++=F1$\"1+++)RRX?$F-7$$\"1+++++++>F1$\"1+++' 4N$4KF-7$$FjjlF1$\"1+++u'oT@$F-7$$\"1+++++++@F1$\"1+++d(R!>KF-7$$\"1++ +++++AF1$\"1+++dz%RA$F-7$$\"1+++++++BF1$\"1+++PG*)GKF-7$$\"1+++++++CF1 $\"1+++SR(QB$F-7$$\"1+++++++DF1$\"1+++)y!*)QKF-7$$\"1+++++++EF1$\"1+++ ZG%RC$F-7$$\"1+++++++FF1$\"1+++p&H!\\KF-7$$\"1+++++++GF1$\"1+++j.:aKF- 7$$\"1+++++++HF1$\"1+++4YIfKF-7$$F_[mF1$\"1+++P;\\kKF-7$$\"1+++++++JF1 $\"1+++V2rpKF-7$$\"1+++++++KF1$\"1+++z6'\\F$F-7$$\"1+++++++LF1$\"1+++] @C!G$F-7$$\"1+++++++MF1$\"1+++?Gb&G$F-7$$\"1+++++++NF1$\"1+++1B*3H$F-7 $$\"1+++++++OF1$\"1+++i'fiH$F-7$$\"1+++++++PF1$\"1+++3Rl,LF-7$$\"1++++ +++QF1$\"1+++7S22LF-7$$\"1+++++++RF1$\"1+++u)=DJ$F-7$$Fd[mF1$\"1+++Xt) zJ$F-7$$\"1+++++++TF1$\"1+++F#yMK$F-7$$\"1+++++++UF1$\"1+++Q-**GLF-7$$ \"1+++++++VF1$\"1+++g?_MLF-7$$\"1+++++++WF1$\"1+++4B2SLF-7$$\"1+++++++ XF1$\"1+++J&RcM$F-7$$\"1+++++++YF1$\"1+++>AA^LF-7$$\"1+++++++ZF1$\"1++ +%z=oN$F-7$$\"1+++++++[F1$\"1+++:wUiLF-7$$\"1+++++++\\F1$\"1+++xp/oLF- 7$$Fi[mF1$\"1+++*4vOP$F-7$$\"1+++++++^F1$\"1+++c,JzLF-7$$\"1+++++++_F1 $\"1+++O-&\\Q$F-7$$\"1+++++++`F1$\"1+++uLf!R$F-7$$\"1+++++++aF1$\"1+++ UvB'R$F-7$$\"1+++++++bF1$\"1+++\\1)=S$F-7$$\"1,++++++cF1$\"1+++R0_2MF- 7$$\"1+++++++dF1$\"1+++/]:8MF-7$$\"1+++++++eF1$\"1+++fG+IMF-7$$\"1+++++++hF1$\"1+++> BfNMF-7$$\"1+++++++iF1$\"1+++KX;TMF-7$$\"1+++++++jF1$\"1+++QprYMF-7$$ \"1+++++++kF1$\"1+++**pC_MF-7$$\"1+++++++lF1$\"1+++L@vdMF-7$$\"1++++++ +mF1$\"1+++Z(HKY$F-7$$\"1+++++++nF1$\"1+++4snoMF-7$$\"1+++++++oF1$\"1+ ++%)=4uMF-7$$\"1**************oF1$\"1+++H6ZzMF-7$$\"1+++++++qF1$\"1+++ \"H7[[$F-7$$\"1+++++++rF1$\"1+++JF6!\\$F-7$$\"1+++++++sF1$\"1+++;)p`\\ $F-7$$\"1+++++++tF1$\"1+++L4e+NF-7$$\"1+++++++uF1$\"1+++#[Vd]$F-7$$\"1 +++++++vF1$\"1+++9\\&3^$F-7$$\"1+++++++wF1$\"1+++2F\"f^$F-7$$\"1++++++ +xF1$\"1+++*Q94_$F-7$$\"1+++++++yF1$\"1+++Xv&e_$F-7$$\"1+++++++zF1$\"1 +++?)R2`$F-7$$Fh\\mF1$\"1+++8*eb`$F-7$$\"1,++++++\")F1$\"1+++6EJSNF-7$ $\"1+++++++#)F1$\"1+++v()*\\a$F-7$$\"1+++++++$)F1$\"1+++\\`h\\NF-7$$\" 1+++++++%)F1$\"1+++o.;aNF-7$$\"1+++++++&)F1$\"1+++i>jeNF-7$$\"1+++++++ ')F1$\"1+++o$GIc$F-7$$\"1+++++++()F1$\"1+++8zMnNF-7$$\"1+++++++))F1$\" 1+++R!*erNF-7$$\"1+++++++*)F1$\"1+++0.vvNF-7$$F]]mF1$\"1+++d.$)zNF-7$$ \"1+++++++\"*F1$\"1+++xz#Qe$F-7$$\"1+++++++#*F1$\"1+++i?u(e$F-7$$\"1++ +++++$*F1$\"1+++9;d\"f$F-7$$\"1*************R*F1$\"1+++ddJ&f$F-7$$\"1+ ++++++&*F1$\"1+++LP(*)f$F-7$$\"1+++++++'*F1$\"1+++(*[a-OF-7$$\"1++++++ +(*F1$\"1+++8(Ggg$F-7$$\"1+++++++)*F1$\"1+++hZU4OF-7$$\"1+++++++**F1$ \"1+++KFt7OF-7$$\"\"\"F*$\"1+++;C&fh$F--%'COLOURG6&%$RGBGF*F*F*-F$6$7e w7$F($\"1+++$*****\\iF-7$F/$\"1+++w19)*fF-7$F5$\"1+++5%yz*eF-7$F:$\"1+ ++%GKD#eF-7$F?$\"1+++AnyfdF-7$FD$\"1+++!\\*40dF-7$FI$\"1+++b[4ccF-7$FN $\"1+++;WO6cF-7$FS$\"1+++bA**pbF-7$FX$\"1+++EIMJbF-7$Fgn$\"1+++j^&\\\\ &F-7$F\\o$\"1+++!\\![gaF-7$Fao$\"1+++(**[wU&F-7$Ffo$\"1+++hdC'R&F-7$F[ p$\"1+++iq4m`F-7$F`p$\"1+++8&fqL&F-7$Fep$\"1+++'\\8!4`F-7$Fjp$\"1+++([ d=G&F-7$F_q$\"1+++s\\]b_F-7$Fdq$\"1+++M6))H_F-7$Fiq$\"1+++#p?\\?&F-7$F ^r$\"1+++:kc!=&F-7$Fcr$\"1+++cwwc^F-7$Fhr$\"1+++'QzM8&F-7$F]s$\"1+++)= h16&F-7$Fbs$\"1+++gnF)3&F-7$Fgs$\"1+++'=$Hm]F-7$F\\t$\"1+++c0oW]F-7$Fa t$\"1+++;:TB]F-7$Fft$\"1+++K4Y-]F-7$F[u$\"1+++-d!=)\\F-7$F`u$\"1+++OVU h\\F-7$Feu$\"1+++NpHT\\F-7$Fju$\"1+++Q\\S@\\F-7$F_v$\"1+++%)4t,\\F-7$F dv$\"1+++7(eA)[F-7$Fiv$\"1+++=F(H'[F-7$F^w$\"1++++%eQ%[F-7$Fcw$\"1+++Z >!\\#[F-7$Fhw$\"1+++X,41[F-7$F]x$\"1+++k.T(y%F-7$Fbx$\"1+++:0&)oZF-7$F gx$\"1+++\\*)R]ZF-7$F\\y$\"1+++.W/KZF-7$Fay$\"1+++afx8ZF-7$Ffy$\"1+++) )He&p%F-7$F[z$\"1+++$>bun%F-7$F`z$\"1+++cBQfYF-7$Fez$\"1+++2XNTYF-7$Fj z$\"1+++U=OBYF-7$F_[l$\"1+++8YR0YF-7$Fd[l$\"1+++%=Vue%F-7$Fi[l$\"1+++' )y\\pXF-7$F^\\l$\"1+++$4\\:b%F-7$Fc\\l$\"1+++YreLXF-7$Fh\\l$\"1+++@Bg: XF-7$F]]l$\"1+++bZe(\\%F-7$Fb]l$\"1+++pW_zWF-7$Fg]l$\"1+++P8ThWF-7$F\\ ^l$\"1+++\"*\\BVWF-7$Fa^l$\"1+++zZ)\\U%F-7$Ff^l$\"1+++w(\\mS%F-7$F[_l$ \"1+++u(=#)Q%F-7$F`_l$\"1+++/+opVF-7$Fe_l$\"1+++x8-^VF-7$Fj_l$\"1+++i- BKVF-7$F_`l$\"1+++:MH8VF-7$Fd`l$\"1+++cp>%H%F-7$Fi`l$\"1+++Oi#\\F%F-7$ F^al$\"1+++*ylaD%F-7$Fcal$\"1+++o\"*zNUF-7$Fhal$\"1+++z)3f@%F-7$F]bl$ \"1+++#=wd>%F-7$Fbbl$\"1+++v4QvTF-7$Fgbl$\"1+++*f,Z:%F-7$F\\cl$\"1+++n XrLTF-7$Facl$\"1+++[VR7TF-7$Ffcl$\"1+++cIr!4%F-7$F[dl$\"1+++-+koSF-7$F `dl$\"1+++(RTh/%F-7$Fedl$\"1+++a'zJ-%F-7$Fjdl$\"1+++9Gr**RF-7$F_el$\"1 +++HPpvRF-7$Fdel$\"1+++V*o5&RF-7$Fiel$\"1+++euxDRF-7$F^fl$\"1+++`!\\(* *QF-7$Fcfl$\"1+++H?!H(QF-7$Fhfl$\"1+++=/9XQF-7$F]gl$\"1+++-*\\j\"QF-7$ Fbgl$\"1+++4DR'y$F-7$Fggl$\"1+++%y)4bPF-7$F\\hl$\"1+++ilDAPF-7$F`hl$\" 1+++&z$f(o$F-7$Fdhl$\"1+++p8v]OF-7$Fihl$\"1+++6wB6OF-7$F]il$\"1+++6pMo NF-7$Fail$\"1+++5\"**4_$F-7$Feil$\"1+++.0QnMF-7$Fiil$\"1+++TD'RS$F-7$F ]jl$\"1+++\"Qa;K$F-7$F*$\"1+++yd7DJF-FcjlFhjlF][mFb[mFg[mF\\\\mFa\\mFf \\mF[]mF`]mFd]mFi]mF^^mFc^mFh^mF]_mFb_mFg_mF\\`mFa`mFe`mFj`mF_amFdamFi amF^bmFcbmFhbmF]cmFbcmFfcmF[dmF`dmFedmFjdmF_emFdemFiemF^fmFcfmFgfmF\\g mFagmFfgmF[hmF`hmFehmFjhmF_imFdimFhimF]jmFbjmFgjmF\\[nFa[nFf[nF[\\nF` \\nFe\\nFi\\nF^]nFc]nFh]nF]^nFb^nFg^nF\\_nFa_nFf_nF[`nF``nFe`nFj`nF_an FdanFianF^bnFcbnFhbnF\\cnFacnFfcnF[dnF`dnFednFjdnF_enFdenFienF]fnFbfnF gfnF\\gnFagnFfgnF[hnF`hnFehnFjhnF_in-F$6$7ewFfin7$F/$\"1+++?\\8>lF-7$F 5$\"1,++v(fmj'F-7$F:$\"1+++5lbHnF-7$F?$\"1+++f>%)4oF-7$FD$\"1+++z%f@)o F-7$FI$\"1******3g))[pF-7$FN$\"1+++d6V6qF-7$FS$\"1+++$=722(F-7$FX$\"1+ ++;dOFrF-7$Fgn$\"1+++!fa==(F-7$F\\o$\"1+++$RGXB(F-7$Fao$\"1+++,&ecG(F- 7$Ffo$\"1+++z8YNtF-7$F[p$\"1+++=@6%Q(F-7$F`p$\"1+++ucvJuF-7$Fep$\"1+++ IK^yuF-7$Fjp$\"1+++%)z[CvF-7$F_q$\"1+++!*zwpvF-7$Fdq$\"1+++^(HWh(F-7$F iq$\"1*****zRS&ewF-7$F^r$\"1+++?*e@q(F-7$Fcr$\"1+++nyLXxF-7$Fhr$\"1+++ ZS7)y(F-7$F]s$\"1+++L+cIyF-7$Fbs$\"1*****z1%osyF-7$Fgs$\"1+++t6`9zF-7$ F\\t$\"1,++xL8czF-7$Fat$\"1*****\\J?v*zF-7$Fft$\"1+++*[>(Q!)F-7$F[u$\" 1+++'Hc(z!)F-7$F`u$\"1+++^Zl?\")F-7$Feu$\"1*****HUP9;)F-7$Fju$\"1***** RUD@?)F-7$F_v$\"1+++*zQFC)F-7$Fdv$\"1******=oH$G)F-7$Fiv$\"1+++Uy\"QK) F-7$F^w$\"1,++/'>VO)F-7$Fcw$\"1+++u*=[S)F-7$Fhw$\"1+++VDLX%)F-7$F]x$\" 1,++ak(e[)F-7$Fbx$\"1+++H!)*)F-7$F^\\l$\"1***** *zegA!*F-7$Fc\\l$\"1+++`J;l!*F-7$Fh\\l$\"1+++@-)z5*F-7$F]]l$\"1+++#Qu5 :*F-7$Fb]l$\"1,++UKY%>*F-7$Fg]l$\"1+++sZ;Q#*F-7$F\\^l$\"1+++ky>#G*F-7$ Fa^l$\"1,++xCeE$*F-7$Ff^l$\"1,++m\"R8P*F-7$F[_l$\"1*****>**)[;%*F-7$F` _l$\"1+++@]0i%*F-7$Fe_l$\"1*****z*413&*F-7$Fj_l$\"1+++l=`a&*F-7$F_`l$ \"1+++EU\\,'*F-7$Fd`l$\"1+++4h(*['*F-7$Fi`l$\"1+++\")y+(p*F-7$F^al$\"1 *****\\W@cu*F-7$Fcal$\"1+++97&[z*F-7$Fhal$\"1+++5RtW)*F-7$F]bl$\"1+++[ *4`*)*F-7$Fbbl$\"1*****H3Am%**F-7$Fgbl$\"1+++hor)***F-7$F\\cl$\"1+++[X ;05F17$Facl$\"1+++zja55F17$Ffcl$\"1+++7L-;5F17$F[dl$\"1+++W@g@5F17$F`d l$\"1+++G.HF5F17$Fedl$\"1+++-j4L5F17$Fjdl$\"1+++g%H!R5F17$F_el$\"1+++G .5X5F17$Fdel$\"1+++O4K^5F17$Fiel$\"1+++R^qd5F17$F^fl$\"1+++)yoU1\"F17$ Fcfl$\"1+++h/.r5F17$Fhfl$\"1+++)z6!y5F17$F]gl$\"1+++'**Q_3\"F17$Fbgl$ \"1+++HLu#4\"F17$Fggl$\"1+++^Lc+6F17$F\\hl$\"1+++0su36F17$F`hl$\"1+++m nN<6F17$Fdhl$\"1+++0NZE6F17$Fihl$\"1+++,)3i8\"F17$F]il$\"1+++<;sY6F17$ Fail$\"1+++RaDe6F17$Feil$\"1+++%*e@r6F17$Fiil$\"1+++1tQ'=\"F17$F]jl$\" 1+++7Px07F17$F*$\"1+++/`X[7F17$Fdjl$\"1+++lP=O7F17$Fijl$\"1+++T7ZB7F17 $F^[m$\"1+++CC]67F17$Fc[m$\"1+++:))3+7F17$Fh[m$\"1+++oX6*=\"F17$F]\\m$ \"1+++!R.&y6F17$Fb\\m$\"1+++a2?o6F17$Fg\\m$\"1+++^p;e6F17$F\\]m$\"1+++ Y1P[6F17$Fa]m$\"1+++*H(yQ6F17$Fe]m$\"1+++_pRH6F17$Fj]m$\"1+++tO=?6F17$ F_^m$\"1+++PN866F17$Fd^m$\"1+++v\\B-6F17$Fi^m$\"1+++>zZ$4\"F17$F^_m$\" 1+++mT&[3\"F17$Fc_m$\"1+++jgNw5F17$Fh_m$\"1+++dq(z1\"F17$F]`m$\"1+++)G 6(f5F17$Fb`m$\"1+++TPb^5F17$Ff`m$\"1+++&z*\\V5F17$F[am$\"1+++>_aN5F17$ F`am$\"1+++;joF5F17$Feam$\"1+++7)>*>5F17$Fjam$\"1+++l/5F17$Fdbm$\"1+++4OXr**F-7$Fibm$\"1+++!>0s*)*F-7$F^cm$\"1+++[JvB) *F-7$Fccm$\"1+++)Ry5v*F-7$Fgcm$\"1+++LJ;z'*F-7$F\\dm$\"1+++<=*zg*F-7$F adm$\"1,++`&\\v`*F-7$Ffdm$\"1+++vG#yY*F-7$F[em$\"1+++f#*z)R*F-7$F`em$ \"1+++\\\"o/L*F-7$Feem$\"1+++;)=GE*F-7$Fjem$\"1+++><%e>*F-7$F_fm$\"1** ***fbG&H\"*F-7$Fdfm$\"1+++A;(Q1*F-7$Fhfm$\"1,++7R'))**)F-7$F]gm$\"1+++ =$*\\M*)F-7$Fbgm$\"1,++.>xq))F-7$Fggm$\"1+++(*on2))F-7$F\\hm$\"1,++;)4 _u)F-7$Fahm$\"1+++nkO$o)F-7$Ffhm$\"1+++QM9A')F-7$F[im$\"1+++6!Q:c)F-7$ F`im$\"1+++Iva,&)F-7$Feim$\"1+++_+HyF-7$F_]n$\"1+++q]7xxF-7$Fd]n$\"1+++%Rjcs(F-7$Fi]n$\"1,++mQ\"[n(F-7 $F^^n$\"1+++xvdCwF-7$Fc^n$\"1,++X`&\\d(F-7$Fh^n$\"1+++q#[f_(F-7$F]_n$ \"1+++,tbxuF-7$Fb_n$\"1+++@KyHuF-7$Fg_n$\"1+++Re#oF-7$F\\dn$\"1+++dkI(y 'F-7$Fadn$\"1++++yP\\nF-7$Ffdn$\"1++++*\\?r'F-7$F[en$\"1+++f#>`n'F-7$F `en$\"1+++j?=RmF-7$Feen$\"1*****>1MOg'F-7$Fjen$\"1+++?3nolF-7$F^fn$\"1 +++3vGMlF-7$Fcfn$\"1,++mI@F17$Fgn$\"1+++ F`f8@F17$F\\o$\"1+++3'\\t4#F17$Fao$\"1+++:$3=3#F17$Ffo$\"1+++'fzo1#F17 $F[p$\"1+++(f*[_?F17$F`p$\"1+++xydQ?F17$Fep$\"1+++rU4D?F17$Fjp$\"1+++s l*>,#F17$F_q$\"1+++.*[#**>F17$Fdq$\"1+++^/#o)>F17$Fiq$\"1+++?Xou>F17$F ^r$\"1+++2y\"G'>F17$Fcr$\"1+++:)*>^>F17$Fhr$\"1+++@C\")R>F17$F]s$\"1++ +x%R'G>F17$Fbs$\"1+++Xlm<>F17$Fgs$\"1+++T1)o!>F17$F\\t$\"1+++K+F'*=F17 $Fat$\"1+++\"3Ce)=F17$Fft$\"1+++(3Lb(=F17$F[u$\"1+++*>)Ql=F17$F`u$\"1+ ++&G\"Qb=F17$Feu$\"1+++i[]X=F17$Fju$\"1+++I?vN=F17$F_v$\"1+++.k6E=F17$ Fdv$\"1+++Q?f;=F17$Fiv$\"1+++8M<2=F17$F^w$\"1+++i`&yz\"F17$Fcw$\"1+++o Ij)y\"F17$Fhw$\"1+++(*>]z,%o\"F17$Fd[l$\"1+++A'*pv;F17$Fi[l $\"1+++5qVn;F17$F^\\l$\"1+++(z@#f;F17$Fc\\l$\"1+++&o^5l\"F17$Fh\\l$\"1 +++WW#Hk\"F17$F]]l$\"1+++sy$[j\"F17$Fb]l$\"1+++@)*yE;F17$Fg]l$\"1+++(> y(=;F17$F\\^l$\"1+++@4!3h\"F17$Fa^l$\"1+++9f&Gg\"F17$Ff^l$\"1+++H6%\\f \"F17$F[_l$\"1+++_X0(e\"F17$F`_l$\"1+++;T>z:F17$Fe_l$\"1+++kxNr:F17$Fj _l$\"1+++UMaj:F17$F_`l$\"1+++[!\\db\"F17$Fd`l$\"1+++gC(za\"F17$Fi`l$\" 1+++y9@S:F17$F^al$\"1+++pQYK:F17$Fcal$\"1+++?tsC:F17$Fhal$\"1+++X%**p^ \"F17$F]bl$\"1+++,xF4:F17$Fbbl$\"1+++([f:]\"F17$Fgbl$\"1+++Z?%Q\\\"F17 $F\\cl$\"1+++HC7'[\"F17$Facl$\"1+++%\\(Ry9F17$Ffcl$\"1+++qQmq9F17$F[dl $\"1+++(*y\"HY\"F17$F`dl$\"1+++Wc:b9F17$Fedl$\"1+++9FPZ9F17$Fjdl$\"1++ +'Hk&R9F17$F_el$\"1+++x]sJ9F17$Fdel$\"1+++V!\\QU\"F17$Fiel$\"1+++=$HfT \"F17$F^fl$\"1+++>\"ezS\"F17$Fcfl$\"1+++nj#**R\"F17$Fhfl$\"1+++rN#=R\" F17$F]gl$\"1+++Cqj$Q\"F17$Fbgl$\"1+++V;Nv8F17$Fggl$\"1+++X)[pO\"F17$F \\hl$\"1+++aaSe8F17$F`hl$\"1+++S;p\\8F17$Fdhl$\"1+++\"4o2M\"F17$Fihl$ \"1+++b4eJ8F17$F]il$\"1+++'>`?K\"F17$Fail$\"1+++Mj178F17$Feil$\"1++++B U,8F17$Fiil$\"1++++jt*G\"F17$F]jl$\"1+++oe&eE\"F17$F*$\"1+++iIx]7F17$F djl$\"1+++=$4PD\"F17$Fijl$\"1+++jn$pD\"F17$F^[m$\"1+++@&\\)f7F17$Fc[m$ \"1+++!HMDE\"F17$Fh[m$\"1+++gR/l7F17$F]\\m$\"1+++5ITn7F17$Fb\\m$\"1+++ 2fmp7F17$Fg\\m$\"1+++d-#=F\"F17$F\\]m$\"1+++5**)QF\"F17$Fa]m$\"1+++?c) eF\"F17$Fe]m$\"1+++tg\"yF\"F17$Fj]m$\"1++++\")oz7F17$F_^m$\"1+++0x]\"G \"F17$Fd^m$\"1+++<)zKG\"F17$Fi^m$\"1+++E(3]G\"F17$F^_m$\"1+++Nyp'G\"F1 7$Fc_m$\"1+++z.N)G\"F17$Fh_m$\"1+++&4p**G\"F17$F]`m$\"1+++Ekb\"H\"F17$ Fb`m$\"1+++LW6$H\"F17$Ff`m$\"1+++G]k%H\"F17$F[am$\"1+++x*\\hH\"F17$F`a m$\"1+++-3j(H\"F17$Feam$\"1+++T))3*H\"F17$Fjam$\"1+++w`_+8F17$F_bm$\"1 +++c:%>I\"F17$Fdbm$\"1+++G%QLI\"F17$Fibm$\"1+++Ipr/8F17$F^cm$\"1+++oz2 18F17$Fccm$\"1+++aBU28F17$Fgcm$\"1+++d3v38F17$F\\dm$\"1+++ST158F17$Fad m$\"1+++nGO68F17$Ffdm$\"1+++Swk78F17$F[em$\"1+++O!>RJ\"F17$F`em$\"1+++ Pv<:8F17$Feem$\"1+++lOU;8F17$Fjem$\"1+++*)yl<8F17$F_fm$\"1+++K1))=8F17 $Fdfm$\"1+++*H#4?8F17$Fhfm$\"1+++uKH@8F17$F]gm$\"1+++4R[A8F17$Fbgm$\"1 +++lXmB8F17$Fggm$\"1+++Xb$[K\"F17$F\\hm$\"1+++dr*fK\"F17$Fahm$\"1+++4( \\rK\"F17$Ffhm$\"1+++vMHG8F17$F[im$\"1+++*pG%H8F17$F`im$\"1+++TcbI8F17 $Feim$\"1+++EXnJ8F17$Fiim$\"1+++1cyK8F17$F^jm$\"1+++l!*)QL\"F17$Fcjm$ \"1+++A^)\\L\"F17$Fhjm$\"1+++xR2O8F17$F][n$\"1+++@e:P8F17$Fb[n$\"1+++I 3BQ8F17$Fg[n$\"1+++t\"*HR8F17$F\\\\n$\"1+++65OS8F17$Fa\\n$\"1+++9lTT8F 17$Ff\\n$\"1+++:eYU8F17$Fj\\n$\"1+++t!4NM\"F17$F_]n$\"1+++OkaW8F17$Fd] n$\"1+++H!ybM\"F17$Fi]n$\"1+++&)RgY8F17$F^^n$\"1+++AWiZ8F17$Fc^n$\"1++ +t%R'[8F17$Fh^n$\"1+++Z#\\'\\8F17$F]_n$\"1+++_Ql]8F17$Fb_n$\"1++++Ml^8 F17$Fg_n$\"1+++)*zk_8F17$F\\`n$\"1+++Zxj`8F17$Fa`n$\"1+++SFia8F17$Ff`n $\"1+++sIgb8F17$F[an$\"1+++Q)ylN\"F17$F`an$\"1+++E,bd8F17$Fean$\"1+++8 q^e8F17$Fjan$\"1+++(ez%f8F17$F_bn$\"1+++JzVg8F17$Fdbn$\"1+++G@Rh8F17$F ibn$\"1+++SAMi8F17$F]cn$\"1+++Y$)Gj8F17$Fbcn$\"1+++>0Bk8F17$Fgcn$\"1++ +H)o^O\"F17$F\\dn$\"1+++QL5m8F17$Fadn$\"1+++6T.n8F17$Ffdn$\"1+++=7'zO \"F17$F[en$\"1+++>Z))o8F17$F`en$\"1+++pY!)p8F17$Feen$\"1+++H6sq8F17$Fj en$\"1+++bTjr8F17$F^fn$\"1+++-Qas8F17$Fcfn$\"1+++@,Xt8F17$Fhfn$\"1+++p JNu8F17$F]gn$\"1+++&*HDv8F17$Fbgn$\"1+++^'\\hP\"F17$Fggn$\"1+++#=VqP\" F17$F\\hn$\"1+++PO$zP\"F17$Fahn$\"1+++i5#)y8F17$Ffhn$\"1+++/bqz8F17$F[ in$\"1++++qe!Q\"F1F_in-F$6$7ewFcbrFfbrFibrF\\crF_crFbcrFecrFhcrF[drF^d rFadrFddrFgdrFjdrF]erF`erFcerFferFierF\\frF_frFbfrFefrFhfrF[grF^grFagr FdgrFggrFjgrF]hrF`hrFchrFfhrFihrF\\irF_irFbirFeirFhirF[jrF^jrFajrFdjrF gjrFjjrF][sF`[sFc[sFf[sFi[sF\\\\sF_\\sFb\\sFe\\sFh\\sF[]sF^]sFa]sFd]sF g]sFj]sF]^sF`^sFc^sFf^sFi^sF\\_sF__sFb_sFe_sFh_sF[`sF^`sFa`sFd`sFg`sFj `sF]asF`asFcasFfasFiasF\\bsF_bsFbbsFebsFhbsF[csF^csFacsFdcsFgcsFjcsF]d sF`dsFcdsFfdsFids7$F]jl$\"1+++*3liG\"F1F_esFbesFeesFhesF[fsF^fsFafsFdf sFgfsFjfsF]gsF`gsFcgsFfgsFigsF\\hsF_hsFbhsFehsFhhsF[isF^isFaisFdisFgis FjisF]jsF`jsFcjsFfjsFijsF\\[tF_[tFb[tFe[tFh[tF[\\tF^\\tFa\\tFd\\tFg\\t Fj\\tF]]tF`]tFc]tFf]tFi]tF\\^tF_^tFb^tFe^tFh^tF[_tF^_tFa_tFd_tFg_tFj_t F]`tF``tFc`tFf`tFi`tF\\atF_atFbatFeatFhatF[btF^btFabtFdbtFgbtFjbtF]ctF `ctFcctFfctFictF\\dtF_dtFbdtFedtFhdtF[etF^etFaetFdetFgetFjetF]ftF`ftFc ftFfftFiftF\\gtF_gtFbgtFegtFhgtF[htF_in-F$6$7ew7$F($\"1+++e6O-DF17$F/$ \"1+++WXfDGF17$F5$\"1+++6q%z#HF17$F:$\"1+++>KA/IF17$F?$\"1+++KwhnIF17$ FD$\"1+++*Q8I7$F17$FI$\"1+++%RVG<$F17$FN$\"1+++\"fA&=KF17$FS$\"1+++-v& 4E$F17$FX$\"1+++k$p2I$F17$Fgn$\"1+++o_SQLF17$F\\o$\"1+++'Q*>uLF17$Fao$ \"1+++L!4%3MF17$Ffo$\"1+++5tBTMF17$F[p$\"1+++^x%GZ$F17$F`p$\"1++++VP.N F17$Fep$\"1+++*>GH`$F17$Fjp$\"1+++5JghNF17$F_q$\"1+++7(y%*e$F17$Fdq$\" 1+++]Mi;OF17$Fiq$\"1+++ol4VOF17$F^r$\"1+++Y)\\*oOF17$Fcr$\"1+++v(GUp$F 17$Fhr$\"1+++DO(*=PF17$F]s$\"1+++[,AVPF17$Fbs$\"1+++s.+nPF17$Fgs$\"1++ +\\IM!z$F17$F\\t$\"1+++OTF8QF17$Fat$\"1+++0r\"e$QF17$Fft$\"1+++GL*z&QF 17$F[u$\"1+++/B#)zQF17$F`u$\"1+++#)=K,RF17$Feu$\"1+++8%3D#RF17$Fju$\"1 +++()pRVRF17$F_v$\"1+++C:+kRF17$Fdv$\"1+++.\\L%)RF17$Fiv$\"1+++W!4W+%F 17$F^w$\"1+++\\]BCSF17$Fcw$\"1+++cK#Q/%F17$Fhw$\"1+++8L=jSF17$F]x$\"1+ ++9UK#3%F17$Fbx$\"1+++4WD,TF17$Fgx$\"1+++;=)*>TF17$F\\y$\"1+++#)Q^QTF1 7$Fay$\"1+++&edo:%F17$Ffy$\"1+++E&>]<%F17$F[z$\"1+++Jf+$>%F17$F`z$\"1+ ++&oA3@%F17$Fez$\"1+++C`ZGUF17$Fjz$\"1+++F\"pfC%F17$F_[l$\"1+++'34LE%F 17$Fd[l$\"1+++\\**\\!G%F17$Fi[l$\"1+++&>YvH%F17$F^\\l$\"1+++>@X9VF17$F c\\l$\"1+++(z@7L%F17$Fh\\l$\"1+++A\"fyM%F17$F]]l$\"1+++#ynVO%F17$Fb]l$ \"1+++@8v!Q%F17$Fg]l$\"1+++BJ,(R%F17$F\\^l$\"1+++Ck:8WF17$Fa^l$\"1+++, V=HWF17$Ff^l$\"1+++A(*4XWF17$F[_l$\"1+++Fb!4Y%F17$F`_l$\"1+++VWgwWF17$ Fe_l$\"1+++o!*>#\\%F17$Fj_l$\"1+++<>p2XF17$F_`l$\"1+++-a3BXF17$Fd`l$\" 1+++Y=QQXF17$Fi`l$\"1+++kMe`XF17$F^al$\"1+++6CpoXF17$Fcal$\"1+++`2r$e% F17$Fhal$\"1+++*[S')f%F17$F]bl$\"1+++EN[8YF17$Fbbl$\"1+++Ct7#\\F17$F]il$\"1+++u'HX$\\F17$Fail$\"1+++\\3tZ\\F17$Feil$\"1+ ++;w(3'\\F17$Fiil$\"1+++R3(R(\\F17$F]jl$\"1+++k8,()\\F17$F*Feim7$Fdjl$ \"1+++]v$H,&F17$Fijl$\"1+++&zCe-&F17$F^[m$\"1+++'\\i'Q]F17$Fc[m$\"1+++ &Q^90&F17$Fh[m$\"1+++!>#>k]F17$F]\\m$\"1+++;c)o2&F17$Fb\\m$\"1+++bB`*3 &F17$Fg\\m$\"1+++rI8-^F17$F\\]m$\"1+++A%)o9^F17$Fa]m$\"1+++^!*>F^F17$F e]m$\"1+++*el'R^F17$Fj]m$\"1+++O')3_^F17$F_^m$\"1+++(zoW;&F17$Fd^m$\"1 +++am!o<&F17$Fi^m$\"1*****py-\"*=&F17$F^_m$\"1+++TxN,_F17$Fc_m$\"1+++q ?d8_F17$Fh_m$\"1+++3juD_F17$F]`m$\"1*****H)4)yB&F17$Fb`m$\"1+++(fw*\\_ F17$Ff`m$\"1+++cO.i_F17$F[am$\"1+++_E0u_F17$F`am$\"1+++qS.'G&F17$Feam$ \"1*****4PyzH&F17$Fjam$\"1+++?g))4`F17$F_bm$\"1+++quv@`F17$Fdbm$\"1+++ mJfL`F17$Fibm$\"1+++MNRX`F17$F^cm$\"1*****R+frN&F17$Fccm$\"1+++$***))o `F17$Fgcm$\"1*****H\"pe!Q&F17$F\\dm$\"1+++d,D#R&F17$Fadm$\"1+++@,)QS&F 17$Ffdm$\"1+++#>xaT&F17$F[em$\"1+++_ Xof&F17$Fiim$\"1+++uc$zg&F17$F^jm$\"1+++2\"***=cF17$Fcjm$\"1+++xd.IcF1 7$Fhjm$\"1+++_f/TcF17$F][n$\"1+++-*H?l&F17$Fb[n$\"1+++\")y)Hm&F17$Fg[n $\"1+++\\,#Rn&F17$F\\\\n$\"1+++ep#[o&F17$Fa\\n$\"1+++h&3dp&F17$Ff\\n$ \"1+++%>llq&F17$Fj\\n$\"1+++-rRrdF17$ Fh^n$\"1+++dB)=y&F17$F]_n$\"1+++Cwa#z&F17$Fb_n$\"1*****H'**=.eF17$Fg_n $\"1+++#e4Q\"eF17$F\\`n$\"1+++#p1W#eF17$Fa`n$\"1*****z[\")\\$eF17$Ff`n $\"1+++tT`XeF17$F[an$\"1+++Z\\1ceF17$F`an$\"1+++1SdmeF17$Fean$\"1+++M: 1xeF17$Fjan$\"1+++Dx_()eF17$F_bn$\"1+++lF(z*eF17$Fdbn$\"1+++SoR3fF17$F ibn$\"1+++C,!)=fF17$F]cn$\"1+++*z#=HfF17$Fbcn$\"1*****4/X&RfF17$Fgcn$ \"1+++Fq))\\fF17$F\\dn$\"1+++@*3-'fF17$Fadn$\"1+++%*3^qfF17$Ffdn$\"1++ +9Jz!)fF17$F[en$\"1+++Zd0\"*fF17$F`en$\"1+++\\*)H,gF17$Feen$\"1+++$)G_ 6gF17$Fjen$\"1+++1xs@gF17$F^fn$\"1+++xN\">.'F17$Fcfn$\"1+++V13UgF17$Fh fn$\"1+++d!HA0'F17$F]gn$\"1+++r*eB1'F17$Fbgn$\"1+++L0ZsgF17$Fggn$\"1++ +%)Qc#3'F17$F\\hn$\"1+++q\"RE4'F17$Fahn$\"1+++Llp-hF17$Ffhn$\"1+++;ht7 hF17$F[in$\"1+++]!eF7'F1F_in-%*AXESTICKSG6%%(DEFAULTGFc^w-%%FONTG6%%&T IMESG%'ITALICG\"#5F_in-%+AXESLABELSG6$%!GF]_w" 1 2 0 1 0 2 6 1 4 2 1.000000 45.000000 45.000000 0 }}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 11 "Eigenvalues" }}{PARA 0 "" 0 "" {TEXT -1 214 "The roots of the char acteristic polynomial of Kobbelt's scheme in general cannot be found e xplicitly. However, we can obtain enough information about the eigenva lues to verify C1-continuity. We proove that for any " }{XPPEDIT 18 0 "m, k" "6$%\"mG%\"kG" }{TEXT -1 3 ", " }{XPPEDIT 18 0 "m = 1..k-1" "6 #/%\"mG;\"\"\",&%\"kG\"\"\"\"\"\"!\"\"" }{TEXT -1 60 " the largest ei genvalue is real and unique, and that for " }{XPPEDIT 18 0 " m <> k- 1,1" "6$0%\"mG,&%\"kG\"\"\"\"\"\"!\"\"\"\"\"" }{TEXT -1 76 " the large st eigenvalue is less than the largest eigenvalue of blocks 1 and " } {XPPEDIT 18 0 "k-1" "6#,&%\"kG\"\"\"\"\"\"!\"\"" }{TEXT -1 109 ". We \+ also show that the unique largest eigenvalue is a single eigenvalue in the interval [0.5..0.613], for " }{XPPEDIT 18 0 "k" "6#%\"kG" } {TEXT -1 22 " > 4. \nFor the value " }{XPPEDIT 18 0 "k = 3" "6#/%\"kG \"\"$" }{TEXT -1 157 ", eigenvalues are examined separately. The proo f is performed in several steps: \n(1). We show that for c < 0, all \+ roots of the characteristic polynomial " }{XPPEDIT 18 0 "P(c,lambda) " "6#-%\"PG6$%\"cG%'lambdaG" }{TEXT -1 169 " are less than 0.51 (actua lly, they are less than 0.5, but due to numerical nature of our calcul ations, we have to relax the upper boundary). \n(2). We show that for \+ any " }{XPPEDIT 18 0 "c = 0..1" "6#/%\"cG;\"\"!\"\"\"" }{TEXT -1 30 ", there is a unique real root " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 47 " in the interval [0.5,0.613], and the function " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 117 " is C1-continuous and incre ases. \n(3). We \"deflate\" the characteristic polynomial (that is, \+ divide by the monomial " }{XPPEDIT 18 0 "lambda - mu" "6#,&%'lambdaG\" \"\"%#muG!\"\"" }{TEXT -1 26 ") in symbolic form, with " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "c" "6#%\"cG" } {TEXT -1 50 " as indeterminates. Next, we verify that for all " } {XPPEDIT 18 0 "mu = 0.5..0.613" "6#/%#muG;$\"\"&!\"\"$\"$8'!\"$" } {TEXT -1 20 ", and corresponding " }{XPPEDIT 18 0 "c(mu)" "6#-%\"cG6#% #muG" }{TEXT -1 152 " , that all roots of the deflated polynomial are \+ inside the circle of radius 0.5 centered at 0 in the complex plane, th at is, have magnitudes less than " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#% \"cG" }{TEXT -1 10 " for any " }{XPPEDIT 18 0 "c > 0" "6#2\"\"!%\"cG " }{TEXT -1 8 ". Using " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 77 " \+ as the primary parameter is important, as c can be explicitly computed from " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 32 ", but not the othe r way.\nAs for " }{XPPEDIT 18 0 " k" "6#%\"kG" }{TEXT -1 6 " > 4, " } {XPPEDIT 18 0 "cos(2*Pi/k) > 0.51" "6#2$\"#^!\"#-%$cosG6#*(\"\"#\"\"\" %#PiGF,%\"kG!\"\"" }{TEXT -1 63 ", the largest eigenvalue cannot possi bly correspond to a block " }{XPPEDIT 18 0 "m" "6#%\"mG" }{TEXT -1 12 ", for which " }{XPPEDIT 18 0 "cos(2*m*Pi/k) <= 0\n" "6#1-%$cosG6#**\" \"#\"\"\"%\"mGF)%#PiGF)%\"kG!\"\"\"\"!" }{TEXT -1 72 " . From (3), i t follows that the largest root has to be the real root " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 10 " for some " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 14 ". As for any " }{XPPEDIT 18 0 "m > 1, m < k-1" "6$2\"\"\"%\"mG2F%,&%\"kG\"\"\"\"\"\"!\"\"" }{TEXT -1 2 ", " } {XPPEDIT 18 0 "cos(2*m*Pi/k) < cos(2*Pi/k)" "6#2-%$cosG6#**\"\"#\"\"\" %\"mGF)%#PiGF)%\"kG!\"\"-F%6#*(\"\"#F)F+F)F,F-" }{TEXT -1 30 " , and w e have shown (1) that " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" } {TEXT -1 25 ", increases, and for any " }{XPPEDIT 18 0 "c" "6#%\"cG" } {TEXT -1 1 " " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 370 " is the largest root, we conclude that the largest eigenvalue always \+ corresponds to m = 1, is real, and is the unique eigenvalue in the ran ge 0.5..0.613.\n\n On steps 1 and 3 we have to show that roots of a po lynomial are inside a circle of radius r in the complex plane. This \+ task is similar to the task of establishing\nstability of a filter wit h the transfer function " }{XPPEDIT 18 0 "1/a(z)" "6#*&\"\"\"\"\"\"-% \"aG6#%\"zG!\"\"" }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "a(z)" "6#-%\" aG6#%\"zG" }{TEXT -1 706 " is a polynomial. Such filter is stable, if \+ all roots of the polynomial are inside the unit circle. \nA variety o f tests exist for this condition; for our purposes, the algebraic Mard en-Jury test is convenient. With aproporiate rescaling of the variable it can be used to prove that all roots of a polynomial are inside the circle of any given radius r. As the test requires only a simple alge braic calculaion on the coefficients of the polynomial, it can be easi ly performed for symbolic and interval coefficients.\nFinally, we comp ute the largest root of the characteristic polynomial numerically for \+ all valences up to some maximum. For each computed root, we verify th at that the precision is at least " }{XPPEDIT 18 0 "epsilon = 1e-11" " 6#/%(epsilonG$\"\"\"!#6" }{TEXT -1 61 ": we use interval arithmetics t o evaluate the polynomial at " }{XPPEDIT 18 0 "lambda[0] - epsilon" " 6#,&&%'lambdaG6#\"\"!\"\"\"%(epsilonG!\"\"" }{TEXT -1 7 " and " } {XPPEDIT 18 0 "lambda[0]+epsilon" "6#,&&%'lambdaG6#\"\"!\"\"\"%(epsilo nGF(" }{TEXT -1 270 " and assert that the sign is guaranteed to change . There may be more than one root: we still have to prove that ther e is only a single root in the computed interval and that the rest of \+ the roots are smaller. The maximal valence N is chosen in such a way \+ that for N " }{XPPEDIT 18 0 "cos(2*Pi/N)" "6#-%$cosG6#*(\"\"#\"\"\"%#P iGF(%\"NG!\"\"" }{TEXT -1 121 "is \"sufficiently close\" to 1. This \+ means that for all K > N corresponding eigenvalue differs from th e limit value " }{XPPEDIT 18 0 "lambda[infinity]" "6#&%'lambdaG6#%)inf inityG" }{TEXT -1 17 " by no more than " }{XPPEDIT 18 0 "epsilon" "6#% (epsilonG" }{TEXT -1 10 " , where " }{XPPEDIT 18 0 "epsilon" "6#%(eps ilonG" }{TEXT -1 190 " is small enough for us to establish, using in terval arithmetics, that the Jacobian of the characteristic map is po sitive for eigenvectors computed using formulas derived below for all " }{XPPEDIT 18 0 "lambda " "6#%'lambdaG" }{TEXT -1 17 " in the inter val " }{XPPEDIT 18 0 "[lambda[infinity]-epsilon,lambda[infinity]]" "6# 7$,&&%'lambdaG6#%)infinityG\"\"\"%(epsilonG!\"\"&F&6#F(" }{TEXT -1 249 " . The actual computation of the Jacobian and evaluation of the n ecessary contraction functions is performed in the C part of the code. We use maximal value 3000 here, just in case ( below we see that 145 0 is sufficient to require the interval for " }{XPPEDIT 18 0 "lambda" "6#%'lambdaG" }{TEXT -1 24 " with the size of only " }{XPPEDIT 18 0 " 1e-6" "6#$\"\"\"!\"'" }{TEXT -1 2 ".\n" }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 16 "Marden-Jury test" }}{PARA 0 "" 0 "" {MPLTEXT 0 21 26 "Mar denJury(a,var,rootrad) " }{TEXT -1 44 "Compute Marden-Jury table for \+ a polynomial " }{MPLTEXT 0 21 1 "p" }{TEXT -1 13 " in variable " } {MPLTEXT 0 21 3 "var" }{TEXT -1 117 ", with variable rescaled by rootr ad.\nUsed to verify that all roots of a polynomial are inside the circ le of radius r." }}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "MardenJury \+ := proc( a::polynom, var::name, rootrad) local i,k,M, acol, tbl, resta ble;\n M := degree(a,var);\n for i from 0 to M do \n tbl[0,i] \+ := coeff(a, var,i)*rootrad^(i-M);\n od;\n for i from 1 to M do\n \+ for k from 0 to M-i do\n tbl[i,k] := tbl[i-1,0]*tbl[i-1,k]- t bl[i-1,M-k-i+1]*tbl[i-1,M-i+1];\n od\n od;\n for i from 1 to M \+ do \n restable[i] := tbl[i,0]; \n od;\n eval(restable);\nend: " "6#>%+MardenJuryGR6%'%\"aG%(polynomG'%$varG%%nameG%(rootradG 7(%\"iG%\"kG%\"MG%%acolG%$tblG%)restableG6\"F5C'>F1-%'degreeG6$F(F+?(F /\"\"!\"\"\"F1%%trueG>&F36$F?(F0F>&F36$F/F0,&*&&F36$,&F/F=\"\"\"FHF >&F46#F/&F 36$F/F<-%%evalG6#F4F5F5F5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 36 "Inte rval version of Marden-Jury test" }}{PARA 0 "> " 0 "" {XPPEDIT 19 1 "I ntervMardenJury := proc( a::polynom(interval), var::name, rootrad::num eric) local i,k,M, acol, tbl, restable;\n # checking correct type of a is too messy\n M := degree(a,var);\n for i from 0 to M do \n \+ tbl[0,i] := Interval_times(coeff(a, var,i), rootrad^(i-M));\n od; \n for i from 1 to M do\n for k from 0 to M-i do\n tbl[i,k] := Interval_add( Interval_times( tbl[i-1,0], tbl[i-1,k]), \n \+ Interval_times(-1, Interval_times( tbl[i-1,M-k-i+1], tbl[i-1,M-i+1]) ));\n od\n od;\n for i from 1 to M do \n restable[i] := t bl[i,0]; \n od;\n eval(restable);\nend: " "6#>%1IntervMarde nJuryGR6%'%\"aG-%(polynomG6#%)intervalG'%$varG%%nameG'%(rootradG%(nume ricG7(%\"iG%\"kG%\"MG%%acolG%$tblG%)restableG6\"F:C'>F6-%'degreeG6$F(F .?(F4\"\"!\"\"\"F6%%trueG>&F86$FAF4-%/Interval_timesG6$-%&coeffG6%F(F. F4)F1,&F4FBF6!\"\"?(F4\"\"\"FBF6FC?(F5FAFB,&F6FBF4FOFC>&F86$F4F5-%-Int erval_addG6$-FH6$&F86$,&F4FB\"\"\"FOFA&F86$,&F4FB\"\"\"FOF5-FH6$,$\"\" \"FO-FH6$&F86$,&F4FB\"\"\"FO,*F6FBF5FOF4FO\"\"\"FB&F86$,&F4FB\"\"\"FO, (F6FBF4FO\"\"\"FB?(F4\"\"\"FBF6FC>&F96#F4&F86$F4FA-%%evalG6#F9F:F:F:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 9 "Deflation" }}{PARA 0 "" 0 "" {MPLTEXT 0 21 23 "deflate(p,v ar,rootval) " }{TEXT -1 43 "compute the coefficients of the polynomial " }{XPPEDIT 18 0 "p(z)/(z -z[0])" "6#*&-%\"pG6#%\"zG\"\"\",&F'F(&F'6# \"\"!!\"\"F-" }{TEXT -1 39 "; it is assumed that p is divisible by " } {XPPEDIT 18 0 "z-z[0]" "6#,&%\"zG\"\"\"&F$6#\"\"!!\"\"" }{TEXT -1 3 ". " }{MPLTEXT 0 21 3 "var" }{TEXT -1 30 " is the name of the variable, " }{MPLTEXT 0 21 7 "rootval" }{TEXT -1 13 " is the root." }}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "deflate := proc(p::polynom, var::name , rootval)\nlocal i, dp,r; \n dp := 0;\n r := lcoeff(p,var);\n \+ for i from degree(p,var)-1 to 0 by -1 do \n dp := dp + r*var^i ;\n r := coeff(p, var, i ) + rootval*r;\n od; \n dp;\nend: \+ \n" "6#>%(deflateGR6%'%\"pG%(polynomG'%$varG%%nameG%(rootvalG7%%\"iG %#dpG%\"rG6\"F2C&>F0\"\"!>F1-%'lcoeffG6$F(F+?(F/,&-%'degreeG6$F(F+\"\" \"\"\"\"!\"\",$\"\"\"FAF5%%trueGC$>F0,&F0F?*&F1F?)F+F/F?F?>F1,&-%&coef fG6%F(F+F/F?*&F-F?F1F?F?F0F2F2F2" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 27 "Analysis of the eigenvalues" }}{PARA 0 "" 0 "" {TEXT -1 41 "Now we perform steps 1-3 described above." }}{SECT 0 {PARA 0 "" 0 "" {TEXT -1 73 "(1). We show that for c < 0, all roots of the characteri stic polynomial " }{XPPEDIT 18 0 "P(c,lambda)" "6#-%\"PG6$%\"cG%'lambd aG" }{TEXT -1 142 " are less than 0.51 (actually, they are less than 0 .5, but due to numerical nature of our calculations, we have to relax \+ the upper boundary). " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "MJt ab := MardenJury(KobFactor6, lambda, 51/100): " }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 66 "MJtabInterv := map(unapply('inapply'( dummy, c ),dummy), MJtab ): " }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "TestN egativeC := proc( cstart::numeric, cend::numeric,cstep::numeric ) \n \+ local MJ,cx; global MJtabInterv;\n for cx from cstart to cend by cst ep do \n MJ := map( unapply( dummy([cx,cx+cstep]),dummy), MJtabInte rv);\n if op(2, MJ[1]) > 0 or op(1, MJ[2]) < 0 or op(1,MJ[3]) < 0 o r op(1, MJ[4]) < 0 \n or op(1,MJ[5]) < 0 or op(1, MJ[6]) < 0 the n ERROR(`test failed for interval = `, [cx,cx+cstep]); \n fi;\n od ;\n print(`All tests passed`);\nend:" "6#>%.TestNegativeCGR6%'%'cstart G%(numericG'%%cendGF)'%&cstepGF)7$%#MJG%#cxG6\"F1C$?(F0F(F-F+%%trueGC$ >F/-%$mapG6$-%(unapplyG6$-%&dummyG6#7$F0,&F0\"\"\"F-FBF>%,MJtabIntervG @$555552\"\"!-%#opG6$\"\"#&F/6#\"\"\"2-FM6$\"\"\"&F/6#\"\"#FK2-FM6$\" \"\"&F/6#\"\"$FK2-FM6$\"\"\"&F/6#\"\"%FK2-FM6$\"\"\"&F/6#\"\"&FK2-FM6$ \"\"\"&F/6#\"\"'FK-%&ERRORG6$% " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 148 "Do tests, ad justing the step in c. This is not really necessary -- we can simply t ake the smallest step, but to save time we use larger steps first." } {MPLTEXT 1 0 32 "\nTestNegativeC(-1.0,-0.65,0.05);" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#%1All~tests~passedG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "TestNegativeC(-0.6,-0.275,0.025);" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#%1All~tests~passedG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "TestNegativeC(-0.25,-0.06,0.01);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%1All~tests~passedG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "TestNegativeC(-0.05,0.,0.005);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%1All~tests~passedG" }}}}{SECT 0 {PARA 0 "" 0 "" {TEXT -1 26 "(2). We show that for any " }{XPPEDIT 18 0 "c = 0..1" "6#/%\"cG ;\"\"!\"\"\"" }{TEXT -1 30 ", there is a unique real root " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 47 " in the interval [0.5,0.613], and th e function " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 32 " i s C1-continuous and increases." }}{PARA 5 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 39 "Solve the characteritic polynoial \+ for c" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "csolutions := [solve( KobF actor6, c)]:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 136 "We are intersted in the first solution only; we will verify later that te second one i s out of the range [-1..1] for relevant values of " }{XPPEDIT 18 0 "la mbda" "6#%'lambdaG" }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "clambda := csolutions[1];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(clambdaG,$*(%'lambdaG!\"#,,*$F'\"\"$!%Wh*$F'\"\"#!%? >F'\"$K%!#=\"\"\"*$,0*$F'\"\"'\"(%=P%**$F'\"\"&\")s#QK\"*$F'\"\"%!(w " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 24 "Compute the derivative. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "clambdadiff := simplify( diff(clambda, lambda)):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "The solution for " }{XPPEDIT 18 0 "lambda = 1/2" "6#/%'lambdaG*&\"\"\"\"\"\"\"\"#!\"\"" }{TEXT -1 28 " can be c omputed explicitly." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "simplify( su bs( lambda = 1/2, clambda)); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"! " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "The solution for " }{XPPEDIT 18 0 "lambda = 0.613" "6#/%'lambdaG$\"$8'!\"$" }{TEXT -1 22 " is outsi de the range." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "clambdaInterv := i napply( clambda, lambda): clambdaInterv(0.613);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$$\"+ToB25!\"*$\"+jrB25F&" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 41 "Show that t he derivative is positive for " }{XPPEDIT 18 0 "lambda" "6#%'lambdaG" }{TEXT -1 61 " in [0.5..0.613] (the upper bound is the upper estimate for " }{XPPEDIT 18 0 "lambda(1)" "6#-%'lambdaG6#\"\"\"" }{MPLTEXT 1 0 47 " \nintervcdiff := inapply( clambdadiff, lambda):" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "for xl from 0.5 by 0.004 to 0.613 do \+ \n res := intervcdiff([xl,xl+0.004]):\n if res[1] < 0 then ERROR(` test failed for interval`, [xl,xl+0.004] );\n fi;\nod:\nprint(`all t ests passed`);" "6#C$?(%#xlG$\"\"&!\"\"$\"\"%!\"$$\"$8'!\"$%%trueGC$>% $resG-%,intervcdiffG6#7$F%,&F%\"\"\"$\"\"%!\"$F8@$2&F26#\"\"\"\"\"!-%& ERRORG6$%9test~failed~for~intervalG7$F%,&F%F8$\"\"%!\"$F8-%&printG6#%1 all~tests~passedG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%1all~tests~passe dG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "We conclude that " } {XPPEDIT 18 0 "c[1](lambda)" "6#-&%\"cG6#\"\"\"6#%'lambdaG" }{TEXT -1 115 " increases from 0 to above 1 on [0.5..0.613]; therefore, the inve rse increases from 0.5 to approx. 0.613 on [0..1]." }{MPLTEXT 1 0 1 " \+ " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 64 "The second solution is outsid e the range of c for this range of " }{XPPEDIT 18 0 "lambda" "6#%'lamb daG" }{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "inapply( cso lutions[2], lambda)(0.5, 0.613);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$ $!+F+]7G!\")$!+w**\\7GF&" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 86 "We co nclude that for c > 0 in the interval 0.5..0.613 there is a unique rea l solution." }}}}{SECT 0 {PARA 0 "" 0 "" {TEXT -1 83 "(3). We \"deflat e\" the characteristic polynomial (that is, divide by the monomial \+ " }{XPPEDIT 18 0 "lambda - mu" "6#,&%'lambdaG\"\"\"%#muG!\"\"" }{TEXT -1 30 ") in the symbolic form, with " }{XPPEDIT 18 0 "mu" "6#%#muG" } {TEXT -1 5 " and " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 53 " as the \+ indeterminates. Next, we verify that for all " }{XPPEDIT 18 0 "c = 0.. 1" "6#/%\"cG;\"\"!\"\"\"" }{TEXT -1 14 ", and for all " }{XPPEDIT 18 0 "lambda = 0.5..0.613" "6#/%'lambdaG;$\"\"&!\"\"$\"$8'!\"$" }{TEXT -1 146 ", all roots of the deflated polynomial are inside the circle o f radius 0.5 centered at 0 in the complex plane, that is, have magnitu des less than " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 10 " for any " }{XPPEDIT 18 0 "c > 0" "6#2\"\"!%\"cG" }{TEXT -1 2 ". " } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "S ymbolic deflation; if we substitute a pair " }{XPPEDIT 18 0 "c, mu(c) " "6$%\"cG-%#muG6#F#" }{TEXT -1 58 " we get the deflated polynomial fo r a specific value of c." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "deflate dKobFactor6 := collect( expand( deflate(KobFactor6, lambda, mu)), lamb da);" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%3deflatedKobFactor6G,D*$%'la mbdaG\"\"&\"\"\"*&,(#!#:\"#;F)%\"cG#!\"$\"#k%#muGF)F)F'\"\"%F)*&,,F/#! \"*\"%C5#\"$(HF9F)F3F,*&F3F)F/F)F0*$F3\"\"#F)F)F'\"\"$F)*&,2F/#\"#@\"% '4%#!$N$\"%#>)F)*$F/F>#!\"(\"&%Q;FF/F)F0*$F3F?F)F)F'F >F)*&,6#\"$$=\"&Ob'F)F/#F8FKFFIFLF7F=F:FMF,*&F3F?F/F) F0*$F3F4F)F)F'F)F)F/#\"\"*\"')GC&#!#XFYF)F=FEFF/F>FI" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 34 "Verify that for all c in 0..1 and " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 64 " in 0.5.. 0.613 deflated polynomial has roots of magnitu de < 0.5" }}{PARA 0 "> " 0 "" {XPPEDIT 19 1 "TestDeflated := proc( ls tart::numeric, lend::numeric,lstep::numeric ) \n local cf,i, MJ,lx, c interv, deflatedinterv; global deflatedMJtabInterv,clambdaInterv; \n \+ for i from 0 to 5 do \n cf[i] := inapply( coeff(deflatedKobFactor 6,lambda,i), c, mu);\n od; \n for lx from lstart to lend by lstep do \n cinterv := clambdaInterv([lx,lx+lstep]);\n deflatedinterv := 0;\n for i from 0 to 4 do \n deflatedinterv := deflatedinterv + eval(cf[i](cinterv, [lx,lx+lstep]))*lambda^i; \n od;\n # f ix for cf[5]: stupid inapply does not always return an interval\n d eflatedinterv := deflatedinterv + [1.0,1.0]*lambda^5;\n MJ := Inter vMardenJury( deflatedinterv, lambda, 0.5);\n #MJ := map( unapply( d ummy(,[lx,lx+lstep]),dummy), deflatedMJtabInterv);\n if op(2, MJ[1] ) > 0 or op(1, MJ[2]) < 0 or op(1,MJ[3]) < 0 or op(1, MJ[4]) < 0 \n \+ or op(1,MJ[5]) < 0 or op(1, MJ[6]) < 0 then ERROR(`test failed for interval = `, [lx,lx+lstep]); fi;\n od;\n print(`Al l tests passed`);\nend:" "6#>%-TestDeflatedGR6%'%'lstartG%(numericG'%% lendGF)'%&lstepGF)7(%#cfG%\"iG%#MJG%#lxG%(cintervG%/deflatedintervG6\" F5C%?(F0\"\"!\"\"\"\"\"&%%trueG>&F/6#F0-%(inapplyG6%-%&coeffG6%%3defla tedKobFactor6G%'lambdaGF0%\"cG%#muG?(F2F(F-F+F;C(>F3-%.clambdaIntervG6 #7$F2,&F2F9F-F9>F4F8?(F0F8F9\"\"%F;>F4,&F4F9*&-%%evalG6#-&F/6#F06$F37$ F2,&F2F9F-F9F9)FFF0F9F9>F4,&F4F9*&7$$\"#5!\"\"$\"#5!\"\"F9*$FF\"\"&F9F 9>F1-%1IntervMardenJuryG6%F4FF$\"\"&!\"\"@$555552F8-%#opG6$\"\"#&F16# \"\"\"2-Ffp6$\"\"\"&F16#\"\"#F82-Ffp6$\"\"\"&F16#\"\"$F82-Ffp6$\"\"\"& F16#\"\"%F82-Ffp6$\"\"\"&F16#\"\"&F82-Ffp6$\"\"\"&F16#\"\"'F8-%&ERRORG 6$% " 0 "" {MPLTEXT 1 0 33 "TestDeflated(0.5, 0.613, 0.0005);" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#%1All~tests~passedG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 107 "We conclude that in the range c = 0..1, all roots of the defla ted polynomial have magnitudes less than 0.5 " }{MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 0 "" 0 " " {TEXT -1 19 "Special case: k = 3" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 99 "All we have to do is to run Marden-Jury test for the deflated char acteristic polynomial for k = 3" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "mu3 := fsolve( subs( c = cos(2*Pi/3), KobFactor6 ), lambda, 0.3..1 ); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$mu3G$\"+\"44LE%!#5" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 19 "Verify precision: " }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 65 "inapply( subs( c = cos(2*Pi/3), KobFactor6), l ambda)(mu3 + 1e-7);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$$\"+******HJ! #>$\"+,++!*QF&" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "inapply( \+ subs( c = cos(2*Pi/3), KobFactor6), lambda)(mu3 - 1e-7);" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#7$$!+,++IR!#>$!+******pJF&" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "intervmu3 := [mu3 - 1e-7, mu3 + 1e-7];" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%*intervmu3G7$$\"+\"43LE%!#5$\"+\"45L E%F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "Construct the interval deflated polynomial" } {MPLTEXT 1 0 299 " \n for i from 0 to 4 do \n cf3[i] := inapply( \+ coeff(subs( c = cos(2*Pi/3), deflate(KobFactor6,lambda,mu)),lambda,i) , mu);\n od: \n cf3[5] := mu -> [1.0,1.0]: \n deflatedinterv3 := 0 :\n for i from 0 to 5 do \n deflatedinterv3:= deflatedinterv3 + \+ eval(cf3[i](intervmu3))*lambda^i; \n od:\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 93 "Verify that all other roots have magnitude < 0.3: th e Marden Jury test for radius 0.3 passes:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "IntervMardenJury(deflatedinterv3, lambda, 0.3);" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#-%&TABLEG6#7'/\"\"\"7$$!+N;******!#5$! +N9******F,/\"\"#7$$\"+o,j$***F,$\"+l#RO***F,/\"\"$7$$\"+VAfj&*F,$\"+2 *pPc*F,/\"\"%7$$\"+\"fsF^&F,$\"+o/99bF,/\"\"&7$$\"+>=9g\\!#6$\"+T))\\( *\\FI" }}}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 67 "Calculation of the l argest eigenvalues with guaranteed precision " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 307 "This function produces a table of approximate values \+ of the eigenvalue with given precision for use with interval arithmeti cs in the C part of the analysis code; to avoid conversion problems, w e write two integers: mantissa + exponent base 10. The last value is \+ the limit value for infinity (computed with " }{XPPEDIT 18 0 "c" "6#% \"cG" }{TEXT -1 300 " set to 1 in the char. polynomial). The result i s a C function written to a file; if the file name is `default`, then \+ the output is written to the standard output.\nThe argument of the fun ction is valence, the function returns the interval value for the lar gest eigenvalue. The body is just a large " }{TEXT 259 6 "switch" } {TEXT -1 122 " statement.\nWe assume that the exponents for eigenvalu es are nonpositive, which is always the case for Kobbelt's scheme.\n" }{XPPEDIT 19 1 "\nComputeEigenvalues := proc(N::integer, eps::numeric, fname::string) \n local K, intervKobFactor6, intervPi, intervc, \n \+ expandedKobFactor6, approxEV, r, deflatedKobFactor6, i,marTable, cK; \+ \n global KobFactor6; \n Digits := 15;\n intervKobFactor6 := inapp ly( KobFactor6, lambda, c);\n # use arccos to make sure we get an int erval for Pi\n intervPi := Interval_times([2.0,2.0],Interval_arccos([ 0,0])); \n intervc := inapply( cos(2*intervPi/K), K); \n fprintf(f name, `virtual Float Eigenvalue(int K) \{\\n`);\n fprintf(fname, ` st atic INTEGER64 EV[] = \{\\n` );\n # do infinty\n expandedKobFactor6 \+ := subs( c = 1, KobFactor6);\n approxEV := fsolve( expandedKobFactor6 , lambda, lambda=0.5..1);\n if op(2, intervKobFactor6(approxEV-eps, 1 )) > 0 \n or op(1,intervKobFactor6(approxEV+eps, 1)) < 0 then \n \+ ERROR(`fsolve precision failure for infinity`);\n fi;\n fprintf( fn ame, `CONST64(%d),CONST64(%d),CONST64(%d),\\n`, op(1,approxEV-eps), op (1,approxEV+eps),10^(-op(2,approxEV)));\n # skip 1,2\n fprintf( fnam e, `CONST64(0),CONST64(0),CONST64(0), CONST64(0),CONST64(0),CONST64(0) ,\\n`);\n for K from 3 to N do \n if (K - 3) mod 100 = 0 then pr int(K);\n fi; \n cK := intervc(K);\n expandedKobFactor6 := subs( c = cos(2*Pi/K), KobFactor6);\n approxEV := fsolve( expanded KobFactor6, lambda, lambda=0.25..1); \n # check that the precision \+ is at least eps \n if op(2, intervKobFactor6(approxEV-eps, cK)) > \+ 0 \n or op(1, intervKobFactor6(approxEV+eps, cK)) < 0 then \n \+ ERROR(`fsolve precision failure for K =`, K);\n fi;\n # writ ing out f.p. numbers means relying on the parser of \n # the compil er, which will probably do rounding incorrectly when converting to the \n # binary format; instead, we write out integers and assume that the interval numbers are\n # intialized correctly from integers, i .e. conversion is done by the processor in the correct \n # roundin g mode;\n fprintf(fname, `CONST64(%d),CONST64(%d),CONST64(%d),\\n`, op(1,eval(approxEV-eps)), op(1,eval(approxEV+eps)),10^(-op(2,approxEV )));\n od;\n fprintf(fname, `CONST64(0)\};\\n return Float(EV[3*K ],EV[3*K+1])/Float(EV[3*K+2]) ;\\n\}\\n\\n`);\n NULL;\nend: " "6#>% 3ComputeEigenvaluesGR6%'%\"NG%(integerG'%$epsG%(numericG'%&fnameG%'str ingG7-%\"KG%1intervKobFactor6G%)intervPiG%(intervcG%3expandedKobFactor 6G%)approxEVG%\"rG%3deflatedKobFactor6G%\"iG%)marTableG%#cKG6\"F%' DigitsG\"#:>F2-%(inapplyG6%%+KobFactor6G%'lambdaG%\"cG>F3-%/Interval_t imesG6$7$$\"#?!\"\"$\"#?!\"\"-%0Interval_arccosG6#7$\"\"!FW>F4-FC6$-%$ cosG6#*(\"\"#\"\"\"F3FjnF1!\"\"F1-%(fprintfG6$F.%Cvirtual~Float~Eigenv alue(int~K)~|fr|+G-F]o6$F.%<~static~INTEGER64~EV[]~=~|fr|+G>F5-%%subsG 6$/FG\"\"\"FE>F6-%'fsolveG6%F5FF/FF;$\"\"&!\"\"\"\"\"@$52FW-%#opG6$\" \"#-F26$,&F6FjnF+F[o\"\"\"2-Fgp6$\"\"\"-F26$,&F6FjnF+Fjn\"\"\"FW-%&ERR ORG6#%Ffsolve~precision~failure~for~infinityG-F]o6'F.%FCONST64(%d),CON ST64(%d),CONST64(%d),|+G-Fgp6$\"\"\",&F6FjnF+F[o-Fgp6$\"\"\",&F6FjnF+F jn)\"#5,$-Fgp6$\"\"#F6F[o-F]o6$F.%_oCONST64(0),CONST64(0),CONST64(0),~ CONST64(0),CONST64(0),CONST64(0),|+G?(F1\"\"$FjnF(%%trueGC(@$/-%$modG6 $,&F1Fjn\"\"$F[o\"$+\"FW-%&printG6#F1>F;-F46#F1>F5-Feo6$/FG-Ffn6#*(\" \"#Fjn%#PiGFjnF1F[oFE>F6-F[p6%F5FF/FF;$\"#D!\"#\"\"\"@$52FW-Fgp6$\"\"# -F26$,&F6FjnF+F[oF;2-Fgp6$\"\"\"-F26$,&F6FjnF+FjnF;FW-Fgq6$%Afsolve~pr ecision~failure~for~K~=GF1-F]o6'F.F\\r-Fgp6$\"\"\"-%%evalG6#,&F6FjnF+F [o-Fgp6$\"\"\"-F[w6#,&F6FjnF+Fjn)\"#5,$-Fgp6$\"\"#F6F[o-F]o6$F.%_oCONS T64(0)|hr;|+~return~Float(EV[3*K],EV[3*K+1])/Float(EV[3*K+2])~;|+|hr|+ |+G%%NULLGF<6#FEF<" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 72 "The deriv ative of the largest eigenvalue with respect to c at infinity. " }} {PARA 0 "" 0 "" {TEXT -1 190 "To establish C1-continuity for all valen ces, we need to analyze behavior of the magnitude of the largest eigen value as the function of the valence, as the valence increases to i nfinity ( " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 53 " approaches 1). We estimate a constant B, such that " }{XPPEDIT 18 0 "abs(lambda - l ambda[infinity]) < B*abs(c-1)" "6#2-%$absG6#,&%'lambdaG\"\"\"&F(6#%)in finityG!\"\"*&%\"BGF)-F%6#,&%\"cGF)\"\"\"F-F)" }{TEXT -1 77 ", suffic iently close to 1. This constant can be taken to be the maximum of " }{XPPEDIT 18 0 "abs(diff(lambda,c))" "6#-%$absG6#-%%diffG6$%'lambdaG% \"cG" }{TEXT -1 36 ", or, equivalently, as maximum of " }{XPPEDIT 18 0 "abs( diff(c,lambda))^(-1)" "6#)-%$absG6#-%%diffG6$%\"cG%'lambdaG ,$\"\"\"!\"\"" }{TEXT -1 52 "; as the characteristic polynomial is qu adratic in " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 50 ", the latter i s relatively easy to compute. Once " }{XPPEDIT 18 0 "B" "6#%\"BG" } {TEXT -1 36 " is known, we can estimate the size " }{XPPEDIT 18 0 "eps ilon" "6#%(epsilonG" }{TEXT -1 22 " of the interval for " }{XPPEDIT 18 0 "lambda" "6#%'lambdaG" }{TEXT -1 7 " near " }{XPPEDIT 18 0 "lamb da[infinity]" "6#&%'lambdaG6#%)infinityG" }{TEXT -1 134 " , such that \+ if the characteristic map is injective and regular for all these value s, it is sufficient to establish C1-continuity for " }{XPPEDIT 18 0 "K > K[0]" "6#2&%\"KG6#\"\"!F%" }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "1 - cos(2*Pi/K[0]) > epsilon/B" "6#2*&%(epsilonG\"\"\"%\"BG!\"\",&\"\" \"F&-%$cosG6#*(\"\"#F&%#PiGF&&%\"KG6#\"\"!F(F(" }{TEXT -1 1 "." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "intervcdiff := inapply( clam bdadiff, lambda):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 105 "Evaluate fo r all lambda in the range 0.7..1; step 0.001 gives reasonable bounds; \+ this may take some time." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 161 "cdiffinterv := []; \nfor i from 0 to 299 do \n cdiffinterv := \+ Interval_union(cdiffinterv, intervcdiff([0.7+i*0.001, 0.7+(i+1)*0.001 ])); \nod: eval(cdiffinterv);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,cd iffintervG7\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$$\"+K43%)**!\"*$\"+ ?T]j:!\")" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "B := op(2, Int erval_reciprocal( op(1, cdiffinterv)));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"BG$\"+YWf,5!#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "For e xample, if we use the interval of size " }{XPPEDIT 18 0 "1e-6" "6#$\" \"\"!\"'" }{TEXT -1 6 " for " }{XPPEDIT 18 0 "lambda[infinity]" "6#&% 'lambdaG6#%)infinityG" }{TEXT -1 65 ", we have to consider all valen ces up to the valence for which " }{XPPEDIT 18 0 "B*abs( cos(2*Pi/K) - 1) < 1e-6" "6#2*&%\"BG\"\"\"-%$absG6#,&-%$cosG6#*(\"\"#F&%#PiGF&%\"KG! \"\"F&\"\"\"F2F&$\"\"\"!\"'" }{TEXT -1 39 " , which turns out to be a pprox. 1450." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 140 "Interval_times( B, Interval_add( Interval_cos( Interval_times( Interval_times( 2, 2*Inte rval_arccos(0.0)), Interval_reciprocal(1450))),-1));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$$!+e%pMS*!#;$!+N\"pKS*F&" }}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 12 "Eigenvectors" }}{PARA 0 "" 0 "" {TEXT -1 229 "Finall y, we derive the expressions for the complex eigenvector of the large st eigenvalue of the subdivision matrix. We use the fact that the larg est eigenvalue has multiplicity 1 and is a real eigenvalue of the fi rst subblock " }{XPPEDIT 18 0 "B[0,0]" "6#&%\"BG6$\"\"!F&" }{TEXT -1 29 " of the subdivision matrix. " }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 30 " First part of the eigenvector" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 116 "KobBlock00 := subs( s^2 = 1-c^2, map( evalc, subs( \+ \{ d = 0, omega = c + I*s\}, submatrix( KobExpanded, 1..6,1..6))));" } }{PARA 11 "" 1 "" {XPPMATH 20 "6#>%+KobBlock00G-%'MATRIXG6#7(7(,&#\"\" &\"\")\"\"\"%\"cG#F.F-,(*$F/\"\"##!\"\"\"#sF/F4*(%\"IGF.F/F.%\"sGF.#F. F6#F5\"#;\"\"!F=F=7(,(#\"#X\"$G\"F.F/F@*&F8F.F9F.F@,&#F,F " 0 "" {MPLTEXT 1 0 113 "collect( simplify( subs( \{ s^3 = s*(1-c^2), s^2 = 1 - c^2\}, charpoly(KobBlock00,lambda))), lambda) - KobFactor6;" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "redBlock00 := submatrix ( evalm( KobBlock00 - lambda \+ * &*()), 2..6, 1..6): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 117 "v0 := map( simplify, subs( \{ s^3 = (1-c^2)*s, s^2 = 1-c^2, _t[1] = 1 \}, linsolve( redBlock00, vector([0,0,0,0,0]) )));" }}{PARA 12 "" 1 " " {XPPMATH 20 "6#>%#v0G-%'VECTORG6#7(%'lambdaG,$*(F)\"\"\",:*&%\"cG\" \"#F)F,F0!\"\"F,F/F1*$F)F0\"$%Q*&%\"IGF,%\"sGF,F1*&F/F,F)F0F3*&F/F,F) \"\"$!%gD*$F)F9F:**F5F,F)F,F/F,F6F,F0*(F5F,F)F9F6F,F:*(F5F,F)F0F6F,F3* &F/F,F)F,F0F,,4F)\"#!)F.\"\"%F?!#;F2!%7@F8!%?^F1F,F7\"$w&*$F)FB!&Ob'F; \"&![?F1\"\"*F,,$**F5F,F)F,,8F6!\"**&F5F,F/F,FJ*&F5F,F/F0!\"#*&F6F,F/F ,F0*(F5F,F)F,F/F,!$G#*(F5F,F/F,F)F0\"%W8*&F6F,F)F,\"$S#*&F5F,F)F,!$W\" *&F5F,F)F9\"%'4%F5\"#5*&F)F0F6F,!%W8F,F@F1FJ,$*&F-F,F@F1FJ,$**F5F,F)F, ,8FO!\")*(F5F,F/F,F)F9!%'4%F6FgnF5FNFWFZFS\"$W\"*(F6F,F)F,F/F,\"#7*&F) F9F6F,FfnFYFX*&F5F,F)F0Fin*(F5F,F)F,F/F0!#7F,F@F1FN" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 40 "Verify agreement with the regular case: " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "subs( \{ s = 1, c = 0, lambda = 1/2 \}, eval(v0) );" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'VECTORG6#7(#\"\" \"\"\"#,&F'F(%\"IGF'F(,&F(F(F+F',&F(F(F+F(,&F'F(F+F(" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 47 " Second part, separate real and imaginary parts" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "Now we compute the secon d part of the vector:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "KobBlock10 := submatrix ( KobExpanded, 7..12, 1..6);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%+KobBlock10G-%'MATRIXG6#7(7(#\"\"*\"#;\"\"!F*F-F-F-7( ,(-%*conjugateG6#%&omegaG#\"\"\"\"$c##\"#\")F6F5F3#!\"*F6,&F0F9F7F5,&F 3F4F7F5F7F9,&F9F5F0F97(,$F3#!\"\"F,F*F-F*F-F-7(,&F9F5F3F9F7FCF7F7F77(F @F*F-F-F-F*7(,(F9F5F3F7*$F3\"\"#F4,&F7F5F3F9,&F4F5F3F7FCF9F7" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "KobBlock11lambda := submatri x ( evalm( KobExpanded - lambda * &*()), 7..12, 7..12);" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%1KobBlock11lambdaG-%'MATRIXG6#7(7(,&#!\"\"\"#; \"\"\"%'lambdaGF,\"\"!F0F0F0F07(#!\"*\"$c#,&F2F.F/F,#F.F4F0F0,$-%*conj ugateG6#%&omegaGF67(F0F+,$F/F,F0F0F07(,&F6F.F;F6F2F2,&F6F.F/F,F2F27(F0 F0F0F0F=F+7(,$F;F2,$F;F6F0F0F6F5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 188 "v1 := map(simplify, subs( \{ s^4 = (1-c^2)^2, s^2 = \+ 1-c^2, s^3 = (1-c^2)*s\}, map( simplify, map( evalc, subs( omega = I* s + c, evalm( - inverse(KobBlock11lambda) &* KobBlock10 &* v0)))))):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 25 "Put together the vector: " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "KobEigenvect := vector( [seq( v0[i] , i = 1..6), seq( v1[i], i = 1..6) ]):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 39 "Verify agreement with the regular case:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 85 "subs( \{ s = 1, c = 0, lambda = 1/2\}, map( sim plify, evalm( eval(KobEigenvect)) ) );" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'VECTORG6#7.#\"\"\"\"\"#,&F'F(%\"IGF'F(,&F(F(F+F',&F(F(F+F(,&F 'F(F+F(#\"\"$F),&F/F(F+F',&F/F(F+F(,&F/F(F+F/,&F(F(F+F/,&F'F(F+F/" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "Separate real and complex parts" } {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "KobEige nvectRe := map( evalc, map( Re, KobEigenvect)):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "In addition, scale imaginary part by by 1/s " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 87 "KobEigenvectIm := map( simplify, e valm( (1/s) * map( evalc, map( Im, KobEigenvect)))):" }}}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 15 "Code generation" }}{PARA 0 "" 0 "" {TEXT -1 53 "A class with three memeber functions are generated:\n " }{TEXT 256 23 "Float Eigenvalue(int K)" }{TEXT -1 28 " computes the eigenvalu es, \n" }{TEXT 257 58 "void EigenvectorReal(Float c, Float lambda, Flo at EvRe[]) " }{TEXT -1 68 "initializes an array for the real part of t he complex eigenvector, " }{TEXT 258 62 "void EigenvectorImaginary(Fl oat c, Float lambda, Float EvIm[])" }{TEXT -1 108 " initializes the ar ray for the complex part. \nMemory for arrays should be allocated by t he calling function." }}{PARA 0 "" 0 "" {TEXT -1 267 "The output is wr itten to a file; if the name is `default`, it is written to the stand ard output (warning: for some reason, writing to standard output is te rribly slow; writing to a file and then looking at it in an editor is much more efficient. All functions use " }{TEXT 260 5 "Float" } {TEXT -1 158 " as the name of the class for the interval numbers.\nIt is assumed to have explicit casts from 64-bit integers, standard ar ithmetics operations, and macros " }{TEXT 261 3 "FR " }{TEXT -1 6 " an d " }{TEXT 262 4 "Fdiv" }{TEXT -1 9 ", (see " }{MPLTEXT 0 21 14 "Co nvertToFloat" }{TEXT -1 15 " for details).\n" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "OutputFile := `kobbelt.cpp`:" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 280 "Generate the header and the constructor for the c lass; specify the name of the class, the inner and outer diameters of \+ the layer for the char. map, the radius of the neighborhood on wehich \+ the eigenvector is defined, and the name of the corresponding scheme\n on the regular complex" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 61 "MakeClassHeader( OutputFile, `Kobbelt`, 2,4,3, TensorFourPt):" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "ComputeEigenvalues(3000, 1 e-10, OutputFile):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 72 "This functi on is needed for computing the eigenvectors \"near infinity\"." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 233 "fprintf( OutputFile, `virtual void EigenvalueRange(Float c, Float& lambdamin, Float& lambdamax) \{\\n`): \nfprintf( OutputFile, `lambdamax = Eigenvalue(0); lambdamin = exactfl oor(FR(CONST64(%d), CONST64(%d))*c);\\n \}`, op(1,B), op(2,B)):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "GenerateEigenvectorCode(KobE igenvect, OutputFile);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "End of \+ the class " }{MPLTEXT 1 0 26 "fprintf(OutputFile, `\};`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "fclose(OutputFile):" }}}}}{MARK "0 \+ 0" 0 }{VIEWOPTS 1 1 0 3 4 1802 }