add u11,u12,u13

Author: heltonmc

src/besselk.jl

@@ -75,7 +75,14 @@
 #    u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 
 #         - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 
 #         - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10
-#   
+#    u11 = p^11 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10
+#        + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11
+#    u12 = p^12 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14
+#        + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12
+#    u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20
+#        - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10
+#        + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13
+#        
 #    For large orders these formulas will converge much faster than using upward recurrence.
 #    If using up to the fifth U polynomial, it requires evaluation of a 15 degree polynomial.
 #    For tenth U polynomial, it requires evaluation of a 30 degree polynomial.
@@ -234,10 +241,19 @@ function besselk_large_orders(v, x::T) where T <: AbstractFloat
     u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8
     u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9
     u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10
-    out1 = ((u10 + u9) + u8) + u7
+    u11 = p^11 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10
+        + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11
+    u12 = p^12 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14
+        + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12
+    u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20
+        - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10
+        + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13
+        
+        
+
+    out1 = ((u13 + u12 + u11 + u10 + u9) + u8) + u7
     out2 = ((out1 + u6 + u5) + u4) + u3
     out = ((out2 + u2) + u1) + u0
 
     return coef*out
  end
-